************************************************************************ * * * E R R A T A * * * * PLEASE NOTE THE FOLLOWING CHANGES TO THE CPC ARTICLE WHICH DESCRIBES * * CONTWV, I.E. CPC, VOL. 66 (1991), 392-402. * * * * 1) THE LINE NUMBERS THAT WERE REFERRED TO IN THE ADAPTATION * * INSTRUCTIONS ARE IN ERROR. THE CORRECT LINE NUMBERS ARE: * * 48, 56, 103, INSTEAD OF 74, 82, 129, RESPECTIVELY. * * * * 2) THE TEST RUN INPUT REFERRED TO OPTIONS WHICH WERE NOT PART OF THE * * ORIGINAL GRANT, ET AL (1980) CODE. THE CORRECT TEST RUN INPUT * * SHOULD THEREFORE HAVE THE FOLLOWING CHANGES: * * A) THE MCDF CARD SHOULD HAVE OPTION 11 INSTEAD OF OPTION 22 * * B) THE CARD IMMEDIATELY FOLLOWING THE MCDF CARD SHOULD BE: * * 8 9 10 * * INSTEAD OF: * * 8 0 0 13 0 * * C) THE NUCLEAR WEIGHT SHOULD BE GIVEN ON A CARD AS: * * NUC 1 200.59 * * INSTEAD OF: * * 200.59 * * D) THE LOAD CARD SHOULD BE DELETED. * * E) THE EAL CARD SHOULD BE: * * EAL 1 * * * * THE AUTHORS APOLOGIZE FOR ANY INCONVENIENCE THAT THESE CORRECTIONS * * MAY HAVE CAUSED THE USER AND SUGGEST THAT ANY QUESTIONS CONCERNING * * THESE CHANGES MAY BE DIRECTED TO: * * * * WARREN F. PERGER * * PHYSICS DEPARTMENT * * MICHIGAN TECH UNIVERSITY * * HOUGHTON, MI 49931-1295 * * PHONE: 906-487-2855 * * INTERNET ADDRESS: wfp@mtu.edu * * * * E N D O F E R R A T A * ************************************************************************ C***************************************************************** C***************************************************************** C***************************************************************** C C SECTION 2 D A T A I N C C 14 ROUTINES C C 1. DATAIN C 2. DATR C 3. DATANG C 4. DATNR C 5. MAN3 C 6. COUP C 7. MANOUT C 8. DATSCF C 9. CFOUT C 10. CARDIN C 11. BLOCK DATA (TERMS) C 12. BLOCK DATA (CONS) C 13. TIMER C 14. CALEN C C***************************************************************** SUBROUTINE DATAIN IMPLICIT REAL*8(A-H,O-Z) C DIMENSION RT(33),IT(33) C C THIS SUBROUTINE CONTROLS THE READING C OF THE INPUT DATA C THE FREE FORMAT CARD READING ROUTINE IS USED C C SUBROUTINES CALLED : CALEN,CARDIN,DATSCF,DATR,DATNR C COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN COMMON/TITL/IHED(20),ITIME(2),IDATE(2) COMMON/CARD/LABEL COMMON/ORB4/NW,NCF,NH(25),NP(25),NAK(25),IQ(25,30) COMMON/DEF2/NCFX,NMAX,NFMAX,NMCPX,NECX,NSCFX COMMON/KNTR/ITC(20) COMMON/INFORM/IRD,IPD COMMON/MCTR/JTC(20) COMMON/MCTA/KA,IOPAR COMMON/MCTB/IO1,IO2,IO3,IO4,IO5,IO6,IO7 COMMON/STAT/JQS(3,25,30),JCUP(8,30),ICHOP(25,30),INCOR,IAB CBBBB COMMON/OPT3/KTC(20) CBBBB COMMON/ATOM/ATW CBBBB COMMON/WFAC/WFACT COMMON/GRIDC/RGRIDC(2000),ECONT,HCONT,MHALF,NCONT,NUMWR C DATA LABEL1,LABEL2,LABEL3/4HMCP ,4HMCDF,4HEND / DATA LABEL4/4HCONT/ DATA LABEL5/4HMCT / DATA LABEL6,LABEL7/4HMCBP,4HBENA/ C C SET THE DATA SET (DS) NUMBER FOR C (1) CARD READER (IRD) C (2) LINEPRINTER (IPD) C IRD=5 IPD=6 C C ************************************************ C ***** READ A TITLE AND ORBITAL DATA ************ C ************************************************ C C CARD 1 TITLE (COL. 1-80) C C THIS READ STATEMENT USES AN END PARAMETER TO FIND C THE END OF DATA C REPLACE IF THIS IS NOT AVAILABLE C C THIS IS THE ONLY CARD WHICH IS NOT READ BY CARDIN C ECONT = -1.0D0 LABEL=1 READ (IRD,300,END=25) IHED WRITE (IPD,301) LABEL,IHED CALL CALEN C C SET MAXIMUM DIMENSIONS FOR ARRAYS USED IN PROGRAMS C C NCFX - MAX. NO. OF CONFIGURATIONS C NWX - MAX. NO. OF ORBITALS C ND - MAX. NO. OF NUMBERS READ BY CARDIN C NCFX=30 NWX=25 ND=NCFX+3 C C SET THE OPTIONS TO 0 C C IF OPTION J IS TO BE SET THEN AT THE APPROPRIATE PLACE C INPUT THE NUMBER J . THIS WILL SET JTC(J)=1 FOR EXAMPLE . C IF JTC(J)=0 THE OPTION IS NOT SET . C C ITC CONTAINS OPTIONS FOR THE MCDF PROGRAM C JTC(J) ,J=1,9 CONTAINS OPTIONS FOR THE ANGULAR PROGRAMS C JTC(J) ,J=10,20 CONTROL WHICH PROGRAMS ARE CALLED C KTC CONTAINS OPTIONS FOR BENA PROGRAM C DO 1 I=1,20 JTC(I)=0 ITC(I)=0 CBBBB KTC(I)=0 1 CONTINUE C C THERE ARE TWO FORMS OF INPUT FOR THE ORBITAL DATA C C (1) THE CONFIGURATIONS CAN BE DEFINED AS DESCRIBED IN THE C CPC WRITE UP - SEE SUBROUTINE DATR FOR A DESCRIPTION C THIS FORM REQUIRES THE DEFINITION OF EACH RELATIVISTIC C CONFIGURATION C (2) ALTERNATIVELY NON-RELATIVISTIC CONFIGURATIONS CAN BE C DEFINED C - SEE SUBROUTINE DATNR FOR A DESCRIPTION C RELATIVISTIC CONFIGURATIONS ARE GENERATED FROM THESE C C CARD 2 C C READ TWO OR THREE INTEGERS C C NMAN - NUMBER OF CONFIGURATIONS C NWM - NUMBER OF ORBITALS C IOP - >0 IF DATNR IS TO BE CALLED C =0 IF DATR IS TO BE CALLED C C IF ONLY TWO NUMBERS ARE INPUT THEN IOP=0 BY DEFAULT C AND DATR FOR RELATIVISTIC INPUT IS CALLED C CALL CARDIN (0,NA,NB,NN,RT,IT,ND) IF (NN.NE.2.AND.NN.NE.3) GO TO 24 NMAN=IT(1) NWM=IT(2) IOP=0 IF (NN.EQ.3) IOP=IT(3) IF (IOP.EQ.0) GO TO 2 CALL DATNR (NMAN,NWM,NWX,ND) GO TO 3 2 CALL DATR (NMAN,NWM,NWX,ND) 3 CONTINUE C C *********** C *** MCP *** C *********** C C FORMAT : LABEL (COL.1-4) , 1 - 3 NUMBERS C C LABEL = 'MCP ' - IF THIS LABEL IS FOUND THEN A CALL IS C MADE TO THE ANGULAR PACKAGE TO OBTAIN C MCP COEFFICIENTS REQUIRED BY THE MCDF C PROGRAM C C THE LABEL IS FOLLOWED BY : C (1) DS NUMBER FOR OUTPUT OF MCP COEFFICIENTS TO FILE . C CARE SHOULD BE TAKEN THAT THIS NUMBER AGREES WITH C THE DS NUMBER USED TO READ THE MCP COEFFICIENTS C IN THE MCDF PROGRAM IF THE PROGRAMS ARE RUN C TOGETHER . C SET TO ZERO IF NO DISC/TAPE FILE REQUIRED. C (2) TWO POSSIBLE OPTIONS C OPTION 1 - INCLUDE CORE/PEEL COEFFICIENTS . C OPTION 2 - DO NOT INCLUDE I(A,B) COEFFICIENTS C CALL CARDIN (2,NA,NB,NN,RT,IT,ND) IF (NB.NE.LABEL1) GO TO 8 IF (NN.EQ.0.OR.NN.GT.3) GO TO 24 JTC(12)=JTC(11) IO1=IT(1) DO 5 I=1,9 5 JTC(I)=0 IF (NN.EQ.1) GO TO 7 DO 6 I=2,NN J=IT(I) IF (J.LT.1.OR.J.GT.2) GO TO 24 6 JTC(J)=1 7 INCOR=JTC(1) IAB=JTC(2) C C *********** C *** MCDF ** C *********** C C FORMAT : LABEL (COL.1-4) , 0 - 20 NUMBERS C C LABEL = 'MCDF' - IF THIS LABEL IS FOUND THEN A CALL IS C MADE TO DATSCF WHICH READS C THE DATA REQUIRED BY THE MCDF PROGRAM C C THE LABEL IS FOLLOWED BY A LIST OF OPTIONS C UP TO 20 OPTIONS CAN BE SPECIFIED C SEE DATSCF FOR A DESCRIPTION C CALL CARDIN (2,NA,NB,NN,RT,IT,ND) 8 IF (NB.NE.LABEL2) GO TO 11 IF (NN.GT.20) GO TO 24 JTC(13)=1 IF (NN.EQ.0) GO TO 10 DO 9 I=1,NN J=IT(I) IF (J.LT.1.OR.J.GT.20) GO TO 24 9 ITC(J)=1 IF (ITC(20).EQ.1) JTC(13)=0 10 CALL DATSCF (ND) * * * * C ************ C *** CONT *** C ************ * * * FORMAT : LABEL (COL.1-4) , 0 - 20 NUMBERS * * * * LABEL = 'CONT' - IF THIS LABEL IS FOUND THEN CONTINUUM ORBITALS * * ARE CALCULATED, FIXING THE CORE. * * * * THE LABEL IS FOLLOWED BY THE LOCATION OF WHERE THE GRID SWITCHES * * OVER FROM EXPONENTIAL TO LINEAR (DEFAULT=200), FOLLOWED BY THE * * FILE NUMBER TO WRITE THE CONTINUUM ORBITAL TO (DEFAULT=15). * * * * CALL CARDIN (2,NA,NB,NN,RT,IT,ND) 40 IF (NB .NE. LABEL4) GO TO 11 IF (RT(1) .LT. 0.0D0) THEN WRITE (IPD,*) ' ERROR IN DATAIN--CONTINUUM E MUST BE > 0' STOP END IF ECONT = RT(1) IF (IT(2) .EQ. 0) THEN MHALF = 200 ELSE MHALF = IT(2) END IF IF (IT(3) .EQ. 0) THEN NUMWR = 15 ELSE NUMWR = IT(3) END IF C C *********** C *** MCT *** C *********** C C FORMAT : LABEL (COL.1-4) , 2 - 3 NUMBERS C C LABEL = 'MCT ' - IF THIS LABEL IS FOUND THEN A CALL IS C MADE TO THE ANGULAR PACKAGE TO OBTAIN C MCT COEFFICIENTS . C C FOLLOWING THE LABEL : C (1) DS NUMBER FOR OUTPUT OF MCT COEFFICIENTS TO FILE C SET TO ZERO IF NO DISC/TAPE FILE REQUIRED. C IF RUNNING MCT AND OSCL TOGETHER THEN ENSURE C THAT THE DS NUMBERS ARE COMPATIBLE C (2) KA - ORDER OF THE TENSOR OPERATOR C (3) IOPAR - PARITY OF THE TENSOR OPERATOR C C IOPAR= 1 FOR EVEN PARITY C =-1 FOR ODD PARITY C = 0 FOR NO PARITY CHECK C C IF ONLY 2 NUMBERS ARE INPUT THEN THE DEFAULT IS IOPAR=0 C CALL CARDIN (2,NA,NB,NN,RT,IT,ND) 11 IF (NB.NE.LABEL5) GO TO 19 IF (NN.LT.2.OR.NN.GT.3) GO TO 24 JTC(14)=JTC(11) IO2=IT(1) KA=IT(2) IOPAR=0 IF (NN.EQ.3) IOPAR=IT(3) C C ************ C *** MCBP *** C ************ C C FORMAT : LABEL (COL.1-4) , 1 NUMBER C C LABEL = 'MCBP' - IF THIS LABEL IS FOUND THEN A CALL IS C MADE TO THE ANGULAR PACKAGE TO OBTAIN C MCBP COEFFICIENTS REQUIRED BY THE BENA C PROGRAM. C C THE LABEL IS FOLLOWED BY THE C DS NUMBER FOR OUTPUT OF MCBP FILE. C SET TO ZERO IF NO DISC/TAPE FILE REQUIRED. C CARE SHOULD BE TAKEN THAT THIS NUMBER AGREES WITH C THE DS NUMBER USED TO READ THE MCBP COEFFICIENTS IN C THE BENA PROGRAM IF THE PROGRAMS ARE RUN TOGETHER. C CALL CARDIN (2,NA,NB,NN,RT,IT,ND) 19 IF (NB.NE.LABEL6) GO TO 20 CBBBB IF (NN.NE.1) GO TO 24 CBBBB JTC(15) = JTC(11) CBBBB IO3 = IT(1) C C ************ C *** BENA *** C ************ C C FORMAT : LABEL (COL.1-4) , 0 - 20 NUMBERS. C C LABEL = 'BENA' - IF THIS LABEL IS FOUND THEN A CALL IS MADE C TO THE BREIT ENERGY ANALYSIS PACKAGE, BENA, C TO OBTAIN THE BREIT AND QED CORRECTIONS TO C THE ENERGY LEVELS. C C THE LABEL IS FOLLOWED BY A LIST OF OPTIONS C UP TO 20 OPTIONS CAN BE SPECIFIED C C OPTION MEANING C C 1 PRINT OUT COULOMB MATRIX C 2 PRINT OUT EIGENVALUES OF COULOMB MATRIX C 3 PRINT OUT EIGENVECTORS OF COULOMB MATRIX C 4 PRINT OUT BREIT MATRIX C 5 PRINT OUT CONTRIBUTION OF BREIT MATRIX ELEMENTS TO C EIGENVALUES C 6 PRINT OUT EIGENVECTORS OF COULOMB + BREIT MATRIX C 8 PRINT OUT EACH CONTRIBUTION TO THE BREIT MATRIX ELEMENTS. C 9 PRINT OUT EACH CONTRIBUTION TO THE ENERGY LEVELS C RESULTING FROM QED CORRECTIONS. C 10 SUPPRESS THE FINAL SUMMARY OF ENERGY LEVEL CONTRIBUTIONS. C 11 GIVE INFORMATION ON BESSEL FUNCTION EVALUATION AND THE C SOURCE OF THE FUNCTION. C 20 SUPPRESS THE CALL TO THE BENA PACKAGE C CALL CARDIN (2,NA,NB,NN,RT,IT,ND) 20 IF (NB.NE.LABEL7) GO TO 23 CBBBB IF (NN.GT.20) GO TO 24 CBBBB JTC(16) = 1 CBBBB IF (NN .EQ. 0) GO TO 22 CBBBB DO 21 I=1,NN CBBBB J=IT(I) CBBBB IF (J.LT.1 .OR. J.GT.20) GO TO 24 CBBBB 21 KTC(J) = 1 CBBBB IF (KTC(20).EQ.1) JTC(16) = 0 C C *********************************** C **** INPUT DATA FOR BREIT **** C **** ENERGY ANALYSIS PACKAGE **** C *********************************** C C FORMAT : 2 - 4 NUMBERS. C C (1) DS NUMBER OF FILE CONTAINING MCBP COEFFICIENTS. C (2) DS NUMBER OF FILE CONTAINING MCDF DUMP. C (3) ATOMIC WEIGHT OF ATOM TO BE USED TO OBTAIN CORRECTIONS C TO THE RYDBERG. IF NO NUMBER OR ZERO IS GIVEN THEN AN C INFINITE ATOMIC MASS IS ASSUMED. C (4) FACTOR MULTIPLYING FREQUENCY TO ENABLE THE LOW C FREQUENCY (LONG WAVELENGTH) CASE TO BE OBTAINED. IF NO C NUMBER OR ZERO IS GIVEN THEN A FACTOR OF 1.0 IS ASSUMED. C (NOTE: A FACTOR OF 1.E-3 IS USUALLY SUFFICIENT TO OBTAIN C THE LOW FREQUENCY LIMIT) C CBBBB 22 CALL CARDIN(0,NA,NB,NN,RT,IT,ND) CBBBB IF (NN.LT.2 .OR. NN.GT.4) GO TO 24 CBBBB IO4=IT(1) CBBBB IO5=IT(2) CBBBB ATW=ZERO CBBBB WFACT=ONE CBBBB IF (NN.GE.3) ATW=RT(3) CBBBB IF (NN.GE.4 .AND. RT(4).GT.ZERO) WFACT=RT(4) C C *********** C *** END *** C *********** C C FORMAT : LABEL (COL.1-4) C C LABEL='END ' - DENOTES THE END OF DATA FOR THIS PROBLEM C C CALL CARDIN(2,NA,NB,NN,RT,IT,ND) 23 IF (NB .EQ. LABEL3) RETURN WRITE (IPD,304) STOP 5 C 24 WRITE (IPD,302) STOP 5 25 WRITE(IPD,303) STOP C 300 FORMAT (20A4) 301 FORMAT (1H1,10(1H+),8H DATAIN ,10(1H+)//11H CARD INPUT// X 1X,I3,2X,1H=,2X,20A4) 302 FORMAT (1H0,24HERROR ON ABOVE DATA CARD) 303 FORMAT (1H0,29H** E N D O F D A T A **) 304 FORMAT (1H0,18HEXPECT "END " CARD) END C***************************************************************** ************************************************************************ * * SUBROUTINE SCF (ERAL,MX) * * * THIS SUBROUTINE CONTROLS THE MAIN SEQUENCE OF THE * * CALCULATION TO PRODUCE A COMPLETE SET OF SELF-CONSISTENT * * WAVE FUNCTIONS. * * * * ERAL IS THE MAXIMUM DISCREPANCY BETWEEN ESTIMATED AND * * CALCULATED WAVE FUNCTIONS. * * M: IS AN ERROR INDICATOR. * * * * SUBROUTINES CALLED : TIMER, LAGR, YPOT, XPOT, EIGEN, SPLINE, SPLIN * * SPLINT, XPOTC, HPOT, SOLV, RINT, DUMP, SOLVC, * * NORMC, MAXMIN. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * DIMENSION ASMT(20),ESMT(20) DIMENSION ER(25) * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /DEF1/Z,RNT,H,EPH,ACCY,C,N,NWA : /DEF2/NCFX,NMAX,NFMAX,NMCPX,NECX,NSCFX : /ORB1/E(25),GAMA(25),XAM(25) : /ORB4/NW,NCF,NH(25),NP(25),NAK(25),IQ(25,30) : /WAVE/PF(350,25),QF(350,25) : /INT2/P(350),Q(350),PC(350),QC(350) : /EXCO/PZ(25),QZ(25) : /KNTR/ITC(20) : /FIXD/JFIX(25) : /INFORM/IRD,IPD : /CONVG/XCAX : /POTE/YP(350),YQ(350),XTP(350),XTQ(350) COMMON/WAVEC/PFC(2000,25),QFC(2000,25) : /POTEC/YPC(2000),YQC(2000),XTPC(2000),XTQC(2000) : /GRIDC/RGRIDC(2000),ECONT,HCONT,MHALF,NCONT,NUMWR : /TATBC/TAC(2000),TBC(2000) : /GRID/RGRID(350) : /OFFD/ECV(50),IECC(50),NEC : /ORB5/NKL(25),NKJ(25) : /HMAT/EAV,EMT(30,30),CCF(30),UCF(25),SA(30) DOUBLE PRECISION PC1(350),QC1(350),PC2(2000),QC2(2000), : CC(2000,3) * DIMENSION ECV1(50),ECV2(50),ORTH1(50),ORTH2(50) DIMENSION ECVTEM(50,50),ORTEM(50,50),ECV3(50),ORTH3(50) * DATA EPS,PEIGHT,PNINE/1.0E-10,0.8,0.9/ PI = 4.0D0*ATAN(ONE) * * IF ECONT > 0, THEN CALCULATE CONTINUUM WAVEFUNCTIONS ASSUMING A * FIXED CORE. * IF (ECONT .GE. ZERO) THEN * * DETERMINE IF LAGRANGE MULTIPLIERS ARE REQUIRED * NEC = 0 IF (ITC(2) .EQ. 1) GO TO 141 IF (NW .EQ. 1) GO TO 142 NWM = NW-1 DO 139 JA = 1,NWM JP = JA+1 DO 138 JB = JP,NW IF (NAK(JA) .NE. NAK(JB)) GO TO 138 IF (NP(JA) .EQ. NP(JB)) GO TO 138 IF (ABS (UCF(JA)-FLOAT (NKJ(JA))-ONE) .GT. EPS) : GO TO 137 IF (ABS (UCF(JB)-FLOAT (NKJ(JB))-ONE) .GT. EPS) : GO TO 137 GO TO 138 137 IF (JFIX(JA)+JFIX(JB) .EQ. 2) GO TO 138 NEC = NEC+1 IF (NEC .GT. NECX) GO TO 143 IECC(NEC) = JA+NWA*JB ECV(NEC) = ZERO 138 CONTINUE 139 CONTINUE IF (NEC .EQ. 0) GO TO 142 WRITE (IPD,315) DO 140 I = 1,NEC MA = IECC(I)/NWA MB = IECC(I)-NWA*MA WRITE (IPD,316) NP(MB),NH(MB),NP(MA),NH(MA) 140 CONTINUE GO TO 136 141 WRITE (IPD,317) GO TO 136 142 WRITE (IPD,318) GO TO 136 * 143 WRITE (IPD,319) STOP 5 * 136 WINECV = 0.5D0 DO 130 I=1,NEC ECV1(I) = ZERO ECV2(I) = ZERO ECV(I) = WINECV ORTH1(I) = ZERO ORTH2(I) = ONE 130 CONTINUE NCONT = 2000 R = RNT DO 50 I = 1,NMAX RGRID(I) = R R = R*EPH IF (R .GT. 10000) GO TO 51 50 CONTINUE GO TO 52 51 CONTINUE NMAX = I 52 CONTINUE HCONT = (RGRID(MHALF)-RGRID(MHALF-1)) WRITE (IPD,*) ' HCONT=',HCONT WRITE (IPD,*) ' RGRID(MHALF)=',RGRID(MHALF) RGRIDC(1) = RGRID(1) DO 700 I = 2,NCONT IF (I .LE. MHALF) THEN RGRIDC(I) = RGRID(I) ELSE RGRIDC(I) = RGRIDC(I-1) + HCONT END IF IF ((RGRIDC(I) .GT. RGRID(N)) .AND. (N .LT. 350))N=N+1 700 CONTINUE DO 710 I = 1,NW IF (JFIX(I) .EQ. 0) THEN DO 758 I1 = 1,NCONT PFC(I1,I) = ZERO QFC(I1,I) = ZERO 758 CONTINUE GO TO 710 END IF DO 720 I1 = 1,N PC1(I1) = PF(I1,I) QC1(I1) = QF(I1,I) 720 CONTINUE CALL SPLINE (RGRID,PC1,N,CC) CALL SPLIN (RGRID,PC1,N,CC,RGRIDC,PC2,NCONT-1) CALL SPLINE (RGRID,QC1,N,CC) CALL SPLIN (RGRID,QC1,N,CC,RGRIDC,QC2,NCONT-1) DO 730 I1 = 1,NCONT-1 PFC(I1,I) = PC2(I1) QFC(I1,I) = QC2(I1) 730 CONTINUE PFC(NCONT,I) = ZERO QFC(NCONT,I) = ZERO 710 CONTINUE END IF MX = 0 ITR = 4*NW J = 0 * XINC = XCAX*TENTH WRITE (IPD,301) * * INSERT ARTIFICIAL VALUES INTO THE ARRAY ER, SO THAT * THE FIRST OPERATION WILL CONSIST OF ONE STAGE OF ITERATION * FOR EACH FUNCTION IN THEIR NATURAL ORDER. * DO 1 I = 1,NW ER(I) = ZERO 1 CONTINUE WA = FLOAT (NW) DO 2 I = 1,NW IF (JFIX(I) .EQ. 0) ER(I) = 1000.0*WA IF ((JFIX(I) .EQ. 0) .AND. (ECONT .GE. ZERO)) NP(I) = -1 WA = WA-ONE 2 CONTINUE * * JL IS THE SERIAL NUMBER OF THE FUNCTION TREATED IN THE * PREVIOUS ITERATION. IN THE FIRST ITERATION IT IS ZERO. * * DETERMINE WA, THE LARGEST DISCREPANCY, AND J, WHICH * SPECIFIES WHICH FUNCTION SHOWS THIS LARGEST VALUE. * 3 WA = -ONE JL = J DO 5 I = 1,NW WB = ABS (ER(I)) IF (WB-WA) 5,5,4 4 WA = WB J = I 5 CONTINUE * * NORMAL EXIT IF REQUIRED ACCURACY HAS BEEN ACHIEVED. * IF (WA .LE. ERAL) RETURN * * ITR COUNTS ITERATIONS DOWNWARDS FROM THE SPECIFIED MAXIMUM. * 90 ITR = ITR-1 IF (ITR) 47,6,6 * * STORE ATOMIC NUMBER IN ZST. IPY INDICATES WHICH STAGE * THE AUTOMATIC RECOVERY PROCESS HAS REACHED: * IPY = 0 NORMAL COMPUTATION OF WAVE FUNCTION * IPY = 1 AFTER INITIAL FAILURE, ATOMIC NUMBER IS * BEING INCREASED * IPY = 2 SOLUTION HAS BEEN FOUND WITH ARTIFICIALLY * INCREASED ATOMIC NUMBER, WHICH IS NOW * BEING REDUCED. * THIS REDUCTION OF THE ATOMIC NUMBER IN THE THIRD STAGE IS * ACHIEVED IN 10 STAGES, WHICH ARE COUNTED IN IPW, STEPPING * DOWN FROM 10 TO ZERO. * 6 ZST = Z ZMX = Z IPW = 10 IPY = 0 IPWL = 10 XCA = ONE ZFAC = ONE * * ADJUSTMENT OF OFF-DIAGONAL PARAMETERS. * IF (ECONT .LT. ZERO) CALL LAGR (J) * * XCA IS THE FACTOR MULTIPLYING EXCHANGE TERMS. * Z IS THE EFFECTIVE ATOMIC NUMBER. * * TABULATE THE DIRECT AND EXCHANGE POTENTIALS. * 8 CONTINUE IF (ECONT .LT. ZERO) THEN CALL YPOT (J,VN,ZFAC) CALL XPOT (J,XCA,XCP,XCQ) END IF IF ((ECONT .GT. ZERO) .AND. (JFIX(J) .EQ. 0)) THEN CS = C*C E(J) = -ECONT XMOMEN = SQRT(((ECONT+CS)**TWO-CS*CS)/CS) WRITE (IPD,*) ' ENERGY=',-E(J),' K-VALUE=',XMOMEN KOUNT = 0 KOUNT2 = 0 VALUE = ZERO ISTOP = 0 ISTP = 1 ISTP3 = 0 ACCYSTP=5.0D-4 ACCY2=2.0D-4 770 KOUNT = KOUNT + 1 IF (KOUNT .GT. 500) THEN WRITE (IPD,*) ' NUMBER OF ITERATIONS IN SCF EXCEEDED' STOP END IF DO 127 I1=1,NCONT TAC(I1)=PFC(I1,J)*PFC(I1,J)+QFC(I1,J)*QFC(I1,J) 127 CONTINUE CALL SIMP (TAC,H,1,MHALF,VALUE1,0) CALL SIMP (TAC,HCONT,MHALF,NCONT,VALUE2,0) VALUE = VALUE1+VALUE2 IF (KOUNT .EQ. 1) GO TO 743 IF (ABS(VALUE) .LT. 1.0E-20) THEN ISTOP = 1 GO TO 743 END IF IF (ABS((OLDVAL-VALUE)/VALUE) .LT. ACCYSTP) ISTOP = 1 IF (ISTP .EQ. 0) THEN ISTOP = 0 ISTP = 1 END IF AG = PZ(J) 743 OLDVAL = VALUE IF (KOUNT .EQ. 2) THEN CALL NORMC (J,XX) AG = PZ(J)*XX END IF DO 783 I1 = 1,NCONT PC2(I1) = PFC(I1,J) QC2(I1) = QFC(I1,J) 783 CONTINUE * IF (ISTOP .EQ. 0) THEN CALL YPOT (J,VN,ZFAC) * * INTERPOLATE THE DIRECT POTENTIAL * CALL SPLINE (RGRID,YP,N,CC) CALL SPLIN (RGRID,YP,N,CC,RGRIDC,YPC,NCONT-1) YPC(NCONT) = YPC(NCONT-1) DO 773 I1 = 1,NCONT YQC(I1) = YPC(I1) 773 CONTINUE XCA=ONE CALL XPOTC (J,XCA,XCP,XCQ) CALL SOLVC (J,VN,XCA,XCP,XCQ,AG,WA,MF) * DO 784 I1 = 1,NCONT PFC(I1,J) = 0.6D0*PFC(I1,J)+0.4D0*PC2(I1) QFC(I1,J) = 0.6D0*QFC(I1,J)+0.4D0*QC2(I1) 784 CONTINUE GO TO 770 END IF IF ((NEC .EQ. 0) .OR. (ITC(2) .EQ. 1) .OR. (ISTP3 .EQ. 1)) : GO TO 120 KOUNT2 = KOUNT2+1 DO 119 LB = 1,NEC JA = IECC(LB)/NWA JB = IECC(LB)-NWA*JA IF (J-JA) 102,101,102 101 JJ = JB GO TO 104 102 IF (J-JB) 119,103,119 103 JJ = JA 104 CONTINUE ORTH = RINTC (J,JJ,0) ORTEM(LB,KOUNT2) = ORTH ECVTEM(LB,KOUNT2) = ECV(LB) IF (KOUNT2 .EQ. 1) THEN ORTH1(LB) = ORTH ECV1(LB) = ECV(LB) ECV2(LB) = -ECV(LB) ECV(LB) = ECV2(LB) GO TO 119 END IF IF (KOUNT2 .EQ. 2) THEN IF (ORTH1(LB)*ORTH .GT. ZERO) THEN WRITE (6,*) ' SOLN NOT BRACKETED' WRITE (6,*) ' INCREASE SIZE OF WINECV IN SR SCF' STOP END IF ECV(LB) = (ECV1(LB)+ECV2(LB))/TWO ORTH2(LB) = ORTH GO TO 119 END IF IF (ORTH2(LB)*ORTH .GT. ZERO) THEN ORTH2(LB) = ORTH ECV2(LB) = ECV(LB) ECV(LB) = (ECV(LB)+ECV1(LB))/TWO ELSE ORTH1(LB) = ORTH ECV1(LB) = ECV(LB) ECV(LB) = (ECV2(LB)+ECV(LB))/TWO END IF 119 CONTINUE ISTOP = 0 ISTP = 0 NPTS = 8 IF (KOUNT2 .EQ. 1) GO TO 157 DO 156 I1=1,NEC IF (ABS((ORTEM(I1,KOUNT2)-ORTEM(I1,KOUNT2-1))/ : ORTEM(I1,KOUNT2)) .GT. 1.0D-3) GO TO 157 156 CONTINUE GO TO 154 157 DO 125 I1=1,NEC IF ((ABS(ORTEM(I1,KOUNT2)) .GT. ACCY2) : .OR. (KOUNT2 .LE. NPTS)) GO TO 770 125 CONTINUE GO TO 155 154 DO 126 LB=1,NEC DO 153 IM1=1,NPTS ECV3(IM1)=ECVTEM(LB,IM1+KOUNT2-NPTS) ORTH3(IM1)=ORTEM(LB,IM1+KOUNT2-NPTS) 153 CONTINUE CALL SORT(NPTS,ORTH3,ECV3) CALL EXTRAP(ORTH3,ECV3,NPTS,ZERO,ECV(LB)) 126 CONTINUE ISTP3=ISTP3+1 GO TO 770 120 CONTINUE * 155 WRITE (IPD,*) ' NUMBER OF ITERATIONS IN OUTER LOOP=',KOUNT-1 CALL NORMC (J,XX) CALL MAXMIN (J,1,PEAK) CALL ORTHOG(0) 810 FORMAT (1X,11(1P,1E11.3)) IF (NUMWR .NE. 0) THEN WRITE (IPD,314) NH(J),NUMWR WRITE (NUMWR) NH(J),NP(J),NAK(J),NCONT,Z,H,RNT,E(J), : PZY,QZY,MHALF,HCONT WRITE (NUMWR) (PFC(L,J),L=1,NCONT), : (QFC(L,J),L=1,NCONT) END IF RETURN END IF IF (XCA .LT. EPS) GO TO 9 IF (E(J)) 9,9,10 * * COMPUTE ESTIMATE OF EIGENVALUE, IF NO OTHER ESTIMATE IS * AVAILABLE, OR IF EXCHANGE TERMS ARE ZERO. * 9 WA = EIGEN (J) GO TO 11 * * USE EIGENVALUE FROM PREVIOUS ITERATION, IF ANY. * 10 WA = E(J) 11 IF (WA) 12,12,13 * * IF ESTIMATE IS NOT POSITIVE, USE E = 0.1. * 12 WA = TENTH * * USE ESTIMATE OF PZ FROM PREVIOUS ITERATION. * 13 AG = PZ(J) * * IF IN FINAL STAGE OF RECOVERY FROM PREVIOUS FAILURE, * EXTRAPOLATE ESTIMATES OF E AND PZ FROM PREVIOUS * TWO ITERATIONS. * IF ((IPW-1)*(IPY-1)) 15,15,14 14 AG = ASMT(IPW)+ASMT(IPW)-ASMT(IPW-1) WA = ESMT(IPW)+ESMT(IPW)-ESMT(IPW-1) * * SOLVE DIRAC RADIAL EQUATIONS. * 15 IF (ITC(4) .EQ. 1) CALL HPOT (J,XCA) CALL SOLV (J,VN,XCA,XCP,XCQ,AG,WA,MF) * * FAILURE 9 INDICATES TABLE TOO SHORT; * EXTEND BY TWO STEPS. * IF (MF-9) 18,16,18 16 N = N+2 IF (N .GT. NMAX) GO TO 48 * * SET ADDITIONAL TABLE ENTRIES TO ZERO. * DO 17 I = 1,NW PF(N,I) = ZERO QF(N,I) = ZERO PF(N-1,I) = ZERO QF(N-1,I) = ZERO 17 CONTINUE * * .... AND TRY AGAIN. * GO TO 8 * * TEST FOR OTHER TYPES OF FAILURE. * 18 IF (MF) 39,19,39 * * FIND LARGEST DIFFERENCE BETWEEN CALCULATED AND ESTIMATED * FUNCTIONS. * 19 WC = ZERO DO 21 I = 1,N WA = P(I)-PF(I,J) WB = ABS (WA)+ABS (Q(I)-QF(I,J)) IF (WB-WC) 21,21,20 20 WC = WB WD = WA IMX = I 21 CONTINUE * * AVOID ADJUSTMENT OF MIXING COEFFICIENT XAM UNLESS PREVIOUS * ITERATION INVOLVED THE SAME FUNCTION. * IF ((JL .NE. J) .OR. (IPY .EQ. 1) .OR. (IPWL .NE. IPW)) GO TO 30 * * INCREASE XAM IF NEW CORRECTION HAS OPPOSITE SIGN TO * PREVIOUS ONE, AND MORE THAN HALF THE MAGNITUDE. * IF (WC-ABS (ER(J))*HALF) 30,30,23 23 IF (WD*ER(J)) 24,27,27 24 IF (XAM(J)-PEIGHT) 26,26,25 * * DO NOT INCREASE XAM BEYOND 0.9 . * 25 XAM(J) = PNINE GO TO 30 26 XAM(J) = XAM(J)+TENTH GO TO 30 * * DECREASE XAM BY 0.1 IF NEW CORRECTION HAS SAME SIGN AS * PREVIOUS ONE, AND MORE THAN HALF THE MAGNITUDE. * 27 IF (XAM(J)-TENTH) 28,28,29 * * DO NOT DECREASE XAM BELOW ZERO. * 28 XAM(J) = ZERO GO TO 30 29 XAM(J) = XAM(J)-TENTH * * PRINT ONE LINE RECORD OF ITERATION. * 30 CONTINUE ER(J) = SIGN(WC,WD) WB = RNT*EPH**(IMX-1) WRITE (IPD,300) J,NP(J),NH(J),E(J),Z,XAM(J),XCA,WB,ER(J) * * TABULATE NEW ESTIMATE AS COMBINATION OF OLD AND COMPUTED * FUNCTIONS. * WB = XAM(J) WC = ONE-WB DO 32 I = 1,N PF(I,J) = WC*P(I)+WB*PF(I,J) QF(I,J) = WC*Q(I)+WB*QF(I,J) IF (ABS (PF(I,J))+ABS (QF(I,J))-EPS) 31,32,32 31 PF(I,J) = ZERO QF(I,J) = ZERO 32 CONTINUE * * RE-NORMALISE THE RESULT. * WA = ONE/SQRT (RINT (J,J,0)) DO 33 I = 1,N PF(I,J) = PF(I,J)*WA QF(I,J) = QF(I,J)*WA 33 CONTINUE * IPWL = IPW * * IF THIS COMPLETES STAGE 1 OF THE RECOVERY PROCESS, * BEGIN STAGE 2, COUNTING UPWARDS FROM ZERO IN IPW. * IF (IPY-1) 45,35,36 35 IPY = 2 ZMX = Z IPW = 0 36 IF (IPW-10) 37,45,45 * * STORE VALUES OF E AND PZ FOR SUBSEQUENT EXTRAPOLATION. * 37 ASMT(IPW+1) = PZ(J) ESMT(IPW+1) = E(J) ITR = ITR-1 * * REPEAT THIS STEP IF CHANGE WAS SIGNIFICANT. * IF (ABS (ER(J))-0.01) 38,38,8 * * ... OTHERWISE REDUCE Z AND REPEAT. * 38 IPW = IPW+1 XCA = XCA+XINC Z = ZST+(ZMX-ZST)*FLOAT (10-IPW) ZFAC = Z/ZST GO TO 8 * * IN CASE OF FAILURE, REPEAT COMPUTATION WITH EXTRA PRINTING * IF OPTION 13 IS SET. * 39 IF (ITC(13)) 41,41,40 40 ITCA = ITC(19) ITCB = ITC(18) ITCC = ITC(17) ITCD = ITC(15) ITC(19) = 1 ITC(18) = 1 ITC(17) = 1 ITC(15) = 1 CALL YPOT (J,VN,ZFAC) CALL XPOT (J,XCA,XCP,XCQ) IF (ITC(4) .EQ. 1) CALL HPOT (J,XCA) CALL SOLV (J,VN,XCA,XCP,XCQ,AG,WA,MF) ITC(19) = ITCA ITC(18) = ITCB ITC(17) = ITCC ITC(15) = ITCD 41 IF (IPY-1) 42,43,44 * * BEGIN RECOVERY BY SETTING EXCHANGE TERMS TO ZERO. * 42 IPY = 1 XCA = ZERO JL = J GO TO 8 * * IF FAILURE WITH IPY = 1 INCREASE Z UNTIL CONVERGENCE. * 43 Z = Z+TENTH ZFAC = Z/ZST IF (Z .LE. ZST+0.1*ZST) GO TO 8 MF = 5 * * IF FAILURE WITH IPY = 2 ABANDON CALCULATION. * 44 WRITE (IPD,305) NP(J),NH(J) IF (MF .EQ. -2) WRITE (IPD,306) IF (MF .EQ. -1) WRITE (IPD,307) IF (MF .EQ. 1) WRITE (IPD,308) IF (MF .EQ. 2) WRITE (IPD,309) IF (MF .EQ. 3) WRITE (IPD,310) IF (MF .EQ. 4) WRITE (IPD,311) IF (MF .EQ. 5) WRITE (IPD,312) CALL TIMER (0) STOP 4 * * ON FINAL SUCCESS, RESET Z TO TRUE VALUE AND TRY ANOTHER * ITERATION. * 45 Z = ZST GO TO 3 * 47 MX = 1 IF (NUMORB .NE. 0) THEN WRITE (IPD,*) ' TOO MANY ITERATIONS IN SR SCF' STOP END IF WRITE (IPD,303) IF (ITC(12) .EQ. 0) RETURN CALL TIMER (0) STOP 2 48 WRITE (IPD,304) CALL TIMER (0) STOP 3 * 300 FORMAT (1X,I2,1X,I2,A2,1P,E18.10,1X,0PF7.2,F5.1,1P,E10.1, : 1P,2E12.4) 301 FORMAT (/1X,10(1H+),5H SCF ,10(1H+)// : 13X,10HEIGENVALUE,8X,1HZ,4X,3HMIX,5X,3HXCA, : 7X,4HRMAX,8X,4HDMAX/) 303 FORMAT (/1X,24HITERATION LIMIT EXCEEDED) 304 FORMAT (/1X,24HINSUFFICIENT GRID POINTS) 305 FORMAT (/1X,27HFAILURE IN SCF PROCESS FOR ,I2,A2) 306 FORMAT (/1X,21HWRONG NUMBER OF ZEROS) 307 FORMAT (/1X,29HNUMBER OF ITERATIONS EXCEEDED) 308 FORMAT (/1X,29HSOLUTION CANNOT BE NORMALISED) 309 FORMAT (/1X,34HNORMALISED SOLUTION HAS WRONG SIGN) 310 FORMAT (/1X,22HNO POSITIVE EIGENVALUE) 311 FORMAT (/1X,29HMATCHING POINT IS TOO FAR OUT) 312 FORMAT (/1X,35HZ INCREASED BEYOND REASONABLE LIMIT) 313 FORMAT (/1X,' USING ESTIMATES FOR SLOPE AND ENERGY FROM ORBITAL', : 2X,1I2,1A2) 314 FORMAT (/1X,' WRITING CONTINUUM ORBITAL',1X,1A2,1X, : 'TO TAPE',1X,1I2,/) 315 FORMAT (//37H INCLUDE LAGRANGE MULTIPLIERS BETWEEN/) 316 FORMAT (16X,I2,A2,2X,I2,A2) 317 FORMAT (//44H LAGRANGE MULTIPLIERS ARE NOT TO BE INCLUDED) 318 FORMAT (//38H LAGRANGE MULTIPLIERS ARE NOT REQUIRED) 319 FORMAT (/1X,23HDIMENSION ERROR FOR NEC) END ************************************************************************ * * SUBROUTINE MAXMIN (J,IOUT,PEAK) * * * THIS SUBROUTINE FINDS THE MAXIMA AND MINIMA OF ARRAY PFC. * * * * INPUTS: * * J : SERIAL NUMBER OF ORBITAL * * IOUT : PRINTING OPTION * * * * OUTPUTS: * * PEAK : MAXIMUM AMPLITUDE OF PFC(I,J) * * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /INFORM/IRD,IPD COMMON/WAVEC/PFC(2000,25),QFC(2000,25) : /GRIDC/RGRIDC(2000),ECONT,HCONT,MHALF,NCONT,NUMWR * IF (IOUT .EQ. 1) WRITE (IPD,300) IBIT = 1 PP2MAX = 0.0 DO 72 M = 1,NCONT-1 IF ((ABS(PFC(M,J)) .GT. PP2MAX) .OR. (ABS(PFC(M+1,J)) .GT. : PP2MAX)) THEN IF (IBIT .EQ. 1) PP2MAX = ABS(PFC(M,J)) IF (PFC(M+1,J)*PFC(M,J) .LT. 0.0) IBIT = 1 ELSE IF ((PFC(M-1,J) .GT. 0.0) .AND. (IBIT .EQ. 1)) THEN IF (IOUT .EQ. 1)WRITE (IPD,301) PFC(M-1,J),RGRIDC(M-1) IBIT = 0 PEAK = PFC(M-1,J) PP2MAX = 0.0 ELSEIF ((PFC(M-1,J) .LT. 0.0) .AND. (IBIT .EQ. 1)) : THEN IF (IOUT .EQ. 1)WRITE (IPD,302) PFC(M-1,J),RGRIDC(M-1) PEAK = ABS(PFC(M-1,J)) IBIT = 0 PP2MAX = 0.0 END IF END IF 72 CONTINUE RETURN 300 FORMAT (/,' MAXIMA',9X,'MINIMA',8X,'R') 301 FORMAT (1X,1P,1E12.6,16X,1E12.6) 302 FORMAT (13X,1P,2E14.6) END ************************************************************************ * * SUBROUTINE XPOTC (J,XCF,ALX,ALY) * * * THIS SUBROUTINE TABULATES THE EXCHANGE TERMS XTPC(R) AND XTQC(R) * * THE CASE OF THE CONTINUUM WAVEFUNCTIONS WITH A LINEAR GRID. * * THESE MODIFICATIONS MADE BY W. F. PERGER, MICHIGAN TECHNOLOGICAL * * UNIVERSITY, AUGUST, 1988. * * * * INPUTS: * * J : SERIAL NUMBER OF ORBITAL * * XCF : SCALE FACTOR FOR MULTIPYING EXCHANGE TERMS * * * * OUTPUTS: * * XTPC,XTQC : COMMON ARRAYS CONTAINING EXCHANGE TERMS * * ALX, ALY : LEADING COEFFICIENTS IN THE EXPANSIONS OF THE * * EXCHANGE TERMS IN POWERS OF R * * * * SUBROUTINES CALLED : XA, XB, QUAD, YZKC. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * COMMON/CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /DEF1/Z,RNT,H,EPH,ACCY,C,N,NWA : /ORB4/NW,NCF,NH(25),NP(25),NAK(25),IQ(25,30) : /ORB5/NKL(25),NKJ(25) : /WAVE/PF(350,25),QF(350,25) : /POTE/YP(350),YQ(350),XTP(350),XTQ(350) : /TATB/TA(360),TB(360) : /HMAT/EAV,EMT(30,30),CCF(30),UCF(25),SA(30) : /FIXD/JFIX(25) : /EXCO/PZ(25),QZ(25) : /MCPA/XSLDR(500),ISLDR(500),NMCP : /OFFD/ECV(50),IECC(50),NEC : /MCPB/NNLDR(30,30),NSLDF(30,30) : /MCPC/CXCF(100),JXCF(100),NSCF,NSCFY : /SEMI/CHK(870),CCR(870) : /DEF7/NCMIN,ICCMIN(30) COMMON/WAVEC/PFC(2000,25),QFC(2000,25) : /POTEC/YPC(2000),YQC(2000),XTPC(2000),XTQC(2000) : /GRIDC/RGRIDC(2000),ECONT,HCONT,MHALF,NCONT,NUMWR : /TATBC/TAC(2000),TBC(2000) : /GRID/RGRID(350) : /KNTR/ITC(20) * DIMENSION CC(2000,3) * EPS = 1.0E-10 * * CLEAR FOR ACCUMULATION OF SUMS. * ALX = ZERO ALY = ZERO DO 1 I = 1,NCONT XTPC(I) = ZERO XTQC(I) = ZERO 1 CONTINUE IF (XCF .LT. EPS) RETURN * * SUM OVER ALL ORBITALS, AVOIDING SERIAL NUMBER J. * DO 10 L = 1,NW IF ((JFIX(J) .EQ. 1) .AND. (JFIX(L) .EQ. 0)) GO TO 10 IF (L .EQ. J) GOTO 10 * * DETERMINE BOUNDS FOR SUMMATION OVER K. * KM = (NKJ(J)+NKJ(L))/2 KN = ABS (NKJ(J)-NKJ(L))/2 IF (NAK(L)*NAK(J) .LT. 0) KN = KN+1 DO 5 K = KN,KM,2 WA = XA (J,L,K) WA = WA*XCF/C IF (ABS (WA) .LT. EPS) GOTO 5 * * TABULATE BASIC INTEGRAND. * DO 2 I = 1,NCONT TAC(I) = (PFC(I,L)*PFC(I,J)+QFC(I,L)*QFC(I,J))*WA/ : RGRIDC(I)**(K+1) 2 CONTINUE * CALL SPLINE (RGRIDC,TAC,NCONT,CC) * CALL SPLINT (RGRIDC,TAC,1,NCONT,NCONT,CC,RESULT,H,RNT) * ALX = ALX+RESULT*QZ(L) * ALY = ALY-RESULT*PZ(L) * * TABULATE FUNCTION YK(R) . * CALL YZKC (L,J,K,TBC,WA) * * ADD CONTRIBUTIONS FROM EXCHANGE TERMS. * DO 4 I = 1,NCONT XTPC(I) = XTPC(I)+TBC(I)*QFC(I,L) XTQC(I) = XTQC(I)-TBC(I)*PFC(I,L) 4 CONTINUE 5 CONTINUE * * ADD CONTRIBUTIONS FROM OFF-DIAGONAL PARAMETERS. * IF (NEC .EQ. 0) GO TO 10 KA = NWA*J+L IF (J-L) 6,10,7 6 KA = NWA*L+J 7 DO 9 I = 1,NEC IF (IECC(I) .NE. KA) GOTO 9 WA = ECV(I)*XCF/(UCF(J)*C) IF (ABS (WA) .LT. EPS) GO TO 9 * * ADD CONTRIBUTION TO LEADING COEFFICIENTS. * ALX = ALX+WA*QZ(L) ALY = ALY-WA*PZ(L) * * ADD CONTRIBUTIONS TO EXCHANGE TERMS. * DO 8 K = 1,NCONT XTPC(K) = XTPC(K)+RGRIDC(K)*WA*QFC(K,L) XTQC(K) = XTQC(K)-RGRIDC(K)*WA*PFC(K,L) 8 CONTINUE 9 CONTINUE 10 CONTINUE * * ADD CONTRIBUTIONS FROM OPEN SHELL ORBITALS. * IF (NCMIN .EQ. 0) RETURN IF (NMCP .EQ. 0) RETURN NSCF = 0 NCFM = NCF-1 DO 15 ITR = 1,NCFM ICFP = ITR+1 DO 15 JTR = ICFP,NCF ILDN = NNLDR(ITR,JTR) IF (ILDN .LE. 0) GOTO 15 JAA = NSLDF(ITR,JTR) JBB = JAA+ILDN-1 DO 14 JJ = JAA,JBB * * EXTRACT QUANTUM NUMBER FROM CODED FORMS. * KX = ISLDR(JJ) JLA = MOD (KX,NWA) KX = KX/NWA JLB = MOD (KX,NWA) KX = KX/NWA JLC = MOD (KX,NWA) IF (JLC .EQ. 0) GO TO 14 KX = KX/NWA JLD = MOD (KX,NWA) K = KX/NWA IF (NCMIN .GT. 1) GO TO 11 WA = XSLDR(JJ)*CCF(ITR)*CCF(JTR) GO TO 13 11 TOT = ZERO IA = ITR IB = JTR DO 12 I = 1,NCMIN TOT = TOT+CCR(IA)*CCR(IB) IA = IA+NCF IB = IB+NCF 12 CONTINUE WA = XSLDR(JJ)*TOT/FLOAT (NCMIN) * * ACCUMULATE COEFFICIENTS FROM EACH OF THE FOUR PARAMETERS * IN TURN. * 13 IF (JLA .EQ. J) CALL XB (JLC,JLB,JLD,K,WA) IF (JLB .EQ. J) CALL XB (JLD,JLA,JLC,K,WA) IF (JLC .EQ. J) CALL XB (JLA,JLB,JLD,K,WA) IF (JLD .EQ. J) CALL XB (JLB,JLA,JLC,K,WA) 14 CONTINUE 15 CONTINUE IF (NSCF .LE. 0) RETURN * * RETRIEVE SUCCESSIVE TERMS. * DO 18 ICF = 1,NSCF KX = JXCF(ICF) JLC = MOD (KX,NWA) KX = KX/NWA JLB = MOD (KX,NWA) KX = KX/NWA JLD = MOD (KX,NWA) K = KX/NWA IF ((JFIX(J) .EQ. 1) .AND. (JFIX(JLB) .EQ. 0)) GO TO 18 IF ((JFIX(J) .EQ. 1) .AND. (JFIX(JLC) .EQ. 0)) GO TO 18 IF ((JFIX(J) .EQ. 1) .AND. (JFIX(JLD) .EQ. 0)) GO TO 18 * * DETERMINE COEFFICIENT. * WA = CXCF(ICF)*XCF/(C*UCF(J)) IF (ABS (WA) .LT. EPS) GOTO 18 * * TABULATE YK(R). * CALL YZKC (JLB,JLD,K,TBC,WA) * * ADD CONTRIBUTIONS TO EXCHANGE TERMS. * DO 17 I = 1,NCONT XTPC(I) = XTPC(I)+TBC(I)*QFC(I,JLC) XTQC(I) = XTQC(I)-TBC(I)*PFC(I,JLC) 17 CONTINUE 18 CONTINUE RETURN END ************************************************************************ * * SUBROUTINE YZKC (I1,I2,K,W3,WA) * * * * * THIS SUBROUTINE RETURNS THE HARTREE Y FUNCTION IN ARRAY W3, FOR * * ORIBITALS I1 AND I2 AND INDEX K. THIS IS AN ADAPTATION OF A SUB- * * ROUTINE ORIGINALLY WRITTEN BY HANS ANDRIESSEN, DELFT, NETHERLANDS. * * * * SUBROUTINES CALLED: NONE. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * COMMON/WAVEC/PFC(2000,25),QFC(2000,25) : /DEF1/Z,RNT,H,EPH,ACCY,C,N,NWA : /GRIDC/RGRIDC(2000),ECONT,HCONT,MHALF,NCONT,NUMWR : /INFORM/IRD,IPD : /KNTR/ITC(20) : /ORB4/NW,NCF,NH(25),NP(25),NAK(25),IQ(25,30) * DIMENSION W1(2000),W3(2000),W4(2000),GRAL(2000), : ORD(2000),PI(2000),QI(2000) * XKA = DBLE (K+1) C1 = SQRT (NAK(I1)*NAK(I1) - (Z/C)*(Z/C)) C2 = SQRT (NAK(I2)*NAK(I2) - (Z/C)*(Z/C)) MM = MHALF DO 622 M = 1,NCONT PI(M) = PFC(M,I1) QI(M) = QFC(M,I1) GRAL(M) = PFC(M,I2) ORD(M) = QFC(M,I2) 622 CONTINUE DO 40 M=1,NCONT W1(M)=(GRAL(M)*PI(M)+ORD(M)*QI(M))*WA 40 CONTINUE DO 50 M=1,MM ORD(M)=(RGRIDC(M)**XKA)*W1(M) 50 CONTINUE GRAL(1)=ORD(1)/(C1+C2+XKA) GRAL(2)=ORD(2)/(C1+C2+XKA) DO 60 M=3,MM GRAL(M)=GRAL(M-2)+H*(ORD(M-2)+4.0D0*ORD(M-1)+ORD(M))/3.0D0 60 CONTINUE DO 70 M=1,MM W1(M)=GRAL(M)/RGRIDC(M)**(XKA-1.0) 70 CONTINUE A1=0.D0 A2=0.D0 GRAL(NCONT)=0.D0 M1=MM+1 DO 80 M=M1,NCONT ORD(M)=RGRIDC(M)**(XKA-1.D0)*W1(M) 80 CONTINUE HW=GRAL(MM) OH=ORD(MM)/RGRIDC(MM) GRAL(MM+1)=GRAL(MM)+HCONT*(5.D0*OH+8.D0*ORD(MM+1)-ORD(MM+2))/ : 12.00D0 M1=MM+2 HW=ORD(MM) ORD(MM)=OH DO 90 M=M1,NCONT GRAL(M)=GRAL(M-2)+HCONT*(ORD(M-2)+4.D0*ORD(M-1)+ORD(M))/3.D0 90 CONTINUE M1=MM+1 DO 100 M=M1,NCONT W1(M)=GRAL(M)/RGRIDC(M)**XKA 100 CONTINUE DO 110 M=MM,NCONT ORD(M)=ORD(M)/RGRIDC(M)**(2.D0*XKA-1.0D0) 110 CONTINUE GRAL(NCONT)=0.D0 GRAL(NCONT-1)=HCONT*(5.D0*ORD(NCONT)+8.D0*ORD(NCONT-1) : -ORD(NCONT-2))/12.0D0 M=NCONT-1 120 M=M-1 GRAL(M)=GRAL(M+2)+HCONT*(ORD(M+2)+4.D0*ORD(M+1)+ORD(M))/3.D0 IF (M .GT. MM) GO TO 120 ORD(MM)=HW DO 140 M=1,MM ORD(M)=ORD(M)/RGRIDC(M)**(2.D0*XKA-1.0) 140 CONTINUE OH=ORD(MM) GRAL(MM-1)=GRAL(MM)+H*(5.D0*OH+8.D0*ORD(MM-1)-ORD(MM-2))/12.0D0 M=MM-1 150 M=M-1 GRAL(M)=GRAL(M+2)+H*(ORD(M+2)+4.D0*ORD(M+1)+ORD(M))/3.D0 IF (M .GT. 1) GO TO 150 DO 160 M=1,MM W4(M)=GRAL(M)*RGRIDC(M)**XKA 160 CONTINUE M1=MM+1 DO 170 M=M1,NCONT W4(M)= RGRIDC(M)**(XKA-1.D0)*GRAL(M) 170 CONTINUE DO 175 M=1,NCONT IF (M .LE. MHALF) THEN W3(M) = W1(M)+W4(M) ELSE W3(M) = (W1(M)+W4(M))*RGRIDC(M) END IF 175 CONTINUE RETURN END ************************************************************************ * * SUBROUTINE SPLINE (X,Y,N,C) * * * THIS ROUTINE IS AN ADAPTATION OF THE CUBIC SPLINES INTERPOLATION * * ROUTINE BY FORSYTHE ET AL., "COMPUTER METHODS FOR MATHEMATICAL * * COMPUTATIONS," (1977) BY W.F. PERGER, MICHIGAN TECH UNIV, DEC. 1988* * * * NO SUBROUTINES CALLED. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N DOUBLE PRECISION X(N),Y(N),C(N,3) * * THE COEFFICIENTS C(I,1), C(I,2), AND C(I,3), I=1,2,...N ARE COMPUTED * FOR A CUBIC INTERPOLATING SPLINE * * S(X)= Y(I) + C(I,1)*(X-X(I)) + C(I,2)*(X-X(I))**2 + C(I,3)*(X-X(I))**3 * * FOR X(I) .LE. X .LE. X(I+1) * * INPUT: * * N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2) * X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER * Y = THE ORDINATES OF THE KNOTS * * OUTPUT: * * C(I,J) = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE. * * USING P TO DENOTE DIFFERENTIATION * * Y(I) = S(X(I)) * C(I,1) = SP(X(I)) * C(I,2) = SPP(X(I))/2 * C(I,3) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT) * * THE ACCOMPANYING FUNCTION SUBPROGRAM SPLIN IS USED TO EVALUATE * THE SPLINE * INTEGER NM1,IB,I * NM1 = N-1 IF (N .LT. 2) RETURN IF (N .LT. 3) GO TO 50 * * SET UP TRIDIAGONAL SYSTEM * * C(I,1) = DIAGONAL, C(I,3) = OFFDIAGONAL, C(I,2) = RIGHT HAND SIDE * C(1,3) = X(2) - X(1) C(2,2) = (Y(2) - Y(1))/C(1,3) DO 10 I = 2,NM1 C(I,3) = X(I+1) - X(I) C(I,1) = 2.0D0*(C(I-1,3) + C(I,3)) C(I+1,2) = (Y(I+1) - Y(I))/C(I,3) C(I,2) = C(I+1,2) - C(I,2) 10 CONTINUE * * END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N) * OBTAINED FROM DIVIDED DIFFERENCES. * C(1,1) = -C(1,3) C(N,1) = -C(N-1,3) C(1,2) = 0.0D0 C(N,2) = 0.0D0 IF (N .EQ. 3) GO TO 15 C(1,2) = C(3,2)/(X(4)-X(2)) - C(2,2)/(X(3)-X(1)) C(N,2) = C(N-1,2)/(X(N)-X(N-2)) - C(N-2,2)/(X(N-1)-X(N-3)) C(1,2) = C(1,2)*C(1,3)**2.0D0/(X(4)-X(1)) C(N,2) = -C(N,2)*C(N-1,3)**2.0D0/(X(N)-X(N-3)) * * FORWARD ELIMINATION * 15 DO 20 I=2,N T = C(I-1,3)/C(I-1,1) C(I,1) = C(I,1) - T*C(I-1,3) C(I,2) = C(I,2) -T*C(I-1,2) 20 CONTINUE * * BACK SUBSTITUTION * C(N,2) = C(N,2)/C(N,1) DO 30 IB =1,NM1 I = N-IB C(I,2) = (C(I,2)-C(I,3)*C(I+1,2))/C(I,1) 30 CONTINUE * * C(I,2) IS NOW THE SIGMA(I) OF THE TEXT * * COMPUTE POLYNOMIAL COEFFICIENTS * C(N,1) = (Y(N) - Y(NM1))/C(NM1,3)+C(NM1,3)*(C(NM1,2)+2.0D0*C(N,2)) DO 40 I = 1,NM1 C(I,1) = (Y(I+1) - Y(I))/C(I,3)-C(I,3)*(C(I+1,2)+2.0D0*C(I,2)) C(I,3) = (C(I+1,2)-C(I,2))/C(I,3) C(I,2) = 3.0D0*C(I,2) 40 CONTINUE C(N,2) = 3.0D0*C(N,2) C(N,3) = C(N-1,3) RETURN * 50 C(1,1) = (Y(2)-Y(1))/(X(2)-X(1)) C(1,2) = 0.0D0 C(1,3) = 0.0D0 C(2,1) = C(1,1) C(2,2) = 0.0D0 C(2,3) = 0.0D0 RETURN END ************************************************************************ * * SUBROUTINE SPLINT (X,Y,JMIN,JMAX,N,YY,RES,HX,RNT) * * * WRITTEN BY WARREN F. PERGER, MICHIGAN TECHNOLOGICAL UNIVERSITY, * * JULY, 1988 * * * * THIS SUBROUTINE IS AN ADAPTATION OF EQUATION 2.3.2.2, "METHODS OF * * NUMERICAL INTEGRATION," BY PHILIP J. DAVIS AND PHILIP RABINOWITZ, * * (1975). IN ADDITION, IT ALSO ESTIMATES THE INTEGRAL BETWEEN R=0 * * AND R=RNT (THE FIRST GRID POINT) BY ASSUMING THAT THE INTEGRAND HAS * * THE FORM R*F(R)=W*(R**PX), WHERE PX IS DETERMINED BY USING THE FIRST * * TWO POINTS IN THE INTEGRAND. THE PROCEDURE IS STARTS BY WRITING THE * * INTEGRAL FROM 0 TO RNT OF W*R**PX = (W*RNT**(PX+1))/(PX+1). THEN, * * F(I=1) = W*RNT**PX, F(I=2) = W*((RNT**PX)*(EXP(H*PX)). TAKING * * F(I=2)/F(I=1) YIELDS PX = (1/H)*LN(F(I=2)/F(I=1)). THE INTEGRAL * * OF F(R)DR FROM 0 TO RNT THEN BECOMES : * * F(I=1)*RNT/(1+(1/H)*LN(F(I=2)/F(I=1))) * * * * INPUTS: X THE ARRAY OF ABSCISSAS. * * Y THE ARRAY OF THE FUNCTION TO BE INTEGRATED AT THE * * LOCATIONS SPECIFIED BY X. * * YY THE ARRAY OF CUBIC SPLINES, WITH YY(I,3) BEING THE * * SECOND DERIVATIVE OF Y(X). * * JMIN THE LOCATION IN X OF THE BEGINNING OF THE INTERVAL. * * JMAX THE LOCATION IN X OF THE END OF THE INTERVAL OF INTEG. * * N THE NUMBER OF POINTS IN X AND Y. * * HX FROM THE EXPONENTIAL GRID, R(I) = RNT*EXP((I-1)*HX) * * RNT (SEE HX) * * * * OUTPUTS: RES THE RESULT OF THE INTEGRATION FROM JMIN TO JMAX * * * * NO SUBROUTINES CALLED. * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * DIMENSION X(N),Y(N),YY(N,3) RES = 0.0D0 DO 10 J = 1,N IF (ABS(Y(J)) .GT. 1.0D-30) GO TO 4 10 CONTINUE RETURN 4 JJ = J+1 DO 1 I = JMIN,JMAX-1 H = X(I+1)-X(I) RES = RES + (Y(I)+Y(I+1))*H/2.0D0 : - (YY(I,3)+YY(I+1,3))*H*H*H/24.D0 1 CONTINUE COP = Y(JJ)/Y(J) IF ((COP .GT. 1.0D0) .AND. (HX .GT. 1.0E-10)) GO TO 6 5 RPX = 0.0D0 * WRITE (6,300) GO TO 7 6 CONTINUE RPX = LOG(COP)/HX 7 RES = RES+Y(J)*RNT/(RPX+1.0D0) * * 300 FORMAT (1X,44HSPLINT : CORRECTION NEAR NUCLEUS SET TO ZERO) * RETURN END ************************************************************************ * * SUBROUTINE SPLIN (XA,YA,N,CC,XB,YB,NB) * * * WRITTEN BY WARREN F. PERGER, MICHIGAN TECHNOLOGICAL UNIVERSITY, * * DECEMBER, 1988 * * * * THIS SUBROUTINE IS AN ADAPTATION OF A ROUTINE FROM THE BOOK BY * * FORSYTHE, ET AL., "COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS." * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * COMMON/INFORM/IRD,IPD DIMENSION XA(N),YA(N),CC(N,3),XB(NB),YB(NB) DATA I/1/ IF (XA(N) .LT. XB(NB)) WRITE (IPD,300) XA(N),XB(NB) DO 40 J=1,NB IF (I .GE. N) I = 1 X = XB(J) IF (X .LT. XA(I)) GO TO 10 IF (X .LE. XA(I+1)) GO TO 30 * * BINARY SEARCH * 10 I = 1 L = N+1 20 K = (I+L)/2 IF (X .LT. XA(K)) L = K IF (X .GE. XA(K)) I = K IF (L .GT. I+1) GO TO 20 * * EVALUATE SPLINE * 30 DX = X-XA(I) YB(J) = YA(I)+DX*(CC(I,1)+DX*(CC(I,2)+DX*CC(I,3))) 40 CONTINUE RETURN 300 FORMAT (/,1X,1P,'WARNING--ERROR IN SR SPLIN: XA(N)=',1X,1E11.3,/, : ' AND XB(NB)=',1X,1E11.3,//) END ************************************************************************ * * SUBROUTINE SOLVC (J,VN,XCF,ALX,ALY,AST,EST,MX) * * * THIS SUBROUTINE SOLVES A SINGLE PAIR OF DIRAC RADIAL EQUATIONS. * * THE DIRECT POTENTIALS ARE TABULATED IN YP AND YQ . * * THE EXCHANGE POTENTIALS ARE TABULATED IN XTP AND XTQ . * * THESE ARRAYS ARE OVERWRITTEN IN SOLV SO THEY MUST BE SET UP * * FOR THE NEXT CALL TO SOLV . * * * * INPUTS: * * * * J : SERIAL NUMBER OF WAVE FUNCTION * * VN : COEFFICIENT OF R IN SERIES EXPANSION OF DIRECT POTENTIAL * * XCF: SCALE FACTOR MULTIPLYING EXCHANGE TERM * * ALX: LEADING COEFFICIENT IN SERIES EXPANSION OF EXCHANGE * * POTENTIAL FOR P EQUATION * * ALY: LEADING COEFFICIENT IN SERIES EXPANSION OF EXCHANGE * * POTENTIAL FOR Q EQUATION * * AST: ESTIMATE OF A * * EST: ESTIMATE OF EIGENVALUE E * * * * OUTPUTS: * * * * MX : SHOULD BE ZERO. A NON-ZERO VALUE INDICATES * * A FAILURE. * * E(J) THE EIGENVALUE * * * * THE WAVE FUNCTIONS WILL BE TABULATED IN ARRAYS * * P AND Q, AND THE COEFFICIENTS IN THE SERIES * * EXPANSIONS WILL BE STORED IN PZ(J), QZ(J). * * * * SUBROUTINES CALLED: * * * * PRWF, OUT, IN, MATCH, QUAD. * * * ************************************************************************ * IMPLICIT DOUBLE PRECISION (A-H,O-Z) * COMMON/INFORM/IRD,IPD : /WAVE/PF(350,25),QF(350,25) : /CONS/ZERO,HALF,TENTH,ONE,TWO,THREE,TEN : /DEF1/Z,RNT,H,EPH,ACCY,C,N,NWA : /ORB4/NW,NCF,NH(25),NP(25),NAK(25),IQ(25,30) : /DEF2/NCFX,NMAX,NFMAX,NMCPX,NECX,NSCFX : /ORB1/E(25),GAMA(25),XAM(25) : /POTE/YP(350),YQ(350),XTP(350),XTQ(350) : /TATB/TA(360),TB(360) : /INT1/TC(350),TD(350),TE(350),TF(350), : TG(350),TH(350) : /INT2/P(350),Q(350),PC(350),QC(350) : /INT3/XU(350),XV(350),XR(350),XS(350),XP(350) : /EXCO/PZ(25),QZ(25) : /KNTR/ITC(20) COMMON/WAVEC/PFC(2000,25),QFC(2000,25) : /POTEC/YPC(2000),YQC(2000),XTPC(2000),XTQC(2000) : /GRIDC/RGRIDC(2000),ECONT,HCONT,MHALF,NCONT,NUMWR : /GRID/RGRID(350) : /TATBC/TAC(2000),TBC(2000) DIMENSION XPTEM(2000),XQTEM(2000),YPTEM(2000),YQTEM(2000), : YP1(350),YQ1(350),PTEM(350),QTEM(350) * IF (ITC(15) .EQ. 1) WRITE (IPD,303) NP(J),NH(J), : VN,XCF,ALX,ALY,AST,EST,Z,N CON = 0.3D0 MAU = 0 * * BEGIN FROM GIVEN ESTIMATES OF A AND E * ACL = AST ECL = E(J) * ELW = -ONE EHG = -ONE FKJ = FLOAT (NAK(J)) HH = -H*HALF HHCONT = -HCONT*HALF WA = -FKJ*HH WAC = ONE+WA WBC = ONE-WA XMOMEN = SQRT((-E(J)+C*C)*(-E(J)+C*C)-C*C*C*C)/C * * SAVE DIRECT & EXCH POTENTIALS FOR REPETITIVE CALLS TO OUT AND OUTC * DO 41 M = 1,N YP1(M) = YP(M) YQ1(M) = YQ(M) 41 CONTINUE DO 40 M = 1,NCONT XPTEM(M) = XTPC(M) XQTEM(M) = XTQC(M) 40 CONTINUE VALUE = ZERO KOUNT = 0 IX = 30 NLOOPS = 5 IF (XCF .LE. ZERO) GO TO 3 * * TABULATE MODIFIED EXCHANGE TERMS, INCLUDING DIFFERENCE * CORRECTIONS. * 42 KOUNT = KOUNT + 1 IF (KOUNT .GT. 50) THEN WRITE (IPD,*) ' NUMBER OF ITERATIONS EXCEEDED IN SOLVC' STOP END IF DO 45 M = 1,NCONT XTPC(M) = XPTEM(M) XTQC(M) = XQTEM(M) 45 CONTINUE DO 50 I = 1,MHALF+IX P(I) = PFC(I,J) Q(I) = QFC(I,J) PC(I) = P(I) QC(I) = Q(I) YP(I) = YP1(I) YQ(I) = YQ1(I) XU(I) = XTPC(I) XV(I) = XTQC(I) 50 CONTINUE DO 55 I = -IX/2,IX/2 PTEM(I+IX/2+1) = P(MHALF+I) QTEM(I+IX/2+1) = Q(MHALF+I) TD(I+IX/2+1) = XTPC(MHALF+I) TE(I+IX/2+1) = XTQC(MHALF+I) TH(I+IX/2+1) = RGRIDC(MHALF+I) 55 CONTINUE DO 60 I = -1,2 RINT = RGRID(MHALF+I) IF (I .LT. 0) RINT = RGRID(MHALF)-HCONT CALL EXTRAP (TH,PTEM,IX,RINT,PC(MHALF+I)) CALL EXTRAP (TH,QTEM,IX,RINT,QC(MHALF+I)) 60 CONTINUE DO 62 I = 0,1 RINT = RGRID(MHALF+I) CALL EXTRAP (TH,TD,IX,RINT,XU(MHALF+I)) CALL EXTRAP (TH,TE,IX,RINT,XV(MHALF+I)) 62 CONTINUE 1 DO 2 I = 4,MHALF+2 XU(I-2) = HH*(XU(I-1)+XU(I-2))+((PC(I-3)-PC(I)) : /12.0D0 -(PC(I-2)-PC(I-1))*0.25)*XCF XV(I-2) = HH*(XV(I-1)+XV(I-2))+((QC(I-3)-QC(I)) : /12.0D0 -(QC(I-2)-QC(I-1))*0.25)*XCF 2 CONTINUE DO 6 I = MHALF+2,NCONT H1 = RGRIDC(I-1) H2 = RGRIDC(I-2) XTPC(I-2) = HHCONT*(XTPC(I-1)/H1+XTPC(I-2)/H2)+((PFC(I-3,J) : -PFC(I,J))/12.0D0 -(PFC(I-2,J)-PFC(I-1,J))*0.25)*XCF XTQC(I-2) = HHCONT*(XTQC(I-1)/H1+XTQC(I-2)/H2)+((QFC(I-3,J) : -QFC(I,J))/12.0D0 -(QFC(I-2,J)-QFC(I-1,J))*0.25)*XCF 6 CONTINUE * * TABULATE MODIFIED DIRECT POTENTIALS. * 3 WD = HH/C DO 4 I = 1,MHALF