c *****************************************************
      PROGRAM PULS
c******************************************************
c This program solves both the nonradial pulsation problem
c in the adiabatic Cowling approximation and the
c purely radial adiabatic problem. We thank
c W. D. Pesnell and P. Jones for help in writing earlier
c versions of this code which solved the full fourth
c order system for adiabatic nonradial oscillations
c (plus quite a few other things).
c The input used is generated by the ZAMS model builder.
c The method is to solve the two Cowling approximation
c equations in their Dziembowski form or the two
c purely radial equations by shooting from
c the center to the surface where, eventually and by
c iteration, the outer boundary conditions on the
c eigenfunctions are satisfied. The only guess you have
c to provide is that for the period (in seconds). You
c must also specify L which is angular momentum
c quantum number of the Spherical Harmonic Function.
c Enter a zero if you want the radial case.
c ******************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON G(500),RHO(500),X(500),YEIG(2,500)
      DOUBLE PRECISION L,LHAT,LINDEX
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      LOGICAL RADIAL
      COMMON/RAD/RADIAL
      COMMON/RS/R(500)
      COMMON/FLAGS/IPRNT
C
C Read in model quantities.
C
      CALL READ
c
   10 WRITE (6,2000)
 2000 FORMAT (' ENTER GUESSED PERIOD (IN SECONDS)')
      WRITE (6,2001)
 2001 FORMAT (' (ENTER A PERIOD OF 0 TO STOP)')
      READ (5,*) PERIOD
      IF (PERIOD.NE.0.0D0) GO TO 20
      WRITE (6,1000)
 1000 FORMAT (' CALCULATION COMPLETED')
      STOP
   20 WRITE (11,1010) PERIOD
 1010 FORMAT (' GUESSED PERIOD=',1PE14.6,4H SEC)
      EIG=(2.D0*PI/PERIOD)**2
c EIG is the first guess at the eigenvalue = sigma**2.
         EIG = (2.0D0*PI/PERIOD)**2
         NCONV = 0
         NSERCH=20
c
c Try to converge to a solution using Newton's method.
c
         DO 70 NTRY=1,NSERCH
            EIGT=EIG
            IF (RADIAL) THEN
               CALL GOOUT0(DISC)
            ELSE
               CALL GOOUTL(DISC)
            ENDIF
c
c Perturb eigenvalue and find change in boundary condition.
c DISC is the boundary condition.
c
            EIGT=(1.D0+EPS)*EIG
            IF (RADIAL) THEN
               CALL GOOUT0(DUM1)
            ELSE
               CALL GOOUTL(DUM1)
            ENDIF
     	    DEIGD=DUM1-DISC
            DEIG=-DISC*EPS*EIG/DEIGD
            WRITE (6,1060) NTRY,EIG,DEIG
            WRITE (11,1060) NTRY,EIG,DEIG
c DEIG is the predicted change in EIG.
 1060 FORMAT(' NTRY=',I3,' EIG=',1PE15.7,' DEIG=',
     *         1PE13.5)
            ADEIG=DABS(DEIG/EIG)
c Check for convergence.
            IF (ADEIG.LT.VERG) NCONV=1
c Compute new guess for EIG.
            EPSEIG=DEIG/EIG
            IF (ADEIG.GT.0.1D0)
     *          EPSEIG=0.1D0*DEIG/DABS(DEIG)
            EIG=EIG*(1.D0+EPSEIG)
            IF (NCONV.EQ.1) GO TO 80
   70    CONTINUE
         WRITE (11,1080)
 1080 FORMAT (14H NOT CONVERGED)
         STOP
c
c Converged.
c
   80    EIGT=EIG
            IF (RADIAL) THEN
               CALL GOOUT0(DISC)
            ELSE
               CALL GOOUTL(DISC)
            ENDIF
         PERIOD=(2.D0*PI)/DSQRT(EIG)
         WRITE ( 6,1090) PERIOD,EIG
         WRITE (11,1090) PERIOD,EIG
 1090 FORMAT (' FINAL VALUES: PERIOD= ',1PE14.6,
     *      ' EIG=',1PE17.9)
