SUBROUTINE POTENL(ICNTRL, MXLMB, MPOTL, LAM, R, P, ITYP) C C THIS IS A REPLACEMENT ROUTINE FOR SUBROUTINE POTENL IN MOLSCAT. C IT IS DESIGNED TO HANDLE CASES WHERE THE POTENTIAL ANGULAR C EXPANSION TERMS, V_SUB_ISYM, ARE GIVEN BY A TABLE OF VALUES C AT RADIAL DISTANCES R(I)=RMIN+(I-1)*DR, I=1,NPT. C N.B. THE SAME, EQUALLY SPACED RADIAL GRID MUST BE AVAILABLE C FOR ALL THE ANGULAR SYMETRIES. C DATA ARE READ FROM UNIT IIN (DEFAULT=12) USING FORMAT C STATEMENTS 501, 502, 503. SEE READ AND FORMAT STATEMENTS C FOR A DESCRIPTION OF THE INPUT. C CONTINUOUS RADIAL FUNCTIONS ARE OBTAINED FROM THE INPUT POINTS C BY CUBIC SPLINE INTERPOLATION, EXPONENTIAL EXTRAPOLATION AT C SMALL DISTANCES, AND INVERSE POWER EXTRAP AT LARGE DISTANCES. C (S. GREEN, J CHEM PHYS 67, 715 (1977). C C ************************************************************** C ** THIS DECK COMPATIBLE WITH MOLSCAT VERSION 14 (FEB 94) ** C ** USES IV INDEXING W/ ITYPE=2 ** C ** USES UPPER /MEMORY/ FOR STORAGE. IT IS POSSIBLE TO ** C ** DESTROY LAM() BY OVERWRITING W/ Y()..SCR() BUT THIS ** C ** SHOULD BE CAUGHT BY CHKSTR AFTER RETURN FROM POTENL ** C ** NB MANY STANDARD &POTL NAMELIST ITEMS HAVE BEEN REMOVED ** C ** DESCRIPTION OF VARIABLES FOLLOWS NAMELIST STATEMENT ** C ************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE IXMEM,RMINSV,DXSV,NPSV,RLSV,RUSV C PARAMETER (IIN=12) PARAMETER (MXLIN=5) C DIMENSION P(MXLMB), LAM(MXLMB) DIMENSION LAMIN(MXLIN) DIMENSION SP(3) INTEGER CFLAG CHARACTER*8 QNAME(10), QTYPE(10) CHARACTER*8 LABEL(10) EQUIVALENCE (MXLAM,MXSYM) COMMON/NPOT/NPTL C COMMON /MEMORY/ MX,IXNEXT,NIPR,IVLFL,X(1) C DATA MXSKPR/50/ DATA QTYPE/'LAMBDA =','ABS(MU)=',' MU = ',' L1 = ', 1 ' L2 = ',' L = ',' V = ','V-PRIME=', 2 ' J = ','J-PRIME='/ C NAMELIST/POTL/RM,EPSIL, MXLAM,MXSYM, NPOTL,CFLAG, 1 LAMMAX,MXV,MXDELV, IPRINT C C RM - LENGTH SCALING FACTOR C EPSIL - ENERGY SCALING FACTOR C MXLAM - NUMBER OF POTENTIAL TERMS RETURNED C MXSYM - SYNONYM FOR MXLAM (EQUIVALENCED VARIALBES) C NPOTL - NUMBER OF ELEMENTS OF VL ARRAY (IN BASE) WHICH C CORRESPOND TO EACH ELEMENT OF P ARRAY C CFLAG - FLAG FOR SCALING POTENTIAL FOR ITYPE = 5 OR 6 C IPRINT - CONTROLS OUTPUT OF POTENTIAL INFORMATION C =0 (MINIMAL), =1 (SUMMARY), =2 (FULL) C ------ VARIABLES TO CONTROL RETAINING RECORDS FROM FILE(IIN) C LAMMAX - HIGHEST LEGENDRE TERM TO KEEP C MXV - HIGHEST VIBRATIONAL LEVEL TO KEEP (ITYPE=2) C MXDELV - HIGHEST DELTA-V TO KEEP (ITYPE=2) C IF (ICNTRL.GE.0.AND.ICNTRL.LE.2) GOTO 1000 IF (ICNTRL.EQ.-1) GOTO 2000 WRITE(6,631) ICNTRL,R 631 FORMAT('0 *** ERROR IN POTENL, ICNTRL =',I6,' R =',E16.8) STOP C C EVALUATE V(R) C 1000 NSTOR=2*NPSV+4 IXSC=IXMEM-NPSV IF (R.LT.RLSV) GO TO 5911 C IF (R.GT.RUSV) GO TO 5912 - BUG FIXED 3 MAR 95, SG IF (R.GE.RUSV) GO TO 5912 C C BELOW FOR CASE THAT POTENTIAL IS INTERPOLATED. XMIN=RMINSV DR=DXSV N=NPSV C FIND I SUCH THAT R IS IN INTERVAL X(I) TO X(I+1). Q=(R-XMIN)/DR IF (Q.LT.0D0) GO TO 9999 I=Q I=I+1 IF (I.GE.N) GO TO 9999 HT1=R-XMIN-(I-1)*DR HT2=HT1-DR Z=HT1*HT2 C LOOP OVER SYMMETRIES DO 1111 J=1,MXLMB IXY=IXSC-(NPSV+4) SSS=(X(IXSC+I+1)-X(IXSC+I))/DR C CALCULATE 2ND DERIVATIVE AT T. SP(3)=X(IXSC+I)+HT1*SSS C CALCULATE FUNCTION AT T. DELSQS=(X(IXSC+I)+X(IXSC+I+1)+SP(3))/6D0 DELY=(X(IXY+I+1)-X(IXY+I))/DR SP(1)=X(IXY+I)+HT1*DELY +Z*DELSQS C FIND 1ST DERIVATIVE AT T. SP(2)=DELY+(HT1+HT2)*DELSQS+Z*SSS/6D0 P(J)=SP(ICNTRL+1) 1111 IXSC=IXSC-NSTOR RETURN 9999 WRITE(6,692) R 692 FORMAT('0 *** POTENL(SPLINE) ERROR INTERP REQUESTED FOR ' 1 ,'ILLEGAL R =',D16.8) STOP C C EXTRAPOLATION AT LEFT 5911 DO 5997 I=1,MXLMB TERM=0.D0 AA=X(IXSC-3) B=X(IXSC-2) IF (AA .EQ.0.D0 .OR. B .EQ.0.D0) GO TO 5897 TERM=AA *EXP(-B *R) IF (ICNTRL.GT.0) TERM=TERM*(-B)**ICNTRL 5897 P(I)=TERM 5997 IXSC=IXSC-NSTOR RETURN C C EXTRAPOLATE AT RIGHT 5912 DO 5996 I=1,MXLMB C=X(IXSC-1) D=X(IXSC) TERM=0.D0 IF (C .EQ.0.D0 .OR. D .EQ.0.D0) GO TO 5896 TERM=C *R**(-D ) IF(ICNTRL.EQ.1) TERM=-D *TERM/R IF(ICNTRL.EQ.2) TERM=D *(D +1.D0)*TERM/(R*R) 5896 P(I)=TERM 5996 IXSC=IXSC-NSTOR RETURN C C POTENTIAL INITIALISATION C 2000 WRITE(6,630) IIN 630 FORMAT('0 POTENL(SPLINE) ROUTINE (FEB 94).'/ 1 '0 TABULATED V(ISYM) WILL BE INPUT FROM UNIT',I4) C RM=0.D0 EPSIL=0.D0 NPOTL=-1 MXLAM=-1 LAMMAX=-1 MXV=-1 MXDELV=-1 CFLAG=0 IPRINT=1 READ(5,POTL) C C REPORT &POTL INPUT AND CHECK CONSISTENCY W/ POTENL(SPLINE) C ITYPE=ITYP-10*(ITYP/10) WRITE(6,634) IPRINT 634 FORMAT('0 PRINT LEVEL FROM &POTL IPRINT =',I3) IF(NPOTL.NE.-1) THEN WRITE(6,*) ' *** POTENL(SPLINE) IGNORING INPUT NPOTL',NPOTL ENDIF IF (RM.NE.0.D0.OR.EPSIL.NE.0.D0) THEN WRITE(6,*) ' *** POTENL(SPLINE) INPUT RM, EPSIL =',RM,EPSIL WRITE(6,*) ' WILL BE IGNORED IN FAVOR OF VALUES ON (IIN).' WRITE(6,*) ' *** BE SURE THIS IS CORRECT.' ENDIF IF (MXLAM.GT.0) THEN WRITE(6,*) ' *** POTENL(SPLINE) INPUT MXLAM =',MXLAM WRITE(6,*) ' MAY LIMIT NUMBER OF SYMMETRIES READ' ENDIF IF (LAMMAX.GE.0) THEN WRITE(6,*) ' *** POTENL(SPLINE) INPUT LAMMAX =',LAMMAX WRITE(6,*) ' MAY LIMIT NUMBER OF SYMMETRIES RETAINED' ENDIF IF (MXV.GE.0.AND.(ITYPE.EQ.2.OR.ITYPE.EQ.7)) THEN WRITE(6,*) ' *** POTENL(SPLINE) INPUT MXV =',MXV WRITE(6,*) ' MAY LIMIT NUMBER OF SYMMETRIES RETAINED' ENDIF IF (MXDELV.GE.0.AND.(ITYPE.EQ.2.OR.ITYPE.EQ.7)) THEN WRITE(6,*) ' *** POTENL(SPLINE) INPUT MXDELV =',MXDELV WRITE(6,*) ' MAY LIMIT NUMBER OF SYMMETRIES RETAINED' ENDIF C C GET AND CHECK HEADER INFORMATION FROM UNIT(IIN) C MXLTP IS NO. SYMMETRIES ON TAPE; NQPLTP IS NO. LABELS/SYMMETRY C EACH SYMMETRY DESCRIBED BY NP VALUES FROM XMIN IN STEPS OF DR C RM, EPSIL ARE LENGTH AND ENERGY SCALING VALUES FOR MOLSCAT C READ(IIN,501,END=9100) LABEL,MXLTP,NQPLTP,NP,XMIN,DR,RM,EPSIL 501 FORMAT(10A8/3I5,4F10.4) WRITE(6,632) LABEL 632 FORMAT(/' LABEL FROM POTENTIAL FILE --'/1X,10A8) NPSV=NP RMINSV=XMIN DXSV=DR RLSV=XMIN RUSV=XMIN+(NP-1)*DR IF (MXLAM.GE.1.AND.MXLAM.LT.MXLTP) THEN WRITE(6,*) ' *** POTENL(SPLINE). INPUT MXLAM =',MXLAM WRITE(6,*) ' IS LESS THAN VALUE FROM TAPE',MXLTP WRITE(6,*) ' AND WILL LIMIT NUMBER OF SYMMETRIES INPUT' MXLTP=MXLAM ENDIF IF (IPRINT.LT.2) WRITE(6,633) NPSV,RMINSV,DXSV,RUSV 633 FORMAT(/' INTERPOLATION ON',I4,' POINTS FROM',F8.3, 1 ' IN STEPS OF',F7.3,' TO',F8.3) IF (NPSV.LT.2) THEN WRITE(6,*) ' *** POTENL(SPLINE). ERROR. TOO FEW POINTS' 1 ,' FOR SPLINE',NPSV STOP ENDIF C C ATTEMPT TO PROCESS ITYPE AND POTENTIAL DESCRIPTION NUMBERS C IF (IPRINT.GE.1) WRITE(6,639) 639 FORMAT(/' ANGULAR DEPENDENCE OF POTENTIAL EXPANDED IN TERMS OF') IF(ITYPE.EQ.1) GOTO 2001 IF(ITYPE.EQ.2) GOTO 2002 IF(ITYPE.EQ.3) GOTO 2003 IF(ITYPE.EQ.5) GOTO 2005 IF(ITYPE.EQ.6) GOTO 2005 IF(ITYPE.EQ.7) GOTO 2002 C WRITE(6,640) ITYPE 640 FORMAT(' *** POTENL(SPLINE). ITYPE =',I4,' IS NOT SUPPORTED.') STOP C 2001 NQPL=1 QNAME(1)=QTYPE(1) IF (IPRINT.GE.1) WRITE(6,641) 641 FORMAT(' LEGENDRE POLYNOMIALS, P(LAMBDA).') GOTO 2100 C 2002 NQPL=3 QNAME(1)=QTYPE(1) QNAME(2)=QTYPE(7) QNAME(3)=QTYPE(8) IV1=2 IV2=3 IF (IPRINT.GE.1) THEN WRITE(6,641) WRITE(6,642) ENDIF 642 FORMAT(' INTEGRATED OVER DIATOM VIBRATIONAL FUNCTIONS') IF(ITYPE.NE.7) GOTO 2100 NQPL=5 QNAME(3)=QTYPE(9) QNAME(4)=QTYPE(8) QNAME(5)=QTYPE(10) IV2=4 IF (IPRINT.GE.1) WRITE(6,643) 643 FORMAT(' FOR EACH PAIR OF V,J LEVELS') GOTO 2100 C 2003 NQPL=3 QNAME(1)=QTYPE(4) QNAME(2)=QTYPE(5) QNAME(3)=QTYPE(6) IF (IPRINT.GE.1) WRITE(6,644) 644 FORMAT(' CONTRACTED NORMALISED SPHERICAL HARMONICS, SUM', 1 '(M1,M2,M) C(L1,M1,L2,M2,L,M) Y(L1,M1) Y(L2,M2) Y(L,M)'/ 2 ' SEE RABITZ, J. CHEM. PHYS. 57, 1718 (1972)') GOTO 2100 C 2005 NQPL=2 QNAME(1)=QTYPE(1) QNAME(2)=QTYPE(2) IF (IPRINT.GE.1) WRITE(6,645) 645 FORMAT(' NORMALISED SPHERICAL HARMONICS: (Y(LAM,MU) + ', 1 '(-)**MU Y(LAM,-MU)) / (1+DELTA(MU,0))') IF(CFLAG.EQ.1) WRITE(6,646) 646 FORMAT(' *** POTENL(SPLINE) CFLAT=1 REQUEST TO SCALE POTENTIAL', 1 ' IGNORED.'/' *** MAKE SURE THIS IS OKEY.') GOTO 2100 C 2100 IF (NQPLTP.NE.NQPL) THEN WRITE(6,*) ' *** POTENL(SPLINE) - POSSIBLE ERROR' WRITE(6,*) ' NQPL APPROPRIATE TO ITYPE',ITYPE WRITE(6,*) ' NOT EQUAL TO VALUE ON TAPE',NQPLTP WRITE(6,*) ' LATTER WILL BE USED' NQPL=NQPLTP IF (NQPL.GT.MXLIN) THEN WRITE(6,*) ' *** ERROR. EXCEEDS DIMENSION MXLIN',MXLIN STOP ENDIF ENDIF C C GET INFORMATION FOR EACH SYMMETRY FROM UNIT(IIN) IQ=0 MXLAM=0 NSKIP=0 C PREPARE TO USE /MEMORY/...,X -- ALLOCATE SPACE FOR X(NPSV) IXMEM=MX MX=MX-NPSV DO 9000 I=1,MXLTP READ(IIN,502,END=9100) (LAMIN(IX),IX=1,NQPL) 502 FORMAT(10I5) C SEE IF THIS SET MEETS ALL REQUIRMENTS OF LAMMAX,MXV,MXDELV IF (LAMMAX.GE.0.AND.LAMIN(1).GT.LAMMAX) THEN NSKIP=NSKIP+1 IF (IPRINT.GE.1.AND.NSKIP.LE.MXSKPR) THEN WRITE(6,*) ' *** POTENL(SPLINE). SKIPPING INPUT SET BECAUSE', 1 ' OF LAMMAX',LAMMAX WRITE(6,*) (LAMIN(IX),IX=1,NQPL) ENDIF GO TO 9001 ENDIF IF (ITYPE.EQ.2.OR.ITYPE.EQ.7) THEN IF (MXV.GE.0.AND.(LAMIN(IV1).GT.MXV.OR.LAMIN(IV2).GT.MXV)) THEN NSKIP=NSKIP+1 IF (IPRINT.GE.1.AND.NSKIP.LE.MXSKPR) THEN WRITE(6,*) ' *** POTENL(SPLINE). SKIPPING INPUT SET', 1 ' BECAUSE OF MXV',MXV WRITE(6,*) (LAMIN(IX),IX=1,NQPL) ENDIF GO TO 9001 ENDIF IF (MXDELV.GE.0.AND.ABS(LAMIN(IV1)-LAMIN(IV2)).GT.MXDELV) THEN NSKIP=NSKIP+1 IF (IPRINT.GE.1.AND.NSKIP.LE.MXSKPR) THEN WRITE(6,*) ' *** POTENL(SPLINE). SKIPPING INPUT SET', 1 ' BECAUSE OF MXDELV',MXDELV WRITE(6,*) (LAMIN(IX),IX=1,NQPL) ENDIF GO TO 9001 ENDIF ENDIF C C IF WE REACH HERE, RETAIN THE INPUT SYMMETRY. CHECK STORAGE MXLAM=MXLAM+1 IF (MXLAM*NQPL.GT.MXLMB) THEN WRITE(6,*) ' *** POTENL(SPLINE). INADEQUATE STORAGE IN EXTERNAL' 1 ,' LAM ARRAY' WRITE(6,*) ' MXLAM, NQPL, MXLMB =',MXLAM,NQPL,MXLMB STOP ENDIF C STORAGE FOR Y(NPSV),AA,B,C,D,SCR(NPSV) NSTOR=2*NPSV+4 MX=MX-NSTOR IF (MX.LT.IXNEXT) THEN WRITE(6,*) ' *** POTENL(SPLINE) INSUFFICIENT SPACE IN /MEMORY/' WRITE(6,*) ' TRYING TO PROCESS SYMMETRY',MXLAM WRITE(6,*) ' ORIGINAL/CURRENT MX, IXNEXT =',IXMEM,MX,IXNEXT STOP ENDIF IXX=MX IXY=MX+NPSV IXSC=IXY+NPSV+4 C C PUT LAMIN() VALUES INTO LAM() AND SET UP, IE, INPUT Y(,MXLAM) DO 9010 J=1,NQPL 9010 LAM(IQ+J)=LAMIN(J) IF (IPRINT.GE.1) THEN WRITE(6,651) MXLAM 651 FORMAT('0 INTERACTION POTENTIAL FOR SYMMETRY TYPE NUMBER',I3) WRITE(6,652) (QNAME(J),LAM(IQ+J),J=1,NQPL) 652 FORMAT(' WHICH HAS ',6(A8,I3,3X)) ENDIF READ(IIN,503,END=9100) (X(IXY+J),J=1,NPSV) 503 FORMAT(4D20.12) DO 1011 J=1,NPSV 1011 X(IXX+J)=RMINSV+(J-1)*DXSV IF (IPRINT.GE.2) THEN WRITE(6,654) 654 FORMAT(1X) WRITE(6,601) NPSV,(X(IXX+J),X(IXY+J),J=1,NPSV) 601 FORMAT('0',5X,'POTENTIAL WILL BE INTERPOLATED FROM FOLLOWING', 1 I5,' POINTS'/ ( 9X,4(F7.4,F15.6,8X))) ENDIF C C GET EXTRAPOLATION TO LEFT AS Y =A(I)*DEXP(-B(I)*R) Y1=X(IXY+1) Y2=X(IXY+2) X1=X(IXX+1) X2=X(IXX+2) IF (Y1*Y2.GT.0.) GO TO 1300 IF (IPRINT.GE.1) WRITE(6,697) 697 FORMAT('0 * * * WARNING. SHORT RANGE CANNOT BE FIT AS EXPONENTIAL &. SET TO ZERO.') AA =0. B =0. SCR1=0.D0 GO TO 1399 1300 B =DLOG(Y1/Y2)/(X2-X1) AA =Y1*DEXP(B*X1) SCR1=AA*B*B*DEXP(-B*X1) 1399 IF (IPRINT.GE.2) WRITE(6,602) RLSV,AA ,B 602 FORMAT('0 FOR R LESS THAN',F7.3,', V =',1P,D12.4,' * DEXP( -', 1 0P,F8.3,' * R ).') C C AND THEN EXTRAPOLATION TO RIGHT AS C*R**(-D) YN=X(IXY+NPSV) YNM1=X(IXY+NPSV-1) XN=X(IXX+NPSV) XNM1=X(IXX+NPSV-1) IF (YNM1*YN.LE.0D0) GOTO 1299 IF (YNM1/YN.GT.1D0) GO TO 1200 1299 IF (IPRINT.GE.1) WRITE(6,698) 698 FORMAT('0 * * * WARNING. LONG-RANGE CANNOT BE FIT AS DECREASING E 1XPONENTIAL. SET TO ZERO.') D =0D0 C =0D0 SCRN=0.D0 X(IXY+NPSV)=0.D0 GO TO 1100 1200 D = 1 DLOG(YNM1/YN)/DLOG(XN/XNM1) C ELIMINATE LONG=RANGE THAT DECREASES EXTREMELY FAST IF (D .GT.8.D1) GO TO 1299 C =YN*XN**D SCRN =D *(D +1.D0)* 1 C *XN**(-D -2.D0) C WEED OUT POSSIBLY SPURIOUS FITS WHICH FALL OFF TOO SLOWLY. IF (D .GT.4.D0) GO TO 1100 IF (IPRINT.GE.1) WRITE(6,686) D 686 FORMAT('0 * * * WARNING. INVERSE POWER =',D14.6,' IS UNPHYSICAL. 1 LONG RANGE FORCED TO ZERO UNLESS SPLNEQ CHANGED.') GO TO 1299 1100 IF (IPRINT.GE.2) WRITE(6,603) RUSV,C ,D 603 FORMAT('0 FOR R GREATER THAN',F7.3,', V =',1P,D12.4,' / R**(', 1 0P,F8.3,' ).') C BUT AA,B,C,D,SCR1,SCRN INTO APPROPRIATE STORAGE X(IXSC-3)=AA X(IXSC-2)=B X(IXSC-1)=C X(IXSC)=D X(IXSC+1)=SCR1 X(IXSC+NPSV)=SCRN C INITIALIZE SCR(,I) FOR SPLINE CALL SPLNXX(RMINSV,DXSV,NPSV,X(IXY+1),X(IXSC+1)) IQ=IQ+NQPL GO TO 9000 C SKIP OVER POINTS IF WE ARE NOT RETAINING THIS SYMMETRY 9001 READ(IIN,503,END=9100) (X(IXX+IX),IX=1,NP) C 9000 CONTINUE C IF (NSKIP.GT.MXSKPR.OR.(NSKIP.GT.0.AND.IPRINT.LE.0)) 1 WRITE(6,604) NSKIP 604 FORMAT('0 *** TOTAL OF',I6,' SYMMETRIES WERE SKIPPED ON &POTL', 1 ' INPUT CRITERIA.') C FREE UP STORAGE USED FOR X(NPSV) AND NO LONGER NEEDED MX=MX+NPSV C SET PARAMETERS INTO CALLING ARGUMENTS NPOTL=MXLAM C BELOW FOR NEW (FEB 94) USE OF IV WITH BOTH ITYPE=2 AND 7 IF (ITYPE.EQ.2.OR.ITYPE.EQ.7) THEN NPOTL=0 DO 2010 I=1,MXLAM 2010 NPOTL=MAX0(NPOTL,LAM((I-1)*NQPL+1)) NPOTL=NPOTL+1 ENDIF C WRITE(6,661) EPSIL,RM,MXLAM,NPOTL 661 FORMAT('0 ENERGY IN UNITS OF EPSILON =',F15.5,' CM-1'/ 1 ' R IN UNITS OF RM =',F15.5,' ANGSTROMS'/ 2 '0 MXLAM =',I5/' NPOTL =',I5) C R=RM P(1)=EPSIL MPOTL=NPOTL MXLMB=MXLAM RETURN C C UNEXPECTED END OF FILE ON (IIN) 9100 WRITE(6,*) ' *** POTENL(SPLINE). UNEXPECTED EOF ON INPUT TAPE.' WRITE(6,*) ' *** TERMINATING.' STOP END SUBROUTINE SPLNXX (XMIN,DR,N,Y,SCR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N), SCR(N) DATA MAXTRY/100/ C C DO INITIALIZATION FOR CUBIC SPLINE FIT. NM1=N-1 DR2=DR*DR H2=2D0*DR2 OMEGA=4D0*(2D0-DSQRT(3D0)) EPS=1D-8 C DEFINITION OF 'NATURAL' SPLINE HAS S2(1) = S2(N) = 0. C NOT USED HERE. ASSUME THAT SCR(1), SCR(N) SET IN CALLING PROG. C SCR(1)=0D0 C SCR(N)=0D0 C DO 52 I=2,NM1 52 SCR(I)=(Y(I+1)-2D0*Y(I)+Y(I-1))/DR2 DO 521 ICNT=1,MAXTRY ETA=0. DO 10 I=2,NM1 DELSQY=(Y(I+1)-2D0*Y(I)+Y(I-1))/H2 W=(3D0*DELSQY-.25D0*(SCR(I-1)+SCR(I+1))-SCR(I))*OMEGA SCR(I)=SCR(I)+W IF (SCR(I).NE.0D0) W=W/SCR(I) 10 ETA=DMAX1(ETA,DABS(W)) C IF (ETA.GE.EPS) GO TO 521 C SCR(), WHICH IS EQUAL TO 2ND DERIV. AT ABSCISSAE, HAS CONVERGED. C WRITE(6,600) ICNT,(SCR(I),I=1,N) 600 FORMAT('0 SPLINE INIT. CONVERGED IN',I4,' ITERATIONS.'/ 1 (' ',8D16.8)) RETURN 521 CONTINUE WRITE(6,601) MAXTRY 601 FORMAT('0 * * * ERROR. SPLINE INIT. DID NOT CONVERGE IN',I4 1 ,' ITERATIONS FOR 2ND DERIVATIVE AT KNOTS.') STOP END