SUBROUTINE PREPOT 8/17R79 C C POTA.FOR HAS ANTI-MORSE BEND PLUS BOWMAN'S COLLINEAR SPLINE ROUTINE. C ANTI-MORSE BENDING POTENTIAL IS ADDED TO COLLINEAR POTENTIAL. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GENERAL INFORMATION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS SUBROUTINE CAN BE USED TO ADD A BENDING CORRECTION TO C ANY COLLINEAR POTENTIAL SUBROUTINE. C C PREPOT MUST BE CALLED BEFORE ANY CALLS TO POT C CAN BE MADE. THE POTENTIAL PARAMETERS FOR THE COL- C LINEAR POTENTIAL ARE READ IN BY A CALL TO PREPOT2 C (FROM OLD PREBEND). THEN THE BENDING PARAMETERS ARE INPUT. C COORDINATES, POTENTIAL ENERGY, AND DERIVATIVES ARE C PASSED THROUGH COMMON /PT31CM/ R(3), ENERGY, DEDR(3). C C NO CALLS TO PREPOT2 OR ORIGINAL POTENTIAL ARE NEEDED IN THE C MAIN PROGRAM. C CCCCCCCCCCCCCCCCCCCCCCCCCCC C DEFINE VARIABLES CCCCCCCCCCCCCCCCCCCCCCCCCCC C C Potential parameters' default settings C Variable Default value C NSURF 0 C NDER 1 C IDEBUG 0 C IPRT 6 C IRD1 4 C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 NAB,NAB1,NBC,NBC1,JUNK,HARTREE,BOHR CHARACTER*80 DATAF 07AUG83 DIMENSION DBCOOR(2),JUNK(4),REQ(3) COMMON /PT31CM/ R(3), ENERGY, DEDR(3) COMMON /PT35CM/ EASYAB, EASYBC, EASYAC COMMON /PT34CM/ IPRT, IRD1 C IRD1 = 4 IPRT = 6 C C READ IN DATA FILE NAME 07AUG83 C READ(9,900,ERR=1110) DATAF 07AUG83 C 900 FORMAT(1X,A80) 07AUG83 C WRITE(IPRT,6101) DATAF 07AUG83 C6101 FORMAT(' POTENTIAL DATA FROM FILE',T40,A80) 07AUG83 C IF(DATAF.EQ.' ') GO TO 1110 07AUG83 C OPEN (UNIT=8,ERR=1100,FILE=DATAF,FORM='FORMATTED',READONLY, 07AUG83 C 1 STATUS='OLD') 07AUG83 OPEN (UNIT=IRD1, FILE='potoh2dim.dat', FORM='FORMATTED', * STATUS='OLD', ERR=1100) GO TO 1110 07AUG83 1100 WRITE(IPRT,6100) 07AUG83 6100 FORMAT(' PROBLEM WITH POTENTIAL DATA FILE, CHECK FILE NAME') 07AUG83 CALL EXIT 07AUG83 1110 CONTINUE 07AUG83 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C READ COLLINEAR POTENTIAL PARAMETERS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C WRITE(IPRT,2) 2 FORMAT(//T30,30HCOLLINEAR POTENTIAL PARAMETERS///) C CALL PREPOT2 8/17R79 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C READ BENDING POTENTIAL PARAMETERS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C WRITE(IPRT,3) 3 FORMAT(//T30,28HBENDING POTENTIAL PARAMETERS///) C C DISSOCIATION ENERGY IN KCAL/MOLE C C READ(8,5) DE READ(IRD1,5) DE 5 FORMAT (3F20.10) WRITE(IPRT,10) DE 10 FORMAT(//1X,32H DISSOCIATION ENERGY (KCAL/MOLE),T50,3F20.10) C C C READ IN MORSE BETAS--RECIPROCAL ANGSTROMS C C READ(8,20) BETA READ(IRD1,20) BETA WRITE(IPRT,21) BETA 20 FORMAT(3F20.10) 21 FORMAT (/1X,33HMORSE BETA (RECIPROCAL ANGSTROMS),T50,3F20.10) C C READ IN THE EQUILIBRIUM DISTANCES IN ANGSTROMS C C READ(8,20) (REQ(IT), IT = 1,3) READ(IRD1,20) (REQ(IT), IT = 1,3) WRITE (IPRT,25) (REQ(IT), IT = 1,3) 25 FORMAT(/1X,33HEQUILIBRIUM DISTANCES (ANGSTROMS),T50,3F20.10) C C READ IN THE PAULING PARAMETER IN ANGSTROMS C C READ(8,5) A READ(IRD1,5) A WRITE (IPRT,30) A 30 FORMAT(/1X,29HPAULING PARAMETER (ANGSTROMS),T50,F20.10) C C READ IN GAMMA FOR THE BEND CORRECTION C GAMMA IS DIMENSIONLESS AND RANGES FROM .5 TO .65 C C READ(8,5) GAMMA READ(IRD1,5) GAMMA WRITE (IPRT,31) GAMMA 31 FORMAT (/1X,16HGAMMA (NO UNITS),T50,F20.10) C C READ IN RBB IN ANGSTROMS--THIS NUMBER IS USED C IN A SCALING CALCULATION FOR THE BENDING CORRECTION C C READ(8,5) RBB READ(IRD1,5) RBB WRITE(IPRT,35) RBB 35 FORMAT(/1X,15HRBB (ANGSTROMS),T50,F20.10) C CCCCCCCCCCCCCCCCCCCCCCCCCCC C USEFUL CONSTANTS CCCCCCCCCCCCCCCCCCCCCCCCCCC C C CONVERSTION FOR ANGSTROMS TO BOHR (ANGSTROM/BOHR) C BOHR = .52917706D0 C C CONVERSION FOR KCAL TO HARTREE (KCAL/HARTREE) C HARTREE = 627.5095D0 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CONVERT TO ATOMIC UNITS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 50 IT = 1,3 REQ(IT) = REQ(IT)/BOHR 50 CONTINUE BETA = BETA * BOHR DE = DE/HARTREE A = A/BOHR RHH = 1.4007977D0 RBB = RBB /BOHR EASYAB = 105.64D0 / HARTREE EASYBC = 109.489D0 / HARTREE EASYAC = 105.64D0 / HARTREE C CCCCCCCCCCCC RETURN CCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C MOVE ONTO THE MAIN PART OF THE ROUTINE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCC ENTRY POT 8/17R79 CCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FIRST CALCULATE SOME USEFUL NUMBERS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C NAB AND NBC ARE THE BOND ORDERS OF THE DIATOMICS C AB AND BC RESPECTIVELY. (EQUATION 2.2) C NAB = EXP((REQ(1) - R(1))/A) C IF (NAB.LT.1.D-15) NAB = 1.D-15 C NBC = EXP((REQ(2) - R(2))/A) C IF (NBC.LT.1.D-15) NBC = 1.D-15 C C CALCULATE BOND ORDERS ALONG THE REACTION COORDINATE C SEE EQUATION 2.7A, 2.7B C C = 1.D0 - NAB/NBC IF (ABS(C).GT.1.D-14) GOTO 90 NAB1=.5D0 NBC1=.5D0 GO TO 95 90 C = C/NAB NAB1 = (2.D0 + C - SQRT(4.D0+C*C))/(2.D0*C) NBC1 = 1.D0 - NAB1 C C ROB IS USED IN THE BENDING POTENTIAL CALCULATIONS C C FIRST TRAP OUT ANY ZERO ARGUEMENTS C 95 IF (NAB1*NBC1.GT.0.0D0) GO TO 103 ROB = 1.D0 GO TO 107 C 103 STUFF = A * DLOG(NAB1*NBC1) ROB = (RBB - STUFF)/(RHH - STUFF) C 107 CONTINUE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCULATE BENDING CORRECTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CALCULATE V(R(3)) C JUNK(3) = EXP(-BETA*(R(3)-REQ(3))/ROB) VAC = GAMMA *1.0D-14 * DE * JUNK(3)* (1.D0 + 0.5D0*JUNK(3)) C C CALCULATE V(R(1)+R(2)) C JUNK(4) = EXP(BETA*(REQ(3)-R(2)-R(1))/ROB) VABBC= GAMMA *1.0D-14 * DE * JUNK(4) * (1.D0 + 0.5D0*JUNK(4)) C C HERE IS THE BENDING CORRECTION C BCOOR = (VAC - VABBC)*1.0D14 C C WHILE IT'S CONVENIENT, CALCULATE THE DERIVATIVE C OF BCOOR WITH RESPECT TO R(1) R(2) AND R(3). C C CALCULATE DAC (THE DERIVATIVE WITH RESPECT TO R(3)) C DAC = -GAMMA*DE*BETA*JUNK(3)*(1.D0+JUNK(3))/ROB DBCOOR(1) = GAMMA * DE * BETA * JUNK(4) * (1.D0 + JUNK(4))/ROB DBCOOR(2)= DBCOOR(1) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NOW CALCULATE THE COLLINEAR ENERGY CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C STORE R(3) IN RACTEMP THEN SET R(3) TO THE COLLINEAR C GEOMETRY AND CALL POT C RACTEMP = R(3) R(3) = R(2) + R(1) C CALL POT2 8/17R79 C DEDR(1) = DEDR(1) + DEDR(3) DEDR(2) = DEDR(2) + DEDR(3) DEDR(3) = DAC R(3) = RACTEMP C ENERGY = ENERGY + BCOOR C DO 300 IT = 1,2 DEDR(IT) = DBCOOR(IT) + DEDR(IT) 300 CONTINUE RETURN END C SUBROUTINE PREPOT2 C C POTENTIAL FUNCTION FOR RATE1D SERIES C RMO-SPLINE TYPE FIT WITH BEBO TYPE CONTINUATION INTO THE ASYMBTOTIC REGION C POTENTIAL ROUTINE INTERACTS WITH REST OF PROGRAM WITH UNITS OF EV/BOHR C POTENTIAL ROUTINE ASSUMES SPLINE INFO IS IN UNITS OF KCAL/ANGSTROM C IMPLICIT REAL*8(A-H,O-Z) DIMENSION V(4),ARRAY(128) COMMON/NSTEP/STEP COMMON/PT31CM/R(3), ENERGY, DEDR(3) COMMON /PT34CM/ IPRT, IRD1 DIMENSION PHI(40),BSP(40,3),XBC(3),XAB(3),CBC(3),ABC(3), 1 CAB(3),AAB(3) DIMENSION CSP(48,3),SCRAP(48) DATA DETRAD/.0174532D0/ C C DIRECT INPUT INTO POTENTIAL ROUTINE C C READ (8,843)TITLE READ (IRD1,843)TITLE 843 FORMAT(5A4) C READ(8,2) STEP READ(IRD1,2) STEP 2 FORMAT(F12.10) WRITE(IPRT,845)TITLE 845 FORMAT(1X,' SPLINE COEF STORED IN FILE:',5A4) WRITE(IPRT,125)STEP 125 FORMAT(1X, ' STEP FOR NUMERICAL DERIVATIVE = ',E15.10) C C INDIRECT DISK INPUT INTO POTENTIAL ROUTINE C C READ(8,401)NPHI READ(IRD1,401)NPHI 401 FORMAT(16I5) C READ(8,403)R1S,R2S,VSP,XBC,XAB,CBC,ABC,CAB,AAB READ(IRD1,403)R1S,R2S,VSP,XBC,XAB,CBC,ABC,CAB,AAB 403 FORMAT(3E20.7) C READ(8,407)(PHI(I),I = 1,NPHI) READ(IRD1,407)(PHI(I),I = 1,NPHI) 407 FORMAT(4E20.7) DO 402 J = 1,3 C READ(8,407)(BSP(I,J),I = 1,NPHI) READ(IRD1,407)(BSP(I,J),I = 1,NPHI) 402 CONTINUE WRITE(IPRT,101) 101 FORMAT(1X, ' POTE(NTIAL INFO FROM THE DISK:'/) WRITE(IPRT,105)R1S,R2S,VSP 105 FORMAT(1X, ' R1S,R2S,VSP = ',3F10.5) WRITE(IPRT,111)(XBC(I),I= 1,3) 111 FORMAT(1X, ' BC: DE,RE,BE = ',10X,3F10.5) WRITE(IPRT,113)(PHI(I),(BSP(I,J),J = 1,3),I = 1,NPHI) 113 FORMAT(1X, ' PHI,DE,RE,BE = ',4F10.5) WRITE(IPRT,115)(XAB(I),I = 1,3) 115 FORMAT(1X, ' AB: DE,RE,BE = ',10X,3F10.5) WRITE(IPRT,117)(CBC(I),ABC(I),I = 1,3),(CAB(I),AAB(I),I = 1,3) 117 FORMAT(1X, ' FOR ASYMPTOTE, MO PARAMETER VALUE = (ASYMP. VALUE) + C*E 1XP(-A*(RI-RIS))'/10X,'(CBC(I),ABC(I),I = 1,3) = ',6F15.5/ 1 10X,'(CAB(I),AAB(I),I = 1,3) = ',6F15.5) RETURN C C COCK SPLINES C ENTRY POT2 IC=1 10 GO TO (81,82,83,84),IC 81 R1=R(1) R2=R(2) R3=R(3) GO TO 85 82 R1=R(1) + STEP R2=R(2) R3=R(3) GO TO 85 83 R1=R(1) R2=R(2) + STEP R3=R(3) GO TO 85 84 R1=R(1) R2=R(2) R3=R(3) + STEP 85 R1=R1*0.529177D0 R2=R2*0.529177D0 R3=R3*0.529177D0 DR1=R1S - R1 DR2=R2S - R2 ANGLE=DATAN(DR1/DR2) ANGLE=ANGLE/DETRAD DO 500 K=1,3 PPP=SPLINE(1,NPHI,PHI,BSP(1,K),CSP(1,K),SCRAP,ANGLE) 500 CONTINUE C C SETUP ANY CONSTANTS C DBC = XBC(1) DERG = XBC(1)/23.061D0 RERG = (R2S-XBC(2))/.529177D0 BETRG = XBC(3)*.529177D0 DEPR = XAB(1)/23.061D0 REPR = (R1S-XAB(2))/.529177D0 BETPR = XAB(3)*.529177D0 C C ENTRY INTO POTENTIAL CALCULATOR C IF(R1.GT.R1S) GO TO 100 IF(R2.GT.R2S) GO TO 102 C C INTERACTION REGION OF POLAR COORDINATES C DR1=R1S-R1 DR2=R2S-R2 ANGLE=DATAN(DR1/DR2) ANGLE=ANGLE/DETRAD RR=DSQRT(DR1*DR1+DR2*DR2) DIS=SPLINE(2,NPHI,PHI,BSP(1,1),CSP(1,1),SCRAP,ANGLE) REQ=SPLINE(2,NPHI,PHI,BSP(1,2),CSP(1,2),SCRAP,ANGLE) BET=SPLINE(2,NPHI,PHI,BSP(1,3),CSP(1,3),SCRAP,ANGLE) C DEDR(1)=SPLINE(3,NPHI,PHI,BSP(1,1),CSP(1,1),SCRAP,ANGLE) C DEDR(2)=SPLINE(3,NPHI,PHI,BSP(1,2),CSP(1,2),SCRAP,ANGLE) C30 DEDR(1)=DEDR(1)/627.5095D0*0.529177D0 C DEDR(2)=DEDR(2)/627.5095D0*0.529177D0 30 X = DIS*((1.D0-DEXP(BET*(RR-REQ)))**2.0D0 -1.D0) + VSP 12AUG84 V(IC) = X IC=IC+1 IF(IC.LT.5) GO TO 10 DO 333 I=1,3 333 V(I)=V(I)/627.5095D0 DEDR(1)=(V(2)-V(1))/STEP DEDR(2)=(V(3)-V(1))/STEP DEDR(3)=0.0D0 ENERGY=V(1) 93 RETURN C C LARGE R1 ASYMPTOTIC REGION C 100 IF(R2.GT.R2S) GO TO 104 RR = R2S-R2 SEP = R1-R1S DIS = XBC(1) + CBC(1)*DEXP(-DMIN1(ABC(1)*SEP,25.D0)) REQ = XBC(2) + CBC(2)*DEXP(-DMIN1(ABC(2)*SEP,25.D0)) BET = XBC(3) + CBC(3)*DEXP(-DMIN1(ABC(3)*SEP,25.D0)) C IF(R1.GT.10.D0)DEDR(1)=1.0D1 C IF(R1.GT.10.D0)DEDR(2)=1.0D-8 CI IF(R1.GT.10.D0)GO TO 30 C DEDR(1)=-CBC(1)*ABC(1)*((1.D0-DEXP(-BET*(RR-REQ)))**2-1.D0) -DIS* C 12.D0*(1.D0-DEXP(-BET*(RR-REQ)))*DEXP(-BET*(RR-REQ))*((-CBC(3)*ABC( C 13)*DEXP(-ABC(3)*SEP)*(RR-REQ) - CBC(2)*ABC(2)*DEXP(ABC(2)*SEP)* C 1BET))*BET*(RR-REQ) C DEDR(2)=2.D0*DIS*(1.D0-DEXP(-BET*(RR-REQ)))*BET*(RR-REQ)*DEXP(-BET* C 1(RR-REQ))*BET GO TO 30 C C LARGE R2 ASYMPTOTIC REGION C 102 CONTINUE RR = R1S-R1 SEP = R2-R2S DIS = XAB(1) + CAB(1)*DEXP(-DMIN1(AAB(1)*SEP,25.D0)) REQ = XAB(2) + CAB(2)*DEXP(-DMIN1(AAB(2)*SEP,25.D0)) BET = XAB(3) + CAB(3)*DEXP(-DMIN1(AAB(3)*SEP,25.D0)) C IF(R2.GT.10.D0)DEDR(2)=1.0D1 C IF(R2.GT.10.D0)DEDR(1)=1.0D-8 C IF(R2.GT.10.D0)GO TO 30 C DEDR(1)=-CAB(1)*AAB(1)*((1.D0-DEXP(-BET*(RR-REQ)))**2-1.D0) -DIS* C 12.D0*(1.D0-DEXP(-BET*(RR-REQ)))*DEXP(-BET*(RR-REQ))*((-CAB(3)*AAB( C 13)*DEXP(-AAC(3)*SEP)*(RR-REQ) - CAB(2)*AAB(2)*DEXP(AAB(2)*SEP)* C 1BET))*BET*(RR-REQ) C DEDR(2)=2.D0*DIS*(1.D0-DEXP(-BET*(RR-REQ)))*BET*(RR-REQ)*DEXP(-BET* C 1(RR-REQ))*BET GO TO 30 C C LARGE R1,R2 3 BODY BREAK-UP REGION C 104 CONTINUE J=J+1 ENERGY=XBC(1)/627.5095D0 DEDR(1)=1.0D-9 DEDR(2)=1.0D-9 DEDR(3)=0.0D0 GO TO 93 END C FUNCTION SPLINE(ISW,NN,X,Y,C,D,XPT) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(1),Y(1),C(1),D(1) C C THIS IS A SUBROUTINE FOR FITTING DATA WITH A CUBIC SPLINE C POLYNOMIAL AND EVALUATING THAT POLYNOMIAL AT A GIVEN POINT C OR ITS DERIVATIVE AT A GIVEN POINT C C CALLING SEQUENCE ....... C C ISW ... CONTROL OPTION C ISW=1 IF A CUBIC SPLINE IS TO BE FITTED TO THE SET OF KNOTS C DEFINED BY THE ARRAYS X AND Y. THE SPLINE COEFFICIENTS C ARE STORED IN THE ARRAY C. C ISW=2 IF THE SPLINE DEFINED BY THE COEFFICIENT ARRAY 'C' IS C TO BE EVALUATED (INTERPOLATED) AT THE POINT DEFINED BY C THE PARAMETER 'XPT'. C ISW=3 AS IN ISW=2, ONLY THE DERIVATIVE/3.D0 IS ALSO CALCULATED AT XPT C ISW=4 THE DERIVATIVE CALCULATED BY THE LAST USE OF SPLINE WITH ISW=3 C IS RETURNED. C C NN ... THE NUMBER OF KNOTS (DATA POINTS) TO WHICH THE SPLINE IS TO C BE FITTED C C X,Y ... THE ARRAYS DEFINING THE KNOTS. THE X-VALUES MUST BE IN C INCREASING ORDER. THE ARRAYS MUST BE DIMENSIONED AT LEAST C NN. C C C ... THE ARRAY THAT CONTAINS THE CUBIC SPLINE COEFFICIENTS. C MUST BE DIMENSIONED AT LEAST NN+2 . C C D ... A WORK SPACE. MUST BE DIMENSIONED AT LEAST NN+2 . C C XPT ... THE POINT AT WHICH THE INTERPOLATION IS DESIRED (IF ISW IS C SET TO 2). THE VALUE OF SPLINE IS SET TO THE C INTERPOLATED VALUE. C C ***** USER NOTES ***** C C INTERPOLATION INVOLVES AT LEAST TWO STEPS ....... C C A. CALL SPLINE WITH THE KNOTS. THIS SETS UP THE C COEFFICIENT ARRAY C. C EG. DUMY=SPLINE(1,NN,X,Y,C,D,XPT) C C B. CALL SPLINE WITH THE ARRAY C WHICH WAS DEFINED BY THE C PREVIOUS CALL AND WILL BE USED TO FIND THE VALUE AT THE C POINT 'XPT' . C EG. VALUE=SPLINE(2,NN,X,Y,C,D,XPT) C C STEP 'A' NEED BE EXECUTED ONLY ONCE FOR A GIVEN SET OF KNOTS. C STEP B MAY BE EXECUTED AS MANY TIMES AS NECESSARY. C 2 N=NN NP1=N+1 NP2=N+2 Z=XPT 24 GO TO (4,5,6,7),ISW 4 C(1)=Y(1) D(1)=1.0D0 C(NP1)=0.0D0 D(NP1)=0.0D0 C(NP2)=0.0D0 D(NP2)=0.0D0 DO 41 I=2,N C(I)=Y(I)-Y(1) 41 D(I)=X(I)-X(1) DO 410 I=3,NP2 IF(D(I-1).NE.0)GO TO 43 WRITE(IPRT,1001) CALL EXIT 43 PIVOT=1.0D0/D(I-1) IF(I.GE.NP2)GO TO 45 SUPD=X(I-1)-X(I-2) IF(SUPD.GE.0)GO TO 44 WRITE(IPRT,1000) CALL EXIT 44 SUPD=SUPD*SUPD*SUPD GO TO 46 45 SUPD=1.0D0 46 DFACT=SUPD*PIVOT CFACT=C(I-1)*PIVOT IF(I.GT.N)GO TO 48 DO 47 J=I,N V=X(J)-X(I-2) C(J)=C(J)-D(J)*CFACT 47 D(J)=V*V*V-D(J)*DFACT 48 CONTINUE IF(I.GE.NP2)GO TO 49 C(NP1)=C(NP1)-D(NP1)*CFACT D(NP1)=1.0D0-D(NP1)*DFACT 49 C(NP2)=C(NP2)-D(NP2)*CFACT 410 D(NP2)=X(I-2)-D(NP2)*DFACT DO 411 I=1,N J=NP2-I IF(J.NE.NP1)GO TO 413 V=1.0D0 GO TO 414 413 V=X(J)-X(J-1) V=V*V*V 414 IF(D(J+1).NE.0)GO TO 415 WRITE(IPRT,1001) CALL EXIT 415 C(J+1)=C(J+1)/D(J+1) 411 C(J)=C(J)-C(J+1)*V IF(D(2).NE.0)GO TO 416 WRITE(IPRT,1001) CALL EXIT 416 C(2)=C(2)/D(2) RETURN 5 SPLINE=C(1)+C(2)*(Z-X(1)) DO 51 I=1,N V=Z-X(I) IF(V.GT.0)GO TO 51 RETURN 51 SPLINE=SPLINE+C(I+2)*V*V*V RETURN 6 CONTINUE SPLINE=C(1)+C(2)*(Z-X(1)) DERIV = C(2)/3.D0 DO 53 I = 1,N V=Z-X(I) IF(V.LE.0) RETURN V2 = V*V SPLINE = SPLINE + C(I+2)*V2*V DERIV = DERIV + C(I+2)*V2 53 CONTINUE RETURN 7 CONTINUE SPLINE = 3.D0*DERIV RETURN 1000 FORMAT(1X,5X,'***** ERROR IN SPLINE ... UNORDERED X-VALUES 1 *****') 1001 FORMAT(1X,5X,' ***** ERROR IN SPLINE ... DIVIDE FAULT *****') END