c Count nodes in Y1.
         CALL MODEID
         IF (RADIAL) THEN
            WRITE (6,1050) NODES1
            WRITE (11,1050) NODES1
 1050 FORMAT (I3,' RADIAL NODES IN Y1')
         ELSE
            WRITE (6,1110) NODES1,NODES2,MODEP
            WRITE (11,1110) NODES1,NODES2,MODEP
 1110 FORMAT(I3,' NODES IN Y1   ',I3,' NODES IN Y2 '/
     @'     PHASE DIAGRAM MODE ',I3)
         ENDIF
         IF (IPRNT.NE.1) GO TO 140
         YSURF=YEIG(1,NSURF)
         DO 600 I=1,NSURF
             YEIG(1,I)=YEIG(1,I)/YSURF
             YEIG(2,I)=YEIG(2,I)/YSURF
  600 CONTINUE
         WRITE (11,1130)
         WRITE (11,1140) (I+1,R(I)/R(NSURF),YEIG(1,I),
     *       YEIG(2,I),I=1,NSURF)
 1130 FORMAT('                   2(N, R/R* ,  Y1  ,  Y2)')
 1140 FORMAT(2(I4,0PF10.6,2(1PE12.4)))
  140 CONTINUE
      GO TO 10
      END
c *****************************************************
c *****************************************************
      SUBROUTINE READ
c  This routine reads in model quantities.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*1 IYORN
      COMMON G(500),RHO(500),X(500),YEIG(2,500)
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      COMMON/RS/R(500)
      COMMON/FLAGS/IPRNT
      COMMON /SAVEL/M,XI(500),C1(4,500),Y1(500),C2(4,500),
     @ Y2(500),C3(4,500),Y3(500),C4(4,500),Y4(500),
     @ C5(4,500),Y5(500),C6(4,500)
      COMMON /SURF/RSURF
      DOUBLE PRECISION L,LHAT,LINDEX
      CHARACTER*20 OUTFILE, INFILE
      LOGICAL RADIAL
      COMMON/RAD/RADIAL
c
      RADIAL=.FALSE.
      WRITE (6,*) ' YOUR INPUT FILENAME IS'
      READ (5,3000) INFILE
 3000 FORMAT (A20)
      OPEN (10,FILE=INFILE,STATUS='UNKNOWN')
      WRITE (6,*)  ' YOUR OUTPUT FILENAME IS '
      READ (5,3000) OUTFILE
      OPEN (11,FILE=OUTFILE,STATUS='UNKNOWN')
      WRITE (6,1080)
 1080 FORMAT (' ENTER L')
      READ (5,*) L
      LHAT=L*(L+1.D0)
      LINDEX=2.0D0-L
c Setting up convergence criteria and EPS for
c Newton's method.
      VERG=1.00D-05
      EPS=1.00D-03
      PI=3.1415926535898D0
      P43=4.D0*PI/3.D0
      PI4=4.D0*PI
      GRAV=6.6723D-08
      READ (10,*) AMASS
      READ (10,*) XX,YY
      WRITE (11,1000) AMASS
      WRITE (6,1000) AMASS
 1000 FORMAT (' MASS (IN MSUN)=',0PF7.3)
      WRITE (11,1001) XX,YY
      WRITE (6,1001) XX,YY
 1001 FORMAT (' X=',0PF6.3,' Y=', F6.3)
      WRITE (11,1090) L
 1090 FORMAT (3H L=,0PF7.0)
c Check to see if a purely radial calculation
c is to be done.
      IF (L .LT. 0.9D0) RADIAL=.TRUE.
c Read in model quantities from ZAMS.
      READ (10,*) NSURF
      READ (10,*) (X(I),R(I),G(I),RHO(I),I=1,NSURF)
      RSURF=R(NSURF)
      READ (10,*) M
      DO 100 I=1,M
           READ (10,*) XI(I),Y1(I),Y2(I),Y3(I)
           READ (10,*) Y4(I),Y5(I)
  100 CONTINUE
      IF (RADIAL) THEN
c For the radial case Y4 contains V (of the U-V plane)
c and Y3 contains Gamma1.
         DO 101 I=1,M
           Y4(I)=1.0D0/Y5(I)-1.0D0
           Y3(I)=Y4(I)/Y2(I)
  101    CONTINUE
      ENDIF
c Set up spline coefficients for use with integrator.
      CALL CONSL(XI,Y1,M,C1)
      CALL CONSL(XI,Y2,M,C2)
      CALL CONSL(XI,Y3,M,C3)
      CALL CONSL(XI,Y4,M,C4)
      CALL CONSL(XI,Y5,M,C5)
      CALL CONSL(XI,R ,M,C6)
c Do you want the eigenfunctions saved to file?
      WRITE (6,2000)
 2000 format (' PRINT OUT EIGENFUNCTIONS (Y/N)?')
      READ (5,1100) IYORN
 1100 FORMAT (A1)
      IF (IYORN.EQ.'Y'.OR.IYORN.EQ.'y') IPRNT=1
      RETURN
      END
c ******************************************************
c ******************************************************
      SUBROUTINE GOOUTL(B2)
c This controls the integration for the Cowling
c approximation.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON G(500),RHO(500),X(500),YEIG(2,500)
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      COMMON/RS/R(500)
      DOUBLE PRECISION L,LHAT,LINDEX
      COMMON/SPLINQ/VQ(6)
      DIMENSION Y(2),WORK(50),IWORK(6)
      EXTERNAL RKFCOW
c
      YEIG(1,1)=1.0D0
      YEIG(2,1)=EIGT/(L*GRAV*P43*RHO(1))
      Y(1)=YEIG(1,1)
      Y(2)=YEIG(2,1)
      IFLAG=1
      NEQN=2
      NSURF1=NSURF-1
   50 DO 20 J=1,NSURF1
         K=J+1
         XSTART=X(J)
         XFIN=X(K)
         RELERR=1.D-8
         ABSERR=DABS(Y(1))
      DO 30 I=2,NEQN
         TEMP=DABS(Y(I))
         ABSERR=DMIN1(ABSERR,TEMP)
   30 CONTINUE
         ABSERR=1.D-8*ABSERR
         ABSERR=DMAX1(ABSERR,1.D-20)
         CALL RKF(RKFCOW,NEQN,Y,XSTART,XFIN,RELERR,ABSERR,
     @        IFLAG,WORK,IWORK)
         IF (IFLAG.EQ.3.OR.IFLAG.EQ.4.OR.IFLAG.EQ.5)
     @              WRITE (11,1000) IFLAG,XFIN
         IFLAG=1
         DO 10 I=1,2
   10         YEIG(I,K)=Y(I)
   20 CONTINUE
      RT=X(NSURF)
      CALL SPLNTL(RT,VQ)
      TEMP=1.D0/VQ(5)-1.D0
c The following is the outer boundary condition. We
c iterate until B2 is close enough to zero that
c corrections to the eigenvalue are within a small
c tolerance.
      B2=YEIG(1,NSURF)*((4.D0+R(NSURF)*EIGT/G(NSURF))
     @    /TEMP-1.D0)+(1.D0-LHAT*G(NSURF)/R(NSURF)
     @   /EIGT/TEMP)*YEIG(2,NSURF)
 1000 FORMAT (17H WATCH OUT,IFLAG=,I4,6H XFIN=,1PE10.2)
      RETURN
      END
c ***************************************************
c ***************************************************
      SUBROUTINE GOOUT0(B2)
c This controls the integration for purely
c radial pulsations.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON G(500),RHO(500),X(500),YEIG(2,500)
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      COMMON/RS/R(500)
      DOUBLE PRECISION L,LHAT,LINDEX
      COMMON/SPLINQ/VQ(6)
      DIMENSION Y(2),WORK(50),IWORK(6)
      EXTERNAL RKFRAD
c
      YEIG(1,1)=1.0D0
      RT=X(1)
      CALL SPLNTL(RT,VQ)
      YEIG(2,1)=-3.0D0*VQ(3)*YEIG(1,1)
      Y(1)=YEIG(1,1)
      Y(2)=YEIG(2,1)
      IFLAG=1
      NEQN=2
      NSURF1=NSURF-1
   50 DO 20 J=1,NSURF1
         K=J+1
         XSTART=X(J)
         XFIN=X(K)
         RELERR=1.D-8
         ABSERR=DABS(Y(1))
      DO 30 I=2,NEQN
         TEMP=DABS(Y(I))
         ABSERR=DMIN1(ABSERR,TEMP)
   30 CONTINUE
         ABSERR=1.D-8*ABSERR
         ABSERR=DMAX1(ABSERR,1.D-20)
         CALL RKF(RKFRAD,NEQN,Y,XSTART,XFIN,RELERR,ABSERR,
     @        IFLAG,WORK,IWORK)
         IF (IFLAG.EQ.3.OR.IFLAG.EQ.4.OR.IFLAG.EQ.5)
     @              WRITE (11,1000) IFLAG,XFIN
         IFLAG=1
         DO 10 I=1,2
   10         YEIG(I,K)=Y(I)
   20 CONTINUE
      RT=X(NSURF)
      CALL SPLNTL(RT,VQ)
c The following is the outer boundary condition. We
c iterate until B2 is close enough to zero that
c corrections to the eigenvalue are within a small
c tolerance.
      B2=YEIG(1,NSURF)*(4.0D0+EIGT*R(NSURF)/G(NSURF))
     @   +YEIG(2,NSURF)
 1000 FORMAT (17H WATCH OUT,IFLAG=,I4,6H XFIN=,1PE10.2)
      RETURN
      END
c ***************************************************
c ***************************************************
      SUBROUTINE RKFCOW(RT,Y,YP)
c
c This routine supplies the derivatives for use in RKF
c for the nonradial case.
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      DOUBLE PRECISION L,LHAT,LINDEX
      DIMENSION Y(2),YP(2)
      COMMON/SPLINQ/VQ(6)
      COMMON /SURF/RSURF
      CALL SPLNTL(RT,VQ)
c
      YP(1)=Y(1)*(VQ(2)-3.D0+LINDEX)+Y(2)*(LHAT*VQ(1)
     @    /EIGT-VQ(2))
      YP(2)=Y(1)*(EIGT/VQ(1)-VQ(3))+Y(2)*(1.D0-VQ(4)
     @    +VQ(3)+LINDEX)
      DO 10 I=1,2
   10      YP(I)=YP(I)*VQ(5)
      RETURN
      END
c ****************************************************
c ****************************************************
      SUBROUTINE RKFRAD(RT,Y,YP)
c
c This routine supplies the derivatives for use in RKF
c for the radial case.
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      DOUBLE PRECISION L,LHAT,LINDEX
      DIMENSION Y(2),YP(2)
      COMMON/SPLINQ/VQ(6)
c
      CALL SPLNTL(RT,VQ)
      YP(1)=-(3.0D0*Y(1)+Y(2)/VQ(3))
      YP(2)=VQ(4)*(4.0D0*Y(1)+EIGT*Y(1)/VQ(1)+Y(2))
      DO 10 I=1,2
   10      YP(I)=YP(I)*VQ(5)
      RETURN
      END
c ****************************************************
c ****************************************************
      SUBROUTINE SPLNTL(XINT,YOUT)
c This subroutine interpolates model quantities for
c intermediate steps in the integration using cubic
c splines.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/SAVNGL/K,KQ
      COMMON /SAVEL/M,X(500),C1(4,500),Y1(500),C2(4,500),
     @   Y2(500),C3(4,500),Y3(500),C4(4,500),Y4(500),
     @   C5(4,500),Y5(500),C6(4,500)
      COMMON/RS/Y6(500)
      DIMENSION YOUT(6)
c
      MM=M-1
      IF (XINT.GE.X(1).AND.XINT.LT.X(M)) GO TO 20
      IF (XINT .LT. (1.+1.D-8)*X(1)) GO TO 60
      IF (XINT .GE. X(M)) GO TO 10
      K=1
      GO TO 50
   10 IF (XINT .GT. (1.+1.D-6)*X(M)) GO TO 60
      K=MM
      GO TO 50
   20 IL=1
      IR=M
   30 K=IL+((IR-IL)/2)
      IF (XINT.GE.X(K)) GO TO 40
      IR=K
      GO TO 30
   40 IF (XINT.LT.X(K+1)) GO TO 50
      IL=K
      GO TO 30
   50 CONTINUE
      X1=X(K+1)-XINT
      XX=XINT-X(K)
      X12=X1*X1
      XX2=XX*XX
      YOUT(1)=X1*(C1(1,K)*X12 +C1(3,K))+XX*(C1(2,K)*XX2
     @    +C1(4,K))
      YOUT(2)=X1*(C2(1,K)*X12 +C2(3,K))+XX*(C2(2,K)*XX2
     @    +C2(4,K))
      YOUT(3)=X1*(C3(1,K)*X12 +C3(3,K))+XX*(C3(2,K)*XX2
     @    +C3(4,K))
      YOUT(4)=X1*(C4(1,K)*X12 +C4(3,K))+XX*(C4(2,K)*XX2
     @    +C4(4,K))
      YOUT(5)=X1*(C5(1,K)*X12 +C5(3,K))+XX*(C5(2,K)*XX2
     @    +C5(4,K))
      YOUT(6)=X1*(C6(1,K)*X12 +C6(3,K))+XX*(C6(2,K)*XX2
     @    +C6(4,K))
      RETURN
   60 CONTINUE
      WRITE (16,1000) XINT
 1000 FORMAT (1H0,25HRANGE ERROR IN SPLINE, X=,E15.6)
      STOP
      END
c ***************************************************
c ***************************************************
      SUBROUTINE CONSL(X,Y,M,C)
c This subroutine sets up coefficient arrays for spline
c interpolation.
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION C(4,M)
      DIMENSION X(M),Y(M)
      DIMENSION A(500,3),D(500),B(500),Z(500),P(500)
c
      MM=M-1
      DO 10 K=1,MM
           D(K)=X(K+1)-X(K)
           P(K)=D(K)/6.D0
   10      Z(K)=(Y(K+1)-Y(K))/D(K)
      DO 20 K=2,MM
           B(K)=Z(K)-Z(K-1)
  20       CONTINUE
      A(1,2)=-1.D0-D(1)/D(2)
      A(1,3)=D(1)/D(2)
      A(2,3)=P(2)-P(1)*A(1,3)
      A(2,2)=2.D0*(P(1)+P(2))-P(1)*A(1,2)
      A(2,3)=A(2,3)/A(2,2)
      B(2)=B(2)/A(2,2)
      DO 30 K=3,MM
           A(K,2)=2.D0*(P(K-1)+P(K))-P(K-1)*A(K-1,3)
           B(K)=B(K)-P(K-1)*B(K-1)
           A(K,3)=P(K)/A(K,2)
           B(K)=B(K)/A(K,2)
   30 CONTINUE
      Q=D(M-2)/D(M-1)
      A(M,1)=1.D0+Q+A(M-2,3)
      A(M,2)=-Q-A(M,1)*A(M-1,3)
      B(M)=B(M-2)-A(M,1)*B(M-1)
      Z(M)=B(M)/A(M,2)
      MN=M-2
      DO 40 I=1,MN
           K=M-I
           Z(K)=B(K)-A(K,3)*Z(K+1)
   40 CONTINUE
      Z(1)=-A(1,2)*Z(2)-A(1,3)*Z(3)
      DO 50 K=1,MM
           Q=1.D0/(6.D0*D(K))
           C(1,K)=Z(K)*Q
           C(2,K)=Z(K+1)*Q
           C(3,K)=Y(K)/D(K)-Z(K)*P(K)
           C(4,K)=Y(K+1)/D(K)-Z(K+1)*P(K)
   50      CONTINUE
      RETURN
      END
c ***************************************************
c ***************************************************
      SUBROUTINE MODEID
c This subroutine identifies the mode of the pulsation
c using a phase diagram method. See Unno et al. (1989).
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON G(500),RHO(500),X(500),YEIG(2,500)
      DOUBLE PRECISION L,LHAT,LINDEX
      COMMON/MISC/L,LHAT,LINDEX,NSURF,PERIOD,GRAV,
     *    PI,PI4,P43,EPS,VERG,EIG,EIGT,NODES1,NODES2,
     *    MODEP
      INTEGER QUAD,QUAD2,QUAD1
c
      MODEP=0
c First count the nodes in Y1 and Y2.
      NODES1=0
      NODES2=0
      DO 10 I=2,NSURF
         FIT1=YEIG(1,I)*YEIG(1,I-1)
         FIT2=YEIG(2,I)*YEIG(2,I-1)
         IF (FIT1.LE.0) NODES1=NODES1+1
         IF (FIT2.LE.0) NODES2=NODES2+1
   10 CONTINUE
c Count crossings in the phase diagram.
      QUAD1=QUAD(YEIG(1,1),YEIG(2,1))
      DO 30 I=2,NSURF
         QUAD2=QUAD(YEIG(1,I),YEIG(2,I))
         IDQ=QUAD2-QUAD1
         IF (IDQ.EQ.0) GO TO 30
c See if quadrant chabge is a crossing in Y1.
         IF (IABS(IDQ).EQ.3) GO TO 20
         IF (QUAD1.EQ.3.AND.QUAD2.EQ.2) GO TO 20
         IF (QUAD1.EQ.2.AND.QUAD2.EQ.3) GO TO 20
c Not a Y1 crossing.
         QUAD1=QUAD2
         GO TO 30
c If crossing is clockwise, increment mode count.
c If crossing is counterclockwise, decrement mode count.
   20    IF (IDQ.EQ.-3.OR.IDQ.EQ.1) KDMOD=-1
         IF (IDQ.EQ.3.OR.IDQ.EQ.-1) KDMOD=1
         MODEP=MODEP+KDMOD
         QUAD1=QUAD2
c Get next point.
   30 CONTINUE
      RETURN
      END
c ****************************************************
c ****************************************************
      FUNCTION QUAD(Y1,Y2)
c Identifies the quadrant in the phase diagram
c used in MODEID.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER QUAD
c
      IF (Y1.GT.0.AND.Y2.GE.0) QUAD=1
      IF (Y1.LE.0.AND.Y2.GT.0) QUAD=2
      IF (Y1.LT.0.AND.Y2.LE.0) QUAD=3
      IF (Y1.GE.0.AND.Y2.LT.0) QUAD=4
      RETURN
      END
c ****************************************************
C ********************************************************
      SUBROUTINE RKF(F,NEQN,Y,T,TOUT,RELERR,ABSERR,
     @        IFLAG,WORK,IWORK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NEQN),WORK(1),IWORK(5)
      EXTERNAL F
      K1M=NEQN+1
      K1=K1M+1
      K2=K1+NEQN
      K3=K2+NEQN
      K4=K3+NEQN
      K5=K4+NEQN
      K6=K5+NEQN
      CALL RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,
     @    WORK(1),WORK(K1M),WORK(K1),WORK(K2),WORK(K3),
     @    WORK(K4),WORK(K5),WORK(K6),WORK(K6+1),
     @    IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5))
      RETURN
      END
C------------------------------------------------------
C
      SUBROUTINE RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,
     @    IFLAG,YP,H,F1,F2,F3,F4,F5,SAVRE,SAVAE,NFE,
     @    KOP,INIT,JFLAG,KFLAG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL HFAILD,OUTPUT
      DIMENSION Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),
     @    F3(NEQN),F4(NEQN),F5(NEQN)
      EXTERNAL F
      DATA U26/2.D-13/ , REMIN/2.D-13/
      DATA MAXNFE/3000/
      IF (NEQN .LT. 1) GO TO 10
      IF ((RELERR .LT. 0.D0)  .OR.  (ABSERR .LT. 0.D0))
     @   GO TO 10
      MFLAG=IABS(IFLAG)
      IF ((MFLAG .GE. 1) .AND. (MFLAG .LE. 7)) GO TO 20
   10 IFLAG=7
      RETURN
   20 IF (MFLAG .EQ. 1) GO TO 50
      IF (T .EQ. TOUT) GO TO 10
      IF(MFLAG .NE. 2) GO TO 25
      IF (INIT .EQ. 0) GO TO 45
      IF (KFLAG .EQ. 3) GO TO 40
      IF ((KFLAG .EQ. 4) .AND.  (ABSERR .EQ. 0.D0))
     @    GO TO 30
      IF ((KFLAG .EQ. 5)  .AND. (RELERR .LE. SAVRE)
     @    .AND. (ABSERR .LE. SAVAE)) GO TO 30
      GO TO 50
   25 IF (IFLAG .EQ. 3) GO TO 40
      IF ((IFLAG .EQ. 4) .AND. (ABSERR .GT. 0.D0))
     @    GO TO 45
   30 WRITE (6,1000) IFLAG,T
 1000 FORMAT (16H RKF SAYS IFLAG=,I5,3H X=,1PD12.4)
      STOP
   40 NFE=0
      IF (MFLAG .EQ. 2) GO TO 50
   45 IFLAG=JFLAG
   50 JFLAG=IFLAG
      KFLAG=0
      SAVRE=RELERR
      SAVAE=ABSERR
      RER=DMAX1(RELERR,REMIN)
      DT=TOUT-T
      IF (MFLAG .EQ. 1) GO TO 60
      IF (INIT .EQ. 0) GO TO 65
      GO TO 80
   60 INIT=0
      KOP=0
      A=T
      CALL F(A,Y,YP)
      NFE=1
      IF (T .NE. TOUT) GO TO 65
      IFLAG=2
      RETURN
   65 INIT=1
      YMAX=0.D0
      YPN=0.D0
      DO 70 K=1,NEQN
        YPN=DMAX1(DABS(YP(K)),YPN)
70      YMAX=DMAX1(DABS(Y(K)),YMAX)
      ETN=RER*YMAX+ABSERR
      H=DABS(DT)
      IF(ETN.GE.YPN*H**5) GO TO 80
      H=DMAX1((ETN/YPN)**0.2D0,U26*DMAX1(DABS(T),H))
80    H=DSIGN(H,DT)
      IF (DABS(H).GE.DABS(DT)) KOP=KOP+1
      IF (KOP.NE.100) GO TO 85
      IFLAG=6
      RETURN
85    IF (DABS(DT).GT.U26*DABS(T))GO TO 95
      DO 90 K=1,NEQN
90      Y(K)=Y(K)+DT*YP(K)
      A=TOUT
      CALL F(A,Y,YP)
      NFE=NFE+1
      GO TO 300
95    OUTPUT=.FALSE.
      SCALE=2.D0/RER
      AE=SCALE*ABSERR
100   HFAILD=.FALSE.
      HMIN=U26*DABS(T)
      DT=TOUT-T
      IF (DABS(DT).GE.2.D0*DABS(H)) GO TO 200
      IF (DABS(DT).GT.DABS(H)/0.9D0) GO TO 150
      OUTPUT=.TRUE.
      H=DT
      GO TO 200
150   H=0.5D0*DT
200   IF (NFE.LE.MAXNFE) GO TO 220
      IFLAG=3
      KFLAG=3
      RETURN
220   CALL FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,F1)
      NFE=NFE+5
      EEOET=0.D0
      DO 250 K=1,NEQN
        ET=DABS(Y(K))+DABS(F1(K))+AE
        IF (ET.GT.0.D0) GO TO 240
        IFLAG=4
        KFLAG=4
        RETURN
240     EE=DABS((-2090.D0*YP(K)+(21970.D0*F3(K)-15048.D0
     @     *F4(K)))+(22528.D0*F2(K)-27360.D0*F5(K)))
250     EEOET=DMAX1(EEOET,EE/ET)
      ESTTOL=DABS(H)*EEOET*SCALE/752400.D0
      IF (ESTTOL.LE.1.D0) GO TO 260
      HFAILD=.TRUE.
      OUTPUT=.FALSE.
      S=0.1D0
      IF (ESTTOL.LT.59049.D0)S=0.9D0/ESTTOL**0.2D0
      H=S*H
      IF (DABS(H).GT.HMIN) GO TO 200
      IFLAG=5
      KFLAG=5
      RETURN
260   T=T+H
      DO 270 K=1,NEQN
270     Y(K)=F1(K)
      A=T
      CALL F(A,Y,YP)
      NFE=NFE+1
      IF (HFAILD) GO TO 290
      S=5.D0
      IF (ESTTOL.GT.1.889568D-04) S=0.9D0/ESTTOL**0.2D0
      H=DSIGN(DMAX1(S*DABS(H),HMIN),H)
290   IF (OUTPUT) GO TO 300
      IF (IFLAG.GT.0) GO TO 100
      IFLAG=-2
      RETURN
300   T=TOUT
      IFLAG=2
      RETURN
      END
C-----------------------------------------------------
C
      SUBROUTINE FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,S)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),
     @   F3(NEQN),F4(NEQN),F5(NEQN),S(NEQN)
      CH=0.25D0*H
      DO 10 K=1,NEQN
   10      F5(K)=Y(K)+CH*YP(K)
      CALL F(T+0.25D0*H,F5,F1)
      CH=0.09375D0*H
      DO 20 K=1,NEQN
   20      F5(K)=Y(K)+CH*(YP(K)+3.D0*F1(K))
      CALL F(T+0.375D0*H,F5,F2)
      CH=H/2197.D0
      DO 30 K=1,NEQN
   30      F5(K)=Y(K)+CH*(1932.D0*YP(K)+(7296.D0*F2(K)
     @     -7200.D0*F1(K)))
      CALL F(T+12.D0/13.D0*H,F5,F3)
      CH=H/4104.D0
      DO 40 K=1,NEQN
   40      F5(K)=Y(K)+CH*((8341.D0*YP(K)-845.D0*F3(K))+
     @        (29440.D0*F2(K)-32832.D0*F1(K)))
      CALL F(T+H,F5,F4)
      CH=H/20520.D0
      DO 50 K=1,NEQN
   50    F1(K)=Y(K)+CH*((-6080.D0*YP(K)+(9295.D0*F3(K)
     @     -5643.D0*F4(K)))+(41040.D0*F1(K)-28352.D0
     @     *F2(K)))
      CALL F(T+0.5D0*H,F1,F5)
      CH=H/7618050.D0
      DO 60 K=1,NEQN
   60     S(K)=Y(K)+CH*((902880.D0*YP(K)+(3855735.D0
     @      *F3(K)-1371249.D0*F4(K)))+(3953664.D0*F2(K)
     @      +277020.D0*F5(K)))
      RETURN
      END
C *************************************************