      SUBROUTINE ALCON2 (FY,JFY,N,X,XW,TAU,TAUMIN,TAUMAX,UMAX,
     &   EPS,INFO,RWORK,LRWORK,IWORK,LIWORK)
C*    Begin Prologue ALCON2
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL FY,JFY
      DIMENSION RWORK(LRWORK),IWORK(LIWORK),X(N),XW(N),INFO(9)
C
C---------------------------------------------------------------------
C
C*  Title
C
C     (Al)gebraic system of equations (Con)tinuation method.
C
C*  Written by        P. Deuflhard, P. Kunkel
C*  Purpose           Solution of parameter dependent systems of 
C                     nonlinear equations.
C*  Method            Numerical pathfollowing with automatic steplength
C                     control
C*  Category          F4 - Parameter Dependent Nonlinear Equation 
C                          Systems
C*  Keywords          Numerical pathfollowing, Homotopy Method
C*  Version           0.9
C*  Revision          September 1985
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C*  Copyright     (c) Konrad-Zuse-Zentrum fuer
C                     Informationstechnik Berlin (ZIB)
C                     Takustrasse 7, D-14195 Berlin-Dahlem
C                     phone : + 49/30/84185-0
C                     fax   : + 49/30/84185-125
C*  Contact           Bodo Erdmann
C                     ZIB, Division Scientific Computing, 
C                          Department Numerical Analysis and Modelling
C                     phone : + 49/30/84185-185
C                     fax   : + 49/30/84185-107
C                     e-mail: erdmann@zib.de
C
C*    Reference:
C
C     /1/ P. Deuflhard, B. Fiedler, P. Kunkel:
C         Efficient Numerical Pathfollowing Beyond Critical Points
C         Univ. Heidelberg, Sfb 123, Tech. Rep. 278 (1984)
C
C  ---------------------------------------------------------------
C
C* Licence
C    You may use or modify this code for your own non commercial
C    purposes for an unlimited time.
C    In any case you should not deliver this code without a special
C    permission of ZIB.
C    In case you intend to use the code commercially, we oblige you
C    to sign an according licence agreement with ZIB.
C
C* Warranty
C    This code has been tested up to a certain level. Defects and
C    weaknesses, which may be included in the code, do not establish
C    any warranties by ZIB. ZIB does not take over any liabilities
C    which may follow from aquisition or application of this code.
C
C* Software status
C    This code is under partial care of ZIB and belongs to ZIB
C    software class 2.
C
C     ------------------------------------------------------------
C
C*   Summary
C
C     Continuation method for systems of algebraic equations  f(x,tau)=0
C     Optional computation of turning and (simple) bifurcation points.
C     Optional automatic construction of complete bifurcation diagrams.
C     ------------------------------------------------------------
C
C  INPUT PARAMETERS  (* MARKS INOUT PARAMETERS)
C
C    - FY (X,TAU,F)
C               EXTERNAL SUBROUTINE DEFINING SYSTEM OF
C               ALGEBRAIC EQUATIONS
C                 X(N)       VECTOR OF STATE VARIABLES    INPUT
C                 TAU        PARAMETER                    INPUT
C                 F(N)       VALUES OF F(X,TAU)           OUTPUT
C
C    - JFY (X,TAU,DFDX,DFDTAU)
C               EXTERNAL SUBROUTINE DEFINING JACOBIAN OF FUNCTION F
C                 X(N)       VECTOR OF STATE VARIABLES    INPUT
C                 TAU        PARAMETER                    INPUT
C                 DFDX(N,N)  DERIVATIVES OF F WRT X       OUTPUT
C                 DFDTAU(N)  DERIVATIVES OF F WRT TAU     OUTPUT
C
C    - N        NUMBER OF ALGEBRAIC EQUATIONS
C
C    - X(N)   * ESTIMATE OF SOLUTION ON CONTINUATION PATH
C               FOR INITIAL TAU
C
C    - XW(N)    INITIAL SCALING QUANTITIES FOR X(N)
C               (ALCON2 USES ADAPTIVE INTERNAL SCALING
C               ALONG THE CONTINUATION PATH)
C
C    - TAU    * INITIAL VALUE OF PARAMETER TAU
C
C    - TAUMIN   MINIMUM ALLOWED VALUE FOR TAU
C
C    - TAUMAX   MAXIMUM ALLOWED VALUE FOR TAU
C
C    - UMAX     MAXIMUM ALLOWED L2-NORM FOR SOLUTION VECTOR X
C               ALONG THE CONTINUATION PATH
C               (SET UMAX=0.D0  IF NO RESTRICTION DESIRED)
C
C    - EPS      REQUIRED RELATIVE ACCURACY
C
C    - INFO(9)  ARRAY USED FOR COMMUNICATION BETWEEN
C               THE PARTICULAR ROUTINES
C                 INFO(1)  PRINT PARAMETER
C                            0  NO PRINT
C                            1  PRINT OF X(N), TAU, AND INFORMATION
C                               ILLUSTRATING THE CONTINUATION PROCESS
C                            2  ADDITIONALLY INFORMATION ABOUT THE
C                               CONVERGENCE BEHAVIOR OF THE
C                               GAUSS-NEWTON METHOD
C                 INFO(2)  MAXIMUM PERMITTED NUMBER OF
C                          CONTINUATION STEPS
C                 INFO(3)  MODE PARAMETER
C                            0  AUTOMATIC CONSTRUCTION OF COMPLETE
C                               BIFURCATION DIAGRAM  (STANDARD MODE)
C                           +1  PATHFOLLOWING ALONG A SINGLE BRANCH
C                               STARTING IN POSITIVE TAU-DIRECTION
C                               (ALCON1 MODE)
C                           -1  PATHFOLLOWING ALONG A SINGLE BRANCH
C                               STARTING IN NEGATIVE TAU-DIRECTION
C                               (ALCON1 MODE)
C                 INFO(4)  CRITICAL POINT OPTIONS
C                            0  NO CRITICAL POINT DETERMINATION
C                            1  TURNING POINT COMPUTATION
C                            2  BIFURCATION POINT COMPUTATION
C                            3  TURNING AND BIFURCATION POINT
C                               COMPUTATION  (RECOMMENDED OPTION)
C                          IF INFO(3)=0, THE COMPUTATION OF BIFURCATION
C                          POINTS IS ANYWAY PERFORMED
C                 INFO(5)  ERROR EXIT PARAMETER  (SEE OUTPUT PARAMETERS)
C                 INFO(6)  INTERNALLY USED
C                 INFO(7)  INTERNALLY USED
C                 INFO(8)  INTERNALLY USED
C                 INFO(9)  INTERNALLY USED
C
C    - RWORK    REAL WORKSPACE
C
C    - LRWORK   LENGTH OF REAL WORKSPACE
C               MUST BE AT LEAST   7 * N**2 + 35 * N + 24
C
C    - IWORK    INTEGER WORKSPACE
C
C    - LIWORK   LENGTH OF INTEGER WORKSPACE
C               MUST BE AT LEAST   3 * N + 2 + NBIF ,
C               WHERE NBIF=0 IN THE CASE OF PURE PATHFOLLOWING
C               OR NBIF=NUMBER OF DETECTED BIFURCATION POINTS
C               IN THE CASE INFO(3)=0
C
C    REMARK: FOR EASE OF IMPLEMENTATION ONLY REGULAR SOLUTIONS ARE
C            ASSUMED NEAR INITIAL TAU, TAUMIN, AND TAUMAX. COMPUTATIONS
C            ARE PERFORMED STRICTLY BETWEEN TAUMIN AND TAUMAX.
C
C
C  OUTPUT PARAMETERS
C
C    - X(N)     FINAL SOLUTION VALUES
C
C    - TAU      FINAL PARAMETER VALUE
C
C    - INFO(5)  ERROR EXIT PARAMETER
C                 0  NO ERROR OCCURRED
C                 1  MORE THAN INFO(2) CONTINUATION STEPS
C                 2* STEPLENGTH IN CONTINUATION PROCESS TOO SMALL
C                    (RELATIVE DIFFERENCE LESS THAN 10*EPS)
C                 3  NOT USED IN ALCON2
C                 4  INITIAL GUESS OF X(N) TOO BAD
C                 5* START-OFF PRECEDURE FAILED
C                 6  WORKSPACE TOO SMALL
C                       * WARNING ONLY FOR INFO(3)=0 -
C                         PROCESSING CONTINUED
C
C  EXTERNAL UNITS
C
C      COMMON /A2UNIT/ UPR,UDIAG,UBIF
C
C      - UPR     PRINT UNIT
C                (STANDARD UNIT 6)
C      - UDIAG   PLOT INFORMATION UNIT, INPUT TO PLOT ROUTINE PLTHM
C                DISC FILE WITH CARD IMAGE  TO BE DECLARED BY THE USER
C                SET UDIAG=0, IF NO PLOT IS DESIRED
C                (STANDARD UNIT 2)
C      - UBIF    BRANCHING INFORMATION UNIT, NECESSARY
C                ONLY IN THE CASE  INFO(3)=0
C                DISC FILE WITH CARD IMAGE  TO BE DECLARED BY THE USER
C                (STANDARD UNIT 1)
C
C     ------------------------------------------------------------
C*    End Prologue
Crr      INTEGER UPR,UDIAG,UBIF
      INTEGER UPR,UDIAG
      COMMON /MACHIN/ EPMACH,SMALL
Crr      COMMON /A2UNIT/ UPR,UDIAG,UBIF
      COMMON /A2UNIT/ UPR,UDIAG
C  MACHINE DEPENDENT CONSTANTS  (ADAPTED TO IBM3081D)
C            EPMACH  UNIT ROUND OFF
C            SMALL   SQRT OF SMALLEST MACHINE NUMBER
      CALL ZIBCONST(EPMACH,SMALL)
C  WORKSPACE REQUIREMENTS
      N1=N+1
      NN1=N*N1
      NH=2*N1
      LRWREQ=7*N**2+35*N+24
      N2=2*N
      LIWREQ=3*N+2
      IF (LRWREQ.GT.LRWORK .OR. LIWREQ.GT.LIWORK) GOTO 9000
C  REAL WORKSPACE DISTRIBUTION
      LY=1
      LYQ=LY+N1
      LYD=LYQ+N1
      LYDH=LYD+N1
      LYH=LYDH+N1
      LDYT=LYH+N1
      LDYTA=LDYT+N1
      LDYTH=LDYTA+N1
      LV=LDYTH+N1
      LD=LV+N1
      LYW=LD+N1
      LETA=LYW+N1
      LDFY=LETA+N1
      LA=LDFY+NN1
      LAH=LA+N*(N2+1)
      LYA=LAH+N1*N1
      LDY=LYA+N1
      LDYQ=LDY+NH
      LW=LDYQ+NH
      LF=LW+N2
      LFH=LF+NH
      LU=LFH+N
      LZ=LU+NH
      LZA=LZ+N
      LZH=LZA+N
      LDFYH=LZH+N
      LB=LDFYH+NN1
      LW1=LB+N*(N2+1)
      LW2=LW1+N1
C  INTEGER WORKSPACE DISTRIBUTION
      LIPIV=1
      LLT=LIPIV+N1
      LINDEX=LLT+2*N
      NINDEX=LIWORK-LINDEX+1
C  REORGANIZE INITIAL VALUE
      DO 1000 I=1,N
      RWORK(I)=X(I)
1000  CONTINUE
      RWORK(N1)=TAU
C  EXTERNAL INITIAL SCALING
      DO 2000 I=1,N
      RWORK(LYW+I-1)=XW(I)
2000  CONTINUE
      RWORK(LYW+N)=DMAX1(DABS(TAUMIN),DABS(TAUMAX))
C  CALL CONTINUATION ROUTINE
      CALL A2HOMQ(FY,JFY,N,N1,NH,TAUMIN,TAUMAX,UMAX,EPS,INFO,
     &   RWORK(LY),RWORK(LYQ),RWORK(LYD),RWORK(LYDH),RWORK(LYH),
     &   RWORK(LDYT),RWORK(LDYTA),RWORK(LDYTH),
     &   RWORK(LV),RWORK(LD),RWORK(LYW),RWORK(LETA),
     &   RWORK(LDFY),RWORK(LA),RWORK(LAH),
     &   RWORK(LYA),RWORK(LDY),RWORK(LDYQ),RWORK(LW),
     &   RWORK(LF),RWORK(LFH),RWORK(LU),
     &   RWORK(LZ),RWORK(LZA),RWORK(LZH),
     &   RWORK(LDFYH),RWORK(LB),
     &   RWORK(LW1),RWORK(LW2),
     &   N2,NINDEX,IWORK(LIPIV),IWORK(LLT),IWORK(LINDEX))
C  RESTORE FINAL VALUES
      DO 3000 I=1,N
      X(I)=RWORK(I)
3000  CONTINUE
      TAU=RWORK(N1)
      RETURN
C  FAIL EXIT
9000  CONTINUE
      INFO(5)=6
      KPRINT=INFO(1)
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0 .AND. LRWREQ.GT.LRWORK) WRITE(UPR,60002) LRWREQ
      IF (KPRINT.GT.0 .AND. LIWREQ.GT.LIWORK) WRITE(UPR,60003) LIWREQ
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      RETURN
60001 FORMAT(///)
60002 FORMAT('   REAL  WORKSPACE NOT SUFFICIENT',4X,'AT LEAST',I5,
     &   ' LOCATIONS NEEDED')
60003 FORMAT(' INTEGER WORKSPACE NOT SUFFICIENT',4X,'AT LEAST',I5,
     &   ' LOCATIONS NEEDED')
      END
C
C
      SUBROUTINE A2HOMQ (FY,JFY,N,N1,NH,TAUMIN,TAUMAX,UMAX,EPS,INFO,
     &   Y,YQ,YD,YDH,YH,DYT,DYTA,DYTH,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,
     &   F,FH,U,Z,ZA,ZH,DFYH,B,W1,W2,N2,INDMAX,IPIV,LT,INDEX)
C
C  SUBROUTINE  A2HOMQ  TO BE USED WITH ROUTINE ALCON2
C
C  DRIVER ROUTINE FOR ALCON2
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),YQ(N1),YD(N1),YDH(N1),YH(N1)
      DOUBLE PRECISION DYT(N1),DYTA(N1),DYTH(N1)
      DOUBLE PRECISION V(N1),D(N1),YW(N1),ETA(N1)
      DOUBLE PRECISION DFY(N,N1),A(1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(NH),DYQ(NH),W(1)
      DOUBLE PRECISION F(NH),FH(N),U(NH)
      DOUBLE PRECISION Z(N),ZA(N),ZH(N)
      DOUBLE PRECISION DFYH(N,N1),B(1)
      DOUBLE PRECISION W1(N1),W2(N1)
      INTEGER INFO(9),IPIV(N1),LT(N2),INDEX(INDMAX)
Crr      INTEGER UPR,UDIAG,UBIF
      INTEGER UPR,UDIAG
      EXTERNAL FY,JFY
      COMMON /COUNT/ IFCTEV,IJACEV,ITER,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
Crr      COMMON /A2UNIT/ UPR,UDIAG,UBIF
      COMMON /A2UNIT/ UPR,UDIAG
      COMMON /A2VALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      COMMON /A2METH/ FACTOR,IEXTR,IJACM
      DATA SUMYH/1.D0/,SUMYQH/1.D0/,VAL1/1.D0/,VAL2/1.D0/
C-----------------------------------------------------------------------
C  PREPARATIONS
      EPDIFF=1.D3*EPMACH
      ETADIF=DSQRT(EPDIFF)
      ITER=0
      IHALT=0
      YHALT = 0.0D0
      KPRINT=INFO(1)
      DO 1010 I=1,N1
      ETA(I)=ETADIF
1010  CONTINUE
      TAUMAH=TAUMAX
      IFCTEV=0
      IJACEV=0
      ISOLS=0
      IDECS=0
      IWARN=0
      INFO(6)=IGNMAX
C  CHECK FOR CONSISTENCY OF INITIAL VALUES
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60201)
      IF (KPRINT.GT.2) WRITE(UPR,60500) Y(N1),(Y(I),I=1,N)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INFO(8)=1
      CALL A2GNHM(N,N1,Y,FY,JFY,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,COND)
      IERR=INFO(9)
      IF (IERR.NE.0) GOTO 9940
      ITS=INFO(7)
      CALL A2TAND(N,N1,JFY,Y,YW,DFY,A,COND,DETFA,IPIV,AH,V,
     &   D,IRANK)
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) GOTO 9940
      DETHMA=DETFA
      IF (IPIVS.NE.N1) CALL A2DETH(N,N1,A,AH,D,N1,IPIV,DETHMA)
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
      IF (KPRINT.GT.0) WRITE(UPR,60500) Y(N1),(Y(I),I=1,N)
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
C.......................................................................
C  START OF BRANCH FOLLOWING
1100  CONTINUE
      INIT=0
      ITER=ITER+1
      IF (ITER.GT.INFO(2)) GOTO 9910
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60202) ITER
      DO 1110 I=1,N1
      YQ(I)=Y(I)
1110  CONTINUE
C.......................................................................
C  INITIAL ESTIMATE FOR STEPLENGTH
      IF (ITER.NE.1) GOTO 2200
      DO 1210 I=1,N
      DYT(IPIV(I))=-V(I)/D(N1)
1210  CONTINUE
      DYT(IPIVS)=1.D0/D(N1)
      SIGMA=SMALL1
      SIGNUM=1.D0
C-    IF (INFO(3).NE.0) SIGNUM=DBLE(INFO(4))   KORR.
      IF (INFO(3).NE.0) SIGNUM=DBLE(INFO(3))
      IF (Y(N1).EQ.TAUMAX) SIGNUM=-1.D0
      IF (DYT(N1)*SIGNUM.GT.0.D0) GOTO 1221
      DO 1220 I=1,N1
      DYT(I)=-DYT(I)
1220  CONTINUE
1221  CONTINUE
      IPIVA=IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20102) IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) Y,DYT,YW
Crr      IF (INFO(3).EQ.0) WRITE(UBIF,10001) Y,YW,ETA
      IF (INFO(3).EQ.0) CALL A2RWR(N1,Y,YW,ETA)
      DO 1230 I=1,N1
      W(I)=-DYT(I)
1230  CONTINUE
Crr      IF (INFO(3).EQ.0) WRITE(UBIF,10001) Y,(W(I),I=1,N1),(W(I),I=1,N1)
Crr      IF (INFO(3).EQ.0) WRITE(UBIF,10002) IPIVA,SIGMA,DETHMA,DETFA
      IF (INFO(3).EQ.0) Call A2RWR(N1,Y,W,W)
      IF (INFO(3).EQ.0) CALL A2RWR4(IPIVA,SIGMA,DETHMA,DETFA)
      IND=1
      INDEX(1)=1
      TAUH=Y(N1)
      IF (TAUH.EQ.TAUMIN .OR. TAUH.EQ.TAUMAX) INDEX(1)=0
      IF (INFO(3).NE.0) INDEX(1)=0
      GOTO 2200
C-----------------------------------------------------------------------
C  STEPLENGTH PREDICTOR
2000  CONTINUE
      ITER=ITER+1
      IF (ITER.GT.INFO(2)) GOTO 9910
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60202) ITER
      IF (ITS.EQ.0) GOTO 2020
      IF (VAL2.LT.THR*VAL1) VAL2=THR*VAL1
      R=DSQRT(VAL1/VAL2)
      SIGMA=R*SIGMA
      GOTO 2030
C  EMPIRICAL STEPLENGTH INCREASE IN NEARLY LINEAR CASE
2020  CONTINUE
      SIGMA=SIGMA/DSQRT(THR)
C  EXTRAPOLATED STEPLENGTH BOUND
2030  CONTINUE
      IF (IHALT.EQ.0) GOTO 2040
      SIGMAQ=FACTOR*(YHALT-YQ(IPIVA))/(DYT(IPIVA)*YW(IPIVA))
      IF (SIGMAQ.LE.10.D0*EPS .OR. SIGMA.LE.SIGMAQ) GOTO 2040
      SIGMA=SIGMAQ
      IF (KPRINT.GT.1) WRITE(UPR,60206)
2040  CONTINUE
      IF (INIT.NE.0) IPIVA=IPIVS
      GOTO 2200
C.......................................................................
C  STEPLENGTH CORRECTOR
2110  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60401)
      TH=DSQRT(SUMYQH/SUMYH)
      R=DSQRT(0.5D0*THMAX/TH)
      IF (R.GT.0.7D0) R=0.7D0
      IF (R.LT.0.1D0) R=0.1D0
      GOTO 2190
2120  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60402)
      R=0.7D0
      GOTO 2190
2130  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60403)
      R=0.7D0
2190  CONTINUE
      IF (SIGMA.LT.EPS) GOTO 9920
      INIT=1
      SIGMA=R*SIGMA
C.......................................................................
C  TRIAL VALUE
2200  CONTINUE
      IF (INDEX(1).EQ.1 .AND. YQ(N1).LT.TAUH) TAUMAX=TAUH
2210  CONTINUE
      SIGMIN=0.D0
      SIGMAX=0.D0
      IF (DYT(N1).EQ.0.D0) GOTO 2211
      SIGMAX=(TAUMAX-YQ(N1))/(DYT(N1)*YW(N1))
      IF (SIGMAX.GT.0.D0 .AND. SIGMA.GT.SIGMAX) SIGMA=SIGMAX
      SIGMIN=(TAUMIN-YQ(N1))/(DYT(N1)*YW(N1))
      IF (SIGMIN.GT.0.D0 .AND. SIGMA.GT.SIGMIN) SIGMA=SIGMIN
2211  CONTINUE
      DO 2220 I=1,N1
      YD(I)=YQ(I)+SIGMA*DYT(I)*YW(I)
      YDH(I)=YD(I)
2220  CONTINUE
C.......................................................................
C  CHECK FOR POSSIBLE CYCLE
      IF (TAUMAX.NE.TAUH .OR. INDEX(1).EQ.0) GOTO 2400
Crr      REWIND UBIF
      CALL A2RREW
cno
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      S=0.D0
      DO 2310 I=1,N1
      S=S+((YH(I)-YD(I))/YW(I))**2
2310  CONTINUE
      IF (IND.EQ.1) GOTO 2321
      DO 2320 L=2,IND
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
2320  CONTINUE
2321  CONTINUE
      IF (S.LE.4.D0*SUMYH) GOTO 2400
      TAUMAX=TAUMAH
      GOTO 2210
C.......................................................................
C  RESCALING
2400  CONTINUE
      DO 2410 I=1,N1
      T=0.5D0*(DABS(Y(I))+DABS(YQ(I)))
      IF (T.LT.YW(I)) T=YW(I)
      DYTH(I)=DYT(I)*YW(I)/T
      YW(I)=T
2410  CONTINUE
      DO 2420 I=1,N1
      Y(I)=YQ(I)
2420  CONTINUE
      IF (SIGMA.EQ.SIGMIN .OR. SIGMA.EQ.SIGMAX) GOTO 3100
C-----------------------------------------------------------------------
C  ITERATION BACK TO CONTINUATION PATH BY GAUSS-NEWTON METHOD
      IF (KPRINT.GT.0) WRITE(UPR,60205) SIGMA,IPIVA
      IF (KPRINT.GT.2) WRITE(UPR,60500) YD(N1),(YD(I),I=1,N)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INFO(8)=2
      CALL A2GNHM(N,N1,YD,FY,JFY,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,COND)
      IERR=INFO(9)
      ITS=INFO(7)
      IF (IERR.NE.0 .AND. INIT.EQ.0) GOTO 2110
      IF (ITS.NE.0) INIT=1
      IF (INIT.EQ.0) GOTO 2020
      GOTO (2110,2110,5110,3110),IERR
C  PREPARE STEPLENGHT PREDICTION FOR NEXT CONTINUATION STEP
      IF (ITS.EQ.0) GOTO 3030
      VAL1=0.D0
      DO 3010 I=1,N1
      VAL1=VAL1+DYT(I)*V(I)
3010  CONTINUE
      TH=DSQRT(SUMYQH/SUMYH)
      VAL1=0.5D0*THMAX*DSQRT(SUMYH)/DABS(VAL1)
      VAL2=0.D0
      DO 3020 I=1,N1
      VAL2=VAL2+((YDH(I)-YD(I))/YW(I))**2
3020  CONTINUE
      VAL2=TH*DSQRT(VAL2)
3030  CONTINUE
C  COMPUTE TANGENT AND DETERMINANT
      CALL A2TAND(N,N1,JFY,YD,YW,DFY,A,COND,DETF,IPIV,AH,V,
     &   D,IRANK)
      IF (IRANK.LT.N .AND. INIT.EQ.0) GOTO 2110
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) GOTO 5110
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
      TAUMAX=TAUMAH
      GOTO 4000
C.......................................................................
C  HIT FINAL VALUE
3100  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60205) SIGMA,IPIVA
      IF (SIGMA.EQ.SIGMIN) YD(N1)=TAUMIN
      IF (SIGMA.EQ.SIGMAX) YD(N1)=TAUMAX
      GOTO 3120
3110  CONTINUE
      IF (YD(N1).LT.TAUMIN) YD(N1)=TAUMIN
      IF (YD(N1).LT.TAUMAX) YD(N1)=TAUMAX
3120  CONTINUE
      IF (KPRINT.GT.2) WRITE(UPR,60500) YD(N1),(YD(I),I=1,N)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INIT=2
      INFO(8)=3
      CALL A2GNHM(N,N1,YD,FY,JFY,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,COND)
      IERR=INFO(9)
      ITS=INFO(7)
      IF (IERR.NE.0) GOTO 2110
C  COMPUTE TANGENT AND DETERMINANT
      CALL A2TAND(N,N1,JFY,YD,YW,DFY,A,COND,DETF,IPIV,AH,V,
     &   D,IRANK)
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) GOTO 2110
      TAUMAX=TAUMAH
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
C-----------------------------------------------------------------------
C  DETERMINANTS
4000  CONTINUE
      DETHM=DETF
      IF (IPIVS.NE.N1) CALL A2DETH(N,N1,A,AH,D,N1,IPIV,DETHM)
      DETFH=DETF
      IF (IPIVS.NE.IPIVA) CALL A2DETH(N,N1,A,AH,D,IPIVA,IPIV,DETFH)
      IF (KPRINT.GT.1) WRITE(UPR,60501) DETFA,DETFH,DETHMA,DETHM
C.......................................................................
C  NEW NORMALIZED TANGENT
      DO 4110 I=1,N
      W(IPIV(I))=-V(I)/D(N1)
4110  CONTINUE
      W(IPIVS)=1.D0/D(N1)
C  CHECK SIGN OF TANGENT
      SIGNUM=DSIGN(1.0D0,YD(IPIVS)-Y(IPIVS))
      IF (W(IPIVS)*SIGNUM.GT.0.D0) GOTO 4121
      DO 4120 I=1,N1
      W(I)=-W(I)
4120  CONTINUE
4121  CONTINUE
C.......................................................................
C  CHECK FOR CONFLICT IN SIGNS
      IF (SIGMA.LE.10.D0*EPS) GOTO 4210
      IF (DYT(IPIVS)*W(IPIVS).LT.0.D0) GOTO 2120
      IF (DYT(IPIVA)*W(IPIVA).LT.0.D0) GOTO 2120
4210  CONTINUE
C.......................................................................
C  SAVE VALUES FOR NEXT ITERATE
      DO 4310 I=1,N1
      DYTA(I)=DYT(I)
      DYT(I)=W(I)
      YQ(I)=YD(I)
4310  CONTINUE
C-----------------------------------------------------------------------
C  CHECK FOR TURNING POINT
      IF (DYTA(N1)*DYT(N1).GT.0.D0) GOTO 5010
      CALL A2TURN(N,N1,IPIVA,EPS,INFO,
     &   Y,W1,YD,YH,DYT,DYTH,W2,V,YW,
     &   IPIV,FY,JFY,TAUMIN,TAUMAX,D,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U)
      IERRQ=INFO(9)
      IF (IERRQ.GT.0 .AND. SIGMA.GE.10.D0*EPS) GOTO 5200
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20103) IPIV(N1)
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0)
     &   WRITE(UDIAG,20000) YH,(W(I),I=1,N1),YW
5010  CONTINUE
C.......................................................................
C  CHECK FOR BIFURCATION POINT
      IF (DETFH*DETFA.GT.0.D0) GOTO 5990
      INFO(8)=1
      GOTO 5111
5110  CONTINUE
      INFO(8)=2
5111  CONTINUE
      DO 5112 I=1,N1
      YD(I)=YQ(I)
5112  CONTINUE
      CALL A2BIFU(N,N1,NH,N2,IPIVA,EPS,INFO,Y,YD,YH,
     &   V,YW,IPIV,LT,FY,JFY,TAUMIN,TAUMAX,D,ETA,DFY,A,AH,YA,DY,DYQ,W,
     &   F,FH,U,Z,ZA,ZH,DFYH,B,DETFH,DETFA,W1,W2)
      IERRQ=INFO(9)
      IF (IERRQ.EQ.0) GOTO 5300
      IF (IERRQ.NE.4) GOTO 5200
C  POSSIBLY NON-SIMPLE BIFURCATION POINT DETECTED
      IF (KPRINT.GT.0) WRITE(UPR,60609)
      GOTO 5990
C.......................................................................
C  RESTORE FORMER VALUES
5200  CONTINUE
      DO 5210 I=1,N1
      YQ(I)=Y(I)
      DYT(I)=DYTA(I)
5210  CONTINUE
      INIT=1
      GOTO 2130
C.......................................................................
C  CHECK IF DETECTED BIFURCATION POINT HAS ALREADY OCCURRED
5300  CONTINUE
      IF (INFO(3).NE.0) GOTO 5990
Crr      REWIND UBIF
      CALL A2RREW
Crr      READ(UBIF,10001) YD,V,(W(I),I=1,N1)
Crr      READ(UBIF,10001) YD,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YD,V,W)
      CALL A2RRD(N1,YD,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      IF (IND.EQ.1) GOTO 5400
      IBIF=0
      DO 5390 L=2,IND
Crr      READ(UBIF,10001) YD,V,(W(I),I=1,N1)
      CALL A2RRD(N1,YD,V,W)
      IF (IBIF.NE.0) GOTO 5370
      VN=0.D0
      DO 5310 I=1,N1
      VN=VN+((YD(I)-YH(I))/YW(I))**2
5310  CONTINUE
      VN=DSQRT(VN)
      IF (VN.GT.10.D0*EPS) GOTO 5370
C  OLD BIFURCATION POINT
      L1=L-1
      IF (KPRINT.GT.0) WRITE(UPR,60601) L1
      IF (INDEX(L).EQ.0) IWARN=IWARN+1
      IF (KPRINT.GT.0 .AND. INDEX(L).EQ.0) WRITE(UPR,60608)
      IBIF=1
      IF (INDEX(L).LE.1) IBIF=2
      INDEX(L)=0
      IF (IBIF.EQ.2) GOTO 5370
      DO 5320 I=1,N1
      YW(I)=V(I)
      ETA(I)=W(I)
5320  CONTINUE
Crr      READ(UBIF,10001) YD,DYTA,DYT
Crr      READ(UBIF,10002) IPIVS,SIGMA,DETF,DETHM
Crr      READ(UBIF,10002) ITS,VAL1,VAL2
      CALL A2RRD(N1,YD,DYTA,DYT)
      CALL A2RRD4(IPIVS,SIGMA,DETF,DETHM)
      CALL A2RRD3(ITS,VAL1,VAL2)
      S=0.D0
      DO 5330 I=1,N1
      S=S+(Y(I)-YH(I))/YW(I)*DYTA(I)
5330  CONTINUE
      IF (S.LE.0.D0) GOTO 5340
Crr      READ(UBIF,10001) YD,DYTA,DYT
Crr      READ(UBIF,10002) IPIVS,SIGMA,DETF,DETHM
Crr      READ(UBIF,10002) ITS,VAL1,VAL2
      CALL A2RRD(N1,YD,DYTA,DYT)
      CALL A2RRD4(IPIVS,SIGMA,DETF,DETHM)
      CALL A2RRD3(ITS,VAL1,VAL2)
      GOTO 5350
5340  CONTINUE
Crr      READ(UBIF,10001) YQ,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YQ,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
5350  CONTINUE
      DO 5360 I=1,N1
      Y(I)=YH(I)
      YQ(I)=YD(I)
5360  CONTINUE
      DETFA=0.D0
      GOTO 5390
C  READ IN DUMMY
5370  CONTINUE
Crr      READ(UBIF,10001) YD,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10001) YD,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YD,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
      CALL A2RRD(N1,YD,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
5390  CONTINUE
      IF (UDIAG.GT.0 .AND. IBIF.EQ.1) WRITE(UDIAG,20101) IPIVS
      IF (UDIAG.GT.0 .AND. IBIF.EQ.2) WRITE(UDIAG,20100) IPIVS
      IF (UDIAG.GT.0 .AND. IBIF.GT.0) WRITE(UDIAG,20000) YH,DYTA,YW
      INIT=1
      GOTO (5990,9000),IBIF
C.......................................................................
C  NEW BIFURCATION POINT
C  COMPUTATION OF TANGENTS
5400  CONTINUE
      CALL A2BIFT(N,N1,INFO,V,W,IPIV,A,B,W1,W2,EPS)
      IF (INFO(9).EQ.0) GOTO 5410
C  POSSIBLY NON-SIMPLE BIFURCATION POINT DETECTED
      IF (KPRINT.GT.0) WRITE(UPR,60609)
      GOTO 5990
C  ORDERING OF TANGENTS
5410  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60601) IND
      IND=IND+1
      IF (IND.GT.INDMAX) GOTO 9960
      INDEX(IND)=2
Crr      WRITE(UBIF,10001) YH,YW,ETA
      CALL A2RWR(N1,YH,YW,ETA)
      INFO(1)=0
      SK1=0.D0
      SK2=0.D0
      VN=0.D0
      DO 5420 I=1,N1
      S=(Y(I)-YH(I))/YW(I)
      SK1=SK1+S*W1(I)
      SK2=SK2+S*W2(I)
      VN=VN+S**2
5420  CONTINUE
      VN=DSQRT(VN)
      IF (DABS(SK1).GE.DABS(SK2)) GOTO 5431
      S=SK1
      SK1=SK2
      SK2=S
      DO 5430 I=1,N1
      S=W1(I)
      W1(I)=W2(I)
      W2(I)=S
5430  CONTINUE
5431  CONTINUE
      IF (SK1.LE.0.D0) GOTO 5441
      SK1=-SK1
      DO 5440 I=1,N1
      W1(I)=-W1(I)
5440  CONTINUE
5441  CONTINUE
C  REJECT COMPUTED BIFURCATION POINT
C  IN CASE OF UNRELIABLE TANGENT ORDERING
      IF (SK1/VN.LT.-0.7D0) GOTO 5500
      IF (KPRINT.GT.0) WRITE(UPR,60607)
      GOTO 5200
C  PREPARE START-OFF PROCEDURE
5500  CONTINUE
      IF (KPRINT.GT.1) WRITE(UPR,60602) W1
      IF (KPRINT.GT.1) WRITE(UPR,60603) W2
      IF (KPRINT.GT.1) WRITE(UPR,60001)
      IF (UDIAG.GT.0) WRITE(UDIAG,20101) IPIV(N1)
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) YH,W1,YW
      IBIF=1
      ICONT=0
      SIGMAH=0.5D0*SIGMA
C.......................................................................
C  FIXING OF LEAVING BRANCHES
C  PREDICT VALUE
5510  CONTINUE
      DO 5511 I=1,N1
      YD(I)=YH(I)-SIGMAH*W2(I)*YW(I)
5511  CONTINUE
      GOTO 5590
5520  CONTINUE
      DO 5521 I=1,N1
      YD(I)=YH(I)+SIGMAH*W2(I)*YW(I)
5521  CONTINUE
      GOTO 5590
5530  CONTINUE
      DO 5531 I=1,N1
      YD(I)=YH(I)+SIGMAH*W1(I)*YW(I)
5531  CONTINUE
C  CORRECTOR
5590  CONTINUE
      ICONT=ICONT+1
      IF (ICONT.GT.20) GOTO 9950
      FMIN=0.D0
      IF (YD(N1).LE.TAUMIN .OR. YD(N1).GE.TAUMAX) GOTO 5800
      DO 5591 I=1,N1
      YDH(I)=YD(I)
5591  CONTINUE
      INFO(8)=2
      CALL A2GNHM(N,N1,YD,FY,JFY,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,COND)
      IERR=INFO(9)
      ITS=INFO(7)
      FMIN=EPMACH*COND
      FMAX=1.D0/DSQRT(THR)
      IF (ITS.NE.0) FMAX=THMAX/DSQRT(SUMYQH/SUMYH)
      IF (IERR.NE.0 .OR. ITS.EQ.0) GOTO 5594
      VAL1=0.D0
      DO 5592 I=1,N1
      VAL1=VAL1+DYT(I)*V(I)
5592  CONTINUE
      VAL1=0.5D0*THMAX*DSQRT(SUMYH)/DABS(VAL1)
      TH=DSQRT(SUMYQH/SUMYH)
      VAL2=0.D0
      DO 5593 I=1,N1
      VAL2=VAL2+((YDH(I)-YD(I))/YW(I))**2
5593  CONTINUE
      VAL2=TH*DSQRT(VAL2)
5594  CONTINUE
      CALL A2TAND(N,N1,JFY,YD,YW,DFY,A,COND,DET,IPIV,AH,V,
     &   D,IRANK)
      IPIVS=IPIV(N1)
      IF (IRANK.LT.N) IERR=3
C  STEPLENGTH CONTROL IN CASE OF FAILURE
      IF (IERR.EQ.0) GOTO 5600
      IF (IERR.EQ.4) GOTO 5800
      IF (FMIN.GT.FMAX) GOTO 9950
      SIGMAH=DMAX1(0.5D0*FMAX,FMIN)*SIGMAH
      GOTO (5510,5520,5530),IBIF
C.......................................................................
C  PROJECTIONS ON TANGENTS OF BIFURCTION POINT
5600  CONTINUE
      SK1=0.D0
      SK2=0.D0
      VN=0.D0
      DO 5610 I=1,N1
      S=(YH(I)-YD(I))/YW(I)
      SK1=SK1+S*W1(I)
      SK2=SK2+S*W2(I)
      VN=VN+S**2
5610  CONTINUE
      VN=DSQRT(VN)
      SK1=SK1/VN
      SK2=SK2/VN
C  DETERMINANTS
      DETHMQ=DET
      IF (IPIVS.NE.N1) CALL A2DETH(N,N1,A,AH,D,N1,IPIV,DETHMQ)
C  NEW NORMALIZED TANGENT
      DO 5620 I=1,N
      W(IPIV(I))=-V(I)/D(N1)
5620  CONTINUE
      W(IPIVS)=1.D0/D(N1)
C  CHECK SIGN OF TANGENT
      SIGNUM=DSIGN(1.0D0,YD(IPIVS)-YH(IPIVS))
      IF (W(IPIVS)*SIGNUM.GT.0.D0) GOTO 5631
      DO 5630 I=1,N1
      W(I)=-W(I)
5630  CONTINUE
5631  CONTINUE
C.......................................................................
C  CHECK BRANCH
      GOTO (5710,5720,5730),IBIF
5710  CONTINUE
      IF (.NOT.(DABS(SK1).LE.DABS(SK2) .AND. SK2.GE.0.7D0)) GOTO 5800
      IF (DETHMA*DETHMQ*DYTA(N1)*W(N1).LT.0.D0) GOTO 5800
      DO 5711 I=1,N1
      V(I)=-W2(I)
5711  CONTINUE
Crr      WRITE(UBIF,10001) YD,V,(W(I),I=1,N1)
Crr      WRITE(UBIF,10002) IPIVS,SIGMAH,DET,DETHMQ
Crr      WRITE(UBIF,10002) ITS,VAL1,VAL2
      CALL A2RWR(N1,YD,V,W)
      CALL A2RWR4(IPIVS,SIGMAH,DET,DETHMQ)
      CALL A2RWR3(ITS,VAL1,VAL2)
      GOTO 5790
5720  CONTINUE
      IF (.NOT.(DABS(SK1).LE.DABS(SK2) .AND. SK2.LE.-0.7D0)) GOTO 5800
      IF (DETHMA*DETHMQ*DYTA(N1)*W(N1).LT.0.D0) GOTO 5800
Crr      WRITE(UBIF,10001) YD,W2,(W(I),I=1,N1)
Crr      WRITE(UBIF,10002) IPIVS,SIGMAH,DET,DETHMQ
Crr      WRITE(UBIF,10002) ITS,VAL1,VAL2
      CALL A2RWR(N1,YD,W2,W)
      CALL A2RWR4(IPIVS,SIGMAH,DET,DETHMQ)
      CALL A2RWR3(ITS,VAL1,VAL2)
      GOTO 5790
5730  CONTINUE
      IF (.NOT.(DABS(SK1).GE.DABS(SK2) .AND. SK1.LE.-0.7D0)) GOTO 5800
      IF (DETHMA*DETHMQ*DYTA(N1)*W(N1).GT.0.D0) GOTO 5800
      IF (DABS(W1(N1)).LT.EPS .AND. DETHMA*DETHMQ.LT.0.D0) GOTO 5800
5790  CONTINUE
      IF (KPRINT.GT.1) WRITE(UPR,60604) IBIF,SIGMAH
      IBIF=IBIF+1
      ICONT=0
      IF (IBIF.EQ.3) SIGMAH=0.5D0*SIGMA
      GOTO (5510,5520,5530,5900),IBIF
C.......................................................................
C  STEPLENGTH REDUCTION
5800  CONTINUE
      IF (FMIN.GE.1.D0 .OR. SIGMAH.LT.EPS) GOTO 9950
      SIGMAH=DMAX1(FMIN,0.5D0)*SIGMAH
      GOTO (5510,5520,5530),IBIF
C.......................................................................
C  EXIT FOR HANDLING OF BIFURCATION POINTS
5900  CONTINUE
      IF (KPRINT.GT.1) WRITE(UPR,60001)
      INFO(1)=KPRINT
      INIT=1
      DO 5910 I=1,N1
      Y(I)=YH(I)
      YQ(I)=YD(I)
      DYTA(I)=W1(I)
      DYT(I)=W(I)
5910  CONTINUE
      DETF=DET
      DETFA=0.D0
      DETHM=DETHMQ
      SIGMA=SIGMAH
5990  CONTINUE
C-----------------------------------------------------------------------
C  FINISH OUTPUT FOR CURRENT CONTINUATION STEP
      IF (KPRINT.GT.0) WRITE(UPR,60500) YQ(N1),(YQ(I),I=1,N)
      IF (UDIAG.GT.0) WRITE(UDIAG,20100) IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) YQ,DYT,YW
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      UNRM=0.D0
      IF (UMAX.EQ.0.D0) GOTO 6011
      DO 6010 I=1,N
      UNRM=UNRM+YQ(I)**2
6010  CONTINUE
      UNRM=DSQRT(UNRM)
6011  CONTINUE
      IF (INIT.EQ.2 .OR. UNRM.GT.UMAX) GOTO 9010
C  EXTRAPOLATED STEPLENGTH BOUND
      IHALT=0
      IF (DETFA*DETFH.LE.0.D0 .OR. DABS(DETFA).LE.DABS(DETFH) .OR.
     &   IEXTR.NE.1) GOTO 6020
      IHALT=1
      S=DETFH/DETFA
      YHALT=(YQ(IPIVA)-S*Y(IPIVA))/(1.D0-S)
6020  CONTINUE
      DETFA=DETF
      DETHMA=DETHM
      GOTO 2000
C-----------------------------------------------------------------------
C  LOOK FOR NEW OFFSET
9000  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
9010  CONTINUE
      DO 9020 L=1,IND
      LQ=L
      IF (INDEX(L).GT.0) GOTO 9100
9020  CONTINUE
      GOTO 9900
C.......................................................................
C  READ IN NEW OFFSET
9100  CONTINUE
      INIT=1
Crr      REWIND UBIF
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
      CALL A2RREW
      CALL A2RRD(N1,YH,V,W)
      ICONT=0
C  CHECK IF CYCLE HAS BEEN FINISHED AT STARTING POINT
      IF (TAUH.EQ.TAUMIN .OR. TAUH.EQ.TAUMAX) GOTO 9140
      IF (YQ(N1).NE.TAUH) GOTO 9140
      VN=0.D0
      DO 9110 I=1,N1
      VN=VN+((YH(I)-YQ(I))/YW(I))**2
9110  CONTINUE
      VN=DSQRT(VN)
      IF (VN.LE.10.D0*EPS) GOTO 9120
      ICONT=1
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
      GOTO 9173
C  CYCLE FINISHED AT STARTING POINT  LOOK FOR NEXT OFFSET
9120  CONTINUE
      INDEX(1)=0
      DO 9130 L=1,IND
      LQ=L
      IF (INDEX(L).GT.0) GOTO 9160
9130  CONTINUE
      GOTO 9900
9140  CONTINUE
      DO 9150 I=1,N1
      Y(I)=YH(I)
      YW(I)=V(I)
      ETA(I)=W(I)
9150  CONTINUE
9160  CONTINUE
Crr      READ(UBIF,10001) YQ,DYTA,DYT
Crr      READ(UBIF,10002) IPIVA,SIGMA,DETHMA,DETFA
      CALL A2RRD(N1,YQ,DYTA,DYT)
      CALL A2RRD4(IPIVA,SIGMA,DETHMA,DETFA)
      INDEX(LQ)=INDEX(LQ)-1
      IF (LQ.LT.2) GOTO 9173
      DETFA=0.D0
      DETHMA=0.D0
      DO 9172 L=2,LQ
Crr      READ(UBIF,10001) Y,YW,ETA
Crr      READ(UBIF,10001) YQ,DYTA,DYT
Crr      READ(UBIF,10002) IPIVS,SIGMA,DETF,DETHM
Crr      READ(UBIF,10002) ITS,VAL1,VAL2
      CALL A2RRD(N1,Y,YW,ETA)
      CALL A2RRD(N1,YQ,DYTA,DYT)
      CALL A2RRD4(IPIVS,SIGMA,DETF,DETHM)
      CALL A2RRD3(ITS,VAL1,VAL2)
      IF (INDEX(LQ).EQ.0) GOTO 9171
Crr      READ(UBIF,10001) YQ,DYTA,DYT
Crr      READ(UBIF,10002) IPIVS,SIGMA,DETF,DETHM
Crr      READ(UBIF,10002) ITS,VAL1,VAL2
      CALL A2RRD(N1,YQ,DYTA,DYT)
      CALL A2RRD4(IPIVS,SIGMA,DETF,DETHM)
      CALL A2RRD3(ITS,VAL1,VAL2)
      GOTO 9172
9171  CONTINUE
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
9172  CONTINUE
9173  CONTINUE
      LQ=LQ+1
      IF (LQ.GT.IND) GOTO 9181
      DO 9180 L=LQ,IND
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10001) YH,V,(W(I),I=1,N1)
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY,DUMMY
Crr      READ(UBIF,10002) IDUMMY,DUMMY,DUMMY
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
      CALL A2RRD(N1,YH,V,W)
      CALL A2RRD4(IDUMMY,DUMMY1,DUMMY2,DUMMY3)
      CALL A2RRD3(IDUMMY,DUMMY1,DUMMY2)
9180  CONTINUE
9181  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60203)
      LQ=LQ-2
      IF (KPRINT.GT.0 .AND. LQ.NE.0) WRITE(UPR,60601) LQ
      IHALT=0
      IF (ICONT.EQ.1) GOTO 5990
      IF (UDIAG.GT.0) WRITE(UDIAG,20102) IPIVS
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) Y,DYTA,YW
      IF (LQ.NE.0) GOTO 9190
      IF (KPRINT.GT.0) WRITE(UPR,60500) Y(N1),(Y(I),I=1,N)
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      GOTO 1100
9190  CONTINUE
      DO 9191 I=1,N1
      YD(I)=YQ(I)
      DYTH(I)=DYTA(I)
9191  CONTINUE
      GOTO 5990
C.......................................................................
C  ALL DONE
9900  CONTINUE
      INFO(5)=0
      IF (KPRINT.GT.0 .AND. IWARN.GT.0) WRITE(UPR,60909) IWARN
      IF (KPRINT.GT.0) WRITE(UPR,60900) IFCTEV,IJACEV
      DO 9901 I=1,N1
      Y(I)=YQ(I)
9901  CONTINUE
      RETURN
C.......................................................................
C  FAIL EXIT
C  MORE THAN ITMAX STEPS
9910  CONTINUE
      INFO(5)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      RETURN
C  STEPLENGTH TOO SMALL
9920  CONTINUE
      INFO(5)=2
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      IF (INFO(3).NE.0) RETURN
      IWARN=IWARN+1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      GOTO 9010
C  INITIAL VALUES TOO BAD
9940  CONTINUE
      INFO(5)=4
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60904)
      RETURN
C  START-OFF PROCEDURE FAILED
9950  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60905) IBIF,SIGMAH,FMIN,FMAX
      IND=IND-1
      INFO(1)=KPRINT
      IWARN=IWARN+1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      GOTO 9010
C  INTEGER WORKSPACE TOO SMALL
9960  CONTINUE
      INFO(5)=6
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60906)
      RETURN
C-----------------------------------------------------------------------
Crr 10001 FORMAT(4D18.10)
Crr 10002 FORMAT(I18,3D18.10)
20000 FORMAT(4D18.10)
20100 FORMAT('0',I4)
20101 FORMAT('1',I4)
20102 FORMAT('2',I4)
20103 FORMAT('3',I4)
60001 FORMAT(/)
60002 FORMAT(///)
60100 FORMAT(' ',132('*'))
60201 FORMAT(' STEP   0       CHECK FOR CONSISTENCY OF INITIAL VALUES'/)
60202 FORMAT(' STEP',I4/)
60203 FORMAT('                NEW OFFSET'/)
60205 FORMAT(' SCALED STEPLENGTH',D11.4/
     &   ' CURRENT CONTINUATION PARAMETER',I4/)
60206 FORMAT(' STEPLENGTH REDUCED BY EXTRAPOLATED STEPLENGTH BOUND')
60300 FORMAT(/' GAUSS-NEWTON METHOD REQUIRED',I3,' ITERATIONS'//)
60401 FORMAT(/' GAUSS-NEWTON METHOD FAILED'/' ',132('.')/)
60402 FORMAT(/' CONFLICT IN SIGNS'/' ',132('.')/)
60403 FORMAT(/' STEP RETRIED WITH REDUCED STEPLENGTH'/' ',132('.')/)
60500 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60501 FORMAT(/' DET WRT ACTUAL CON PARM',2D17.7/
     &   ' DET WRT  REAL  CON PARM',2D17.7/)
60601 FORMAT(' BIFURCATION POINT NO.',I3/)
60602 FORMAT(' TANGENT 1=',5(D18.11,6X)/9(11X,5(D18.11,6X)/))
60603 FORMAT(' TANGENT 2=',5(D18.11,6X)/9(11X,5(D18.11,6X)/))
60604 FORMAT(' START-OFF ON BRANCH',I2,' WITH STEPSIZE',D10.3)
60607 FORMAT(' UNRELIABLE TANGENT ASSIGNMENT'/)
60608 FORMAT(' COMPLETELY INVESTIGATED BIFURCATION POINT REACHED'/
     &   ' POSSIBLE REASON BRANCH SWITCHING DURING PATH FOLLOWING OR',
     &   ' OVERLOOKING BIFURCATION POINT'/
     &   ' BIFURCATION DIAGRAM MAY BE WRONG OR INCOMPLETE'/)
60609 FORMAT(' POSSIBLY NON-SIMPLE BIFURCATION POINT DETECTED'/
     &   ' COMPUTATION IGNORED'/)
60900 FORMAT(' ALCON2 REQUIRED',I5,' FUNCTION EVALUATIONS AND',I5,
     &   ' JACOBIAN EVALUATIONS'//)
60901 FORMAT(' TERMINATION AFTER INFO(2) CONTINUATION STEPS'//)
60902 FORMAT(' TERMINATION SINCE STEPLENGTH TOO SMALL'//)
60904 FORMAT(' TERMINATION SINCE INITIAL GUESS TOO FAR AWAY FROM',
     &   ' A REGULAR POINT OF THE CONTINUATION PATH'/
     &   ' USE NONLINEAR EQUATION SOLVER FOR BETTER INITIAL DATA'//)
60905 FORMAT(' TERMINATION SINCE START-OFF ALGORITHM FAILED AT BRANCH',
     &   I2/' LATEST USED STEPSIZE',D10.3,3X,
     &   ' CONVERGENCE INTERVAL (',D9.3,',',D9.3,')'//)
60906 FORMAT(' TERMINATION SINCE INTEGER WORKSPACE TOO SMALL'//)
60909 FORMAT(' ALCON2 HAS ISSUED',I3,' WARNING MESSAGES'//)
      END
C
C
      SUBROUTINE A2GNHM (N,N1,Y,FY,JFY,EPS,INFO,SUMYH,SUMYQH,
     &  TAUMIN,TAUMAX,PIVOT,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,U,
     &  COND)
C
C  SUBROUTINE  A2GNHM  TO BE USED WITH ROUTINE ALCON2
C
C  GAUSS-NEWTON METHOD AS CORRECTOR IN PATHFOLLOWING PROCEDURE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),V(N1),D(N1),YW(N1),ETA(N1)
      DOUBLE PRECISION DFY(N,N1),A(N,N1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(N1),DYQ(N1),W(N1)
      DOUBLE PRECISION F(N),FH(N),U(N)
      INTEGER PIVOT(N1),INFO(9)
      INTEGER UPR,UDIAG
      EXTERNAL FY,JFY
      COMMON /COUNT/ IFCTEV,IJACEV,IDUMQ,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A2UNIT/ UPR,UDIAG
      COMMON /A2VALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      COMMON /A2METH/ FACTOR,IEXTR,IJACM
C  INTERNAL PARAMETERS
      SIGMA=0.5D0/THMAX
C  INITIAL PREPARATIONS TO START A2GNHM
      SUMY=0.D0
      THMAX2=1.D0
      IF (INFO(8).NE.1) THMAX2=THMAX**2
      KPRINT=INFO(1)
      ITMAX=INFO(6)
      ITER=0
      LEVEL=0
      NEW=0
      NQ=N
      IF (INFO(8).EQ.2) NQ=N1
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y
C-----------------------------------------------------------------------
C  COMPUTATION OF RESIDUAL VECTOR
1     CALL FY(Y,Y(N1),F)
      IFCTEV=IFCTEV+1
      IF (LEVEL.EQ.1 .AND. IJACM.EQ.3 .AND. NEW.GT.0) GOTO 4
      IF (LEVEL.EQ.1) GOTO 43
C-----------------------------------------------------------------------
C  SCALED JACOBIAN MATRIX
2     NEW=0
      CALL JFY(Y,Y(N1),DFY,DFY(1,N1))
      IJACEV=IJACEV+1
      DO 21 K=1,N1
      S=YW(K)
      DO 21 I=1,N
21    DFY(I,K)=DFY(I,K)*S
      GOTO 4
C-----------------------------------------------------------------------
C  RANK-1 APPROXIMATION OF THE JACOBIAN DUE TO BROYDEN
C  (SCALED VERSION)
3     NEW=NEW+1
301   ST=1.D0/SUMY
      DO 31 I=1,N
      T=F(I)*ST
      DO 311 K=1,NQ
311   IF (DFY(I,K).NE.0.D0) DFY(I,K)=DFY(I,K)+T*DY(K)
31    CONTINUE
      IF (IJACM.EQ.3) GOTO 434
C-----------------------------------------------------------------------
C  SOLUTION OF THE LINEAR SYSTEM
4     DO 41 K=1,NQ
      DO 411 I=1,N
411   A(I,K)=-DFY(I,K)
41    CONTINUE
C  HOUSEHOLDER TRIANGULARIZATION
      COND=1.D0/EPMACH
      IRANK=N
      CALL A1DECC(A,N,N1,0,N,NQ,IRANK,COND,D,PIVOT,LEVEL,AH,W)
      IDECS=IDECS+1
      D1=DABS(D(1))
      IF (IRANK.LT.N) GOTO 93
      IF (ITER.NE.0 .OR. LEVEL.EQ.1 .OR. INFO(8).NE.2) GOTO 42
      DO 421 I=1,N
421   V(PIVOT(I))=-W(I)/D(N1)
      V(PIVOT(N1))=1.D0/D(N1)
42    CONTINUE
C  LINEAR LEAST SQUARES SOLUTION
43    DO 431 K=1,N
431   U(K)=F(K)
      CALL A1SOLC(A,N,N1,0,N,NQ,DYQ,U,IRANK,D,PIVOT,LEVEL,AH,W)
      ISOLS=ISOLS+1
      IF (LEVEL.EQ.1) GOTO 44
432   CONV=0.D0
      DO 433 L=1,NQ
      T=DYQ(L)
      CONV=DMAX1(CONV,DABS(T))
      YA(L)=Y(L)
433   Y(L)=Y(L)+T*YW(L)
C  TEST OF ACCURACY
      IF (CONV.LE.EPS) GOTO 9
      IF (IJACM.EQ.3 .AND. NEW.GT.0) GOTO 301
C  ORDINARY GAUSS-NEWTON CORRECTION
434   SUMY=0.D0
      DO 435 L=1,NQ
      T=DYQ(L)
      SUMY=SUMY+T*T
435   DY(L)=T
      LEVEL=1
      GOTO 1
C  SIMPLIFIED GAUSS-NEWTON CORRECTION
44    SUMYQ=0.D0
      DO 441 L=1,NQ
      T=DYQ(L)
441   SUMYQ=SUMYQ+T*T
C-----------------------------------------------------------------------
C  RESTRICTED NATURAL MONOTONICITY TEST
      SUMYQ=DMAX1(SMALL*SUMY,SUMYQ)
      IF (ITER.NE.0) GOTO 51
      SUMYH=SUMY
      SUMYQH=SUMYQ
51    CONTINUE
      IF (SUMYQ.GT.THMAX2*SUMY) GOTO 69
C  A-POSTERIORI ESTIMATE OF RELAXATION FACTOR
      FCH=0.5D0*DSQRT(SUMY/SUMYQ)
C-----------------------------------------------------------------------
C  PREPARATIONS TO START THE FOLLOWING ITERATION STEP
      LEVEL=0
      ITER=ITER+1
      IF (KPRINT.GT.1 .AND. ITER.EQ.1) WRITE(UPR,60003)
      IF (KPRINT.GT.1) WRITE(UPR,60002)
     &   ITER,SUMY,SUMYQ,NEW,IRANK,COND,D1,(PIVOT(I),I=1,NQ)
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y
      IF (Y(N1).LT.TAUMIN .OR. Y(N1).GT.TAUMAX) GOTO 94
      IF (ITER.GE.ITMAX) GOTO 91
      IF (FCH.LT.SIGMA .OR. IRANK.LT.N .OR. IJACM.EQ.2) GOTO 2
      IF (IJACM.EQ.1) GOTO 3
      IF (IJACM.NE.3) GOTO 63
      ALPHA=0.D0
      DO 61 L=1,NQ
61    ALPHA=ALPHA+DY(L)*DYQ(L)
      ALPHA=ALPHA/SUMY
      T=1.D0/(1.D0-ALPHA)
      DO 62 L=1,NQ
62    DYQ(L)=T*DYQ(L)
      NEW=NEW+1
      GOTO 432
63    NEW=NEW-1
      GOTO 43
69    CONTINUE
      IF (INFO(8).NE.1 .OR. NEW.EQ.0) GOTO 92
      DO 691 I=1,NQ
691   Y(I)=YA(I)
      LEVEL=0
      GOTO 1
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9     INFO(7)=ITER
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  TERMINATION AFTER ITMAX ITERATIONS
91    INFO(9)=1
      IF (KPRINT.GT.1) WRITE(UPR,60901)
      RETURN
C  RESTRICTED MONOTONICITY TEST VIOLATED
92    INFO(9)=2
      IF (KPRINT.GT.1) WRITE(UPR,60902)
      SUMYH=SUMY
      SUMYQH=SUMYQ
      RETURN
C  RANK DEFICIENCY OF JACOBIAN
93    INFO(9)=3
      IF (KPRINT.GT.1) WRITE(UPR,60903)
      RETURN
C  TAUMIN OR TAUMAX CROSSED
94    INFO(9)=4
      IF (KPRINT.GT.1) WRITE(UPR,60904)
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60002 FORMAT(5X,I2,4X,D10.4,3X,D10.4,4X,I2,4X,I2,
     &   4X,D8.2,3X,D8.2,3X,20I3)
60003 FORMAT(4X,'ITER',5X,'LEVELX',6X,'LEVELXQ',6X,'NEW',2X,'RANK',
     &   5X,'COND',7X,'SENS',9X,'PIVOTS')
60901 FORMAT(' TERMINATION AFTER ITMAX ITERATIONS')
60902 FORMAT(' TERMINATION SINCE RESTRICTED MONOTONICITY TEST WAS ',
     &   'VIOLATED')
60903 FORMAT(' RANK DEFICIENT JACOBIAN')
60904 FORMAT(' TERMINATION SINCE TAUMIN OR TAUMAX WAS CROSSED')
      END
C
C
      SUBROUTINE A2NEWA (N,N1,NH,NLT,Y,FY,JFY,EPS,INFO,TAUMIN,TAUMAX,
     &   PIVOT,LT,V,D,YW,ETA,DFY,A,AH,YA,DY,DYQ,W,F,FH,
     &   U,Z,ZA,ZH,DFYH,B,W1,W2)
C
C  SUBROUTINE  A2NEWA  TO BE USED WITH ROUTINE ALCON2
C
C  NEWTON METHOD FOR AUGMENTED SYSTEM OF MOORE
C  (SYMMETRIC IMPLEMENTATION)
C  COMPUTATION OF SIMPLE (POSSIBLY IMPERFECT) BIFURCATION POINTS
C
C-
C- CHANGED: CALLS OF HARWELL ROUTINES MA29BD, MA29CD SUBSTITUTED
C-          BY CALLS OF LINPACK ROUTINES SSPFA AND SSPSL .
C- ADAPTED BY: Lutz Weimann    11.12.86
C-
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),V(N1),D(N1),YW(N1),ETA(N1)
      DOUBLE PRECISION DFY(N,N1),A(1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(NH),DYQ(NH),W(1)
      DOUBLE PRECISION F(NH),FH(N),U(NH)
      DOUBLE PRECISION Z(N),ZA(N),ZH(N),DFYH(N,N1)
      DOUBLE PRECISION B(1),W1(N1),W2(N1)
      INTEGER PIVOT(N1),INFO(9),LT(NLT)
      INTEGER UPR,UDIAG
      EXTERNAL FY,JFY
      COMMON /COUNT/ IFCTEV,IJACEV,IDUM,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A2UNIT/ UPR,UDIAG
C-    COMMON /MA29DD/ GR,IRANK,LP
C-    LOC(I,J)=I+(2*NE-J)*(J-1)/2
      LOC(I,J)=I*(I-1)/2+J
C  INTERNAL PARAMETERS
      NE=2*N
      NQ=NE*(NE+1)/2
      LP=0
C  INITIAL PREPARATIONS TO START A2NEWA
      KPRINT=INFO(1)
      ITMAX=INFO(6)
      ITER=0
      LEVEL=0
C-----------------------------------------------------------------------
C  COMPUTATION OF RESIDUAL VECTOR
1     CALL JFY(Y,Y(N1),DFYH,DFYH(1,N1))
      IJACEV=IJACEV+1
C  INITIAL VALUE FOR VARIABLE Z AND INITIAL DECOMPOSITION OF SCALED F'
      IF (LEVEL.NE.0) GOTO 11
      NM1=N-1
      IF (NM1.EQ.0) GOTO 112
      DO 111 I=1,NM1
111   Z(I)=0.D0
112   Z(N)=1.D0
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y,Z
      DO 113 K=1,N1
      S=YW(K)
      DO 113 I=1,N
113   DFY(I,K)=DFYH(I,K)*S
      IRANKH=N
      COND=1.D0/EPMACH
      CALL A1DECC(DFY,N,N1,0,N,N1,IRANKH,COND,D,PIVOT,0,AH,W)
      IDECS=IDECS+1
      D1=DABS(D(1))
11    CONTINUE
C  COMPUTATION OF Z
      DO 12  I=1,N
12    ZH(I)=Z(I)
      DO 13 JJ=1,IRANKH
      J=IRANKH+1-JJ
      S=0.D0
      DO 131 I=J,N
131   S=S+DFY(I,J)*ZH(I)
      S=S/(D(J)*DFY(J,J))
      DO 132 I=J,N
132   ZH(I)=ZH(I)+DFY(I,J)*S
13    CONTINUE
C  COMPUTATION OF F'TRANSPOSED*Z
      DO 14 I=1,N1
      S=0.D0
      DO 141 K=1,N
      S=S+DFYH(K,I)*ZH(K)
141   CONTINUE
      F(I)=S*YW(I)
14    W2(I)=S
C  STORING OF BETA3 AND BETA4
      DO 15 I=1,N1
15    U(I)=F(PIVOT(I))
      DO 16 I=1,N1
16    F(I)=U(I)
C  COMPUTATION OF Q*F
      CALL FY(Y,Y(N1),FH)
      IFCTEV=IFCTEV+1
      DO 17 I=1,N
17    W1(I)=FH(I)
      DO 18 J=1,IRANKH
      S=0.D0
      DO 181 I=J,N
      S=S+DFY(I,J)*W1(I)
181   CONTINUE
      S=S/(D(J)*DFY(J,J))
      DO 182 I=J,N
      W1(I)=W1(I)+DFY(I,J)*S
182   CONTINUE
18    CONTINUE
C  INITIAL ESTIMATE FOR ALPHA
      IF (ITER.EQ.0 .AND. LEVEL.EQ.0) ALPHA=W1(N)
C  STORING OF BETA1BAR
      IF (NM1.EQ.0) GOTO 191
      DO 19 I=1,NM1
      F(N1+I)=(W1(I)-ALPHA*Z(I))-(W1(N)-ALPHA*Z(N))*Z(I)/Z(N)
      U(N1+I)=F(N1+I)
19    CONTINUE
191   CONTINUE
      F(N1+N)=W1(N)
      U(N1+N)=F(N1+N)
      IF (LEVEL.EQ.1) GOTO 43
C-----------------------------------------------------------------------
C  REDUCED JACOBIAN OF AUGMENTED SYSTEM IN SCALED VARIABLES
2     NEW=0
      IF (ITER.EQ.0) GOTO 21
C  BACKPERMUTATION OF BETA3 AND BETA4
      DO 211 I=1,N1
211   U(PIVOT(I))=F(I)
C  DECOMPOSITION OF SCALED F' FOR NEW PERMUTATION AND
C  NEW ORTHOGONAL TRANSFORMATION Q
      DO 212 K=1,N1
      S=YW(K)
      DO 212 I=1,N
212   DFY(I,K)=DFYH(I,K)*S
      IRANKH=N
      COND=1.D0/EPMACH
      CALL A1DECC(DFY,N,N1,0,N,N1,IRANKH,COND,D,PIVOT,0,AH,W)
      IDECS=IDECS+1
      D1=DABS(D(1))
C  UPDATING OF BETA3 AND BETA4
      DO 213 I=1,N1
213   F(I)=U(PIVOT(I))
      DO 214 I=1,N1
214   U(I)=F(I)
C  COMPUTATION OF Z AND Q*F WITH NEW Q
      DO 215 I=1,N
215   Z(I)=ZH(I)
      DO 216 J=1,IRANKH
      S1=0.D0
      S2=0.D0
      DO 2161 I=J,N
      S1=S1+DFY(I,J)*FH(I)
      S2=S2+DFY(I,J)*Z(I)
2161  CONTINUE
      S1=S1/(D(J)*DFY(J,J))
      S2=S2/(D(J)*DFY(J,J))
      DO 2162 I=J,N
      FH(I)=FH(I)+DFY(I,J)*S1
      Z(I)=Z(I)+DFY(I,J)*S2
2162  CONTINUE
216   CONTINUE
C  UPDATING OF BETA1BAR DUE TO NEW Z AND Q*F
      IF (NM1.EQ.0) GOTO 2171
      DO 217 I=1,NM1
      F(N1+I)=(FH(I)-ALPHA*Z(I))-(FH(N)-ALPHA*Z(N))*Z(I)/Z(N)
      U(N1+I)=F(N1+I)
217   CONTINUE
2171  CONTINUE
      F(N1+N)=FH(N)
      U(N1+N)=F(N1+N)
C  REDUCED JACOBIAN OF AUGMENTED SYSTEM IN SYMMETRIC STORAGE MODE
21    CONTINUE
C  MATRICES R AND SBAR
      IF (NM1.EQ.0) GOTO 221
      DO 22 I=1,NM1
      B(LOC(N1+I,N))=DFY(I,N)-D(N)*Z(I)/Z(N)
      B(LOC(N1+I,N1))=DFY(I,N1)-DFY(N,N1)*Z(I)/Z(N)
      DO 22 K=1,NM1
      IF (I.LT.K) B(LOC(N1+I,K))=DFY(I,K)
      IF (I.EQ.K) B(LOC(N1+I,I))=D(I)
      IF (I.GT.K) B(LOC(N1+I,K))=0.D0
22    CONTINUE
221   CONTINUE
C  MATRIX T
      IF (NM1.EQ.0) GOTO 231
      DO 23 I=1,NM1
      DO 23 K=1,I
      S=-ALPHA*Z(I)*Z(K)/Z(N)**2
      IF (I.EQ.K) S=S-ALPHA
23    B(LOC(N1+I,N1+K))=S
231   CONTINUE
C  MATRIX C
      DO 24 I=1,N1
      S=0.D0
      DO 241 K=1,N
241   S=S+DFYH(K,I)*ZH(K)
24    W2(I)=S
      CALL A2DERI(JFY,N,N1,Y,DFYH,ZH,W2,W,YW,ETA,AH)
      DO 25 I=1,N1
      II=PIVOT(I)
      S=YW(II)
      DO 25 K=1,I
      KK=PIVOT(K)
25    B(LOC(I,K))=0.5D0*(AH(II,KK)+AH(KK,II))*YW(II)*YW(KK)
C-----------------------------------------------------------------------
C  SOLUTION OF THE LINEAR SYSTEM
      DO 41 I=1,NQ
      A(I)=-B(I)
41    CONTINUE
      DO 42 I=1,NE
      U(I)=F(I)
42    CONTINUE
C  RATIONAL CHOLESKY DECOMPOSITION WITH PIVOTING
C-    CALL MA29BD(NE,A,LT)
      CALL DSPFA(A,NE,LT,IRINFO)
      IF (IRINFO.NE.0) GOTO 93
C  SOLUTION OF TRIANGULAR SYSTEM
C-43    CALL MA29CD(NE,A,LT,U,W)
43    CALL DSPSL(A,NE,LT,U)
ctest+
c-      write(6,*) 'u:',(u(i),i=1,n1)
ctest-
      DO 431 I=1,N1
431   DYQ(PIVOT(I))=U(I)
      IF (NM1.EQ.0) GOTO 4321
      DO 432 I=1,NM1
432   DYQ(N1+I)=U(N1+I)
4321  CONTINUE
      IF (LEVEL.EQ.1) GOTO 44
C  ORDINARY NEWTON CORRECTION
      S=0.D0
      IF (NM1.EQ.0) GOTO 4331
      DO 433 I=1,NM1
433   S=S+Z(I)*DYQ(N1+I)
4331  CONTINUE
      DYQ(NE+1)=-S/Z(N)
      DYQ(NE+2)=(D(N)*U(N)+DFY(N,N1)*U(N1)-ALPHA*
     &   DYQ(NE+1)+(F(N1+N)-ALPHA*Z(N)))/Z(N)
      SUMY=0.D0
      CONV=0.D0
      DO 434 L=1,NH
      T=DYQ(L)
      SUMY=SUMY+T*T
      CONV=DMAX1(CONV,DABS(T))
434   DY(L)=T
      DO 435 L=1,N1
      T=DYQ(L)
      YA(L)=Y(L)
435   Y(L)=Y(L)+T*YW(L)
      S=0.D0
      DO 436 L=1,N
      T=DYQ(N1+L)
      ZA(L)=Z(L)
      Z(L)=Z(L)+T
436   S=S+Z(L)**2
      S=DSQRT(S)
      DO 437 L=1,N
437   Z(L)=Z(L)/S
      T=DYQ(NH)
      ALPHAA=ALPHA
      ALPHA=ALPHA+T
C  TEST OF ACCURACY
      IF (CONV.LE.EPS) GOTO 9
      LEVEL=1
      GOTO 1
C  SIMPLIFIED NEWTON CORRECTION
44    CONTINUE
      S=0.D0
      IF (NM1.EQ.0) GOTO 4411
      DO 441 I=1,NM1
441   S=S+ZA(I)*DYQ(N1+I)
4411  CONTINUE
      DYQ(NE+1)=-S/ZA(N)
      DYQ(NE+2)=(D(N)*U(N)+DFY(N,N1)*U(N1)-ALPHAA*
     &   DYQ(NE+1)+(F(N1+N)-ALPHA*Z(N)))/ZA(N)
      SUMYQ=0.D0
      DO 442 L=1,NH
      T=DYQ(L)
442   SUMYQ=SUMYQ+T*T
C-----------------------------------------------------------------------
C  NATURAL MONOTONICITY TEST
      IF (SUMYQ.GT.SUMY) GOTO 92
C-----------------------------------------------------------------------
C  PREPARATIONS TO START THE FOLLOWING ITERATION STEP
      LEVEL=0
      ITER=ITER+1
      IF (KPRINT.GT.1 .AND. ITER.EQ.1) WRITE(UPR,60003)
      IF (KPRINT.GT.1) WRITE(UPR,60002)
     &   ITER,SUMY,SUMYQ,NEW,IRANKH,COND,D1,(PIVOT(I),I=1,N1)
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y,Z
      IF (Y(N1).LT.TAUMIN .OR. Y(N1).GT.TAUMAX) GOTO 94
      IF  (ITER.GT.ITMAX) GOTO 91
      GOTO 2
C-----------------------------------------------------------------------
C  SOLUTION EXIT
C-9     IF (IRANK.LT.NE) GOTO 93
9     IF (IRINFO.NE.0) GOTO 93
      INFO(9)=0
      INFO(7)=ITER
      DO 901 I=1,N
901   Z(I)=ZH(I)
      IF (DABS(DYQ(NH)).LT.0.5D0*DABS(ALPHAA)) GOTO 902
      IF (KPRINT.GT.1) WRITE(UPR,60091)
      RETURN
902   CONTINUE
      IF (KPRINT.GT.1) WRITE(UPR,60092) ALPHA
      RETURN
C  FAIL EXIT
C  TERMINATION AFTER ITMAX ITERATIONS
91    INFO(9)=1
      IF (KPRINT.GT.1) WRITE(UPR,60901)
      RETURN
C  ITERATION DIVERGES
92    INFO(9)=2
      IF (KPRINT.GT.1) WRITE(UPR,60902)
      RETURN
C  RANK DEFICIENCY OF JACOBIAN
93    INFO(9)=3
      IF (KPRINT.GT.1) WRITE(UPR,60903)
      RETURN
C  TAUFIN CROSSED
94    INFO(9)=4
      IF (KPRINT.GT.1) WRITE(UPR,60904)
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60002 FORMAT(5X,I2,4X,D10.4,3X,D10.4,4X,I2,4X,I2,4X,
     &   D8.2,3X,D8.2,3X,20I3)
60003 FORMAT(4X,'ITER',5X,'LEVELX',6X,'LEVELXQ',6X,'NEW',2X,'RANK',
     &   5X,'COND',7X,'SENS',9X,'PIVOTS')
60091 FORMAT(' NUMERICALLY PERFECT BIFURCATION'/)
60092 FORMAT(' NUMERICALLY IMPERFECT BIFURCATION'/
     &   ' IMPERFECTION PARAMETER IS',D12.3/)
60901 FORMAT(' TERMINATION AFTER ITMAX ITERATIONS')
60902 FORMAT(' TERMINATION SINCE ITERATION DIVERGES')
C-60903 FORMAT(' RANK DEFICIENT SOLUTION OBTAINED')
60903 FORMAT(' TERMINATION SINCE RANK DEFICIENT SOLUTION OBTAINED')
60904 FORMAT(' TERMINATION SINCE TAUFIN WAS CROSSED')
      END
C
C
      SUBROUTINE A2TURN (N,N1,IPIVA,EPS,INFO,
     &   Y,Y0,Y1,YH,DYT,DYT0,DYT1,V,YW,
     &   IPIV,FY,JFY,TAUMIN,TAUMAX,D,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U)
C
C  SUBROUTINE  A2TURN  TO BE USED WITH ROUTINE ALCON2
C
C  DETERMINATION OF TURNING POINTS
C  INTERVAL METHOD WITH CUBIC HERMITE INTERPOLATION AS PREDICTOR
C  AND GAUSS-NEWTON METHOD AS CORRECTOR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),Y0(N1),Y1(N1),YH(N1),DYT(N1),DYT0(N1)
      DOUBLE PRECISION DYT1(N1),V(N1),YW(N1)
      DOUBLE PRECISION D(N1),ETA(N1),DFY(N,N1),AQ(1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(N1),DYQ(N1),W(N1)
      DOUBLE PRECISION F(N),FH(N),U(N)
      INTEGER INFO(9),IPIV(N1)
      INTEGER UPR,UDIAG
      EXTERNAL FY,JFY
      COMMON /A2UNIT/ UPR,UDIAG
      COMMON /A2VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C  PREPARATIONS
      KPRINT=INFO(1)
      INFO(9)=-1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INFO(4).EQ.0 .OR. INFO(4).EQ.2) GOTO 9990
      INFO(1)=0
      IPIVH=IPIVA
      IF (IPIVH.EQ.N1) IPIVH=IPIV(N1)
      DO 1000 I=1,N1
      Y0(I)=Y(I)
      DYT1(I)=DYT(I)
1000  CONTINUE
      F0P=DYT0(N1)/DYT0(IPIVH)
      F1P=DYT1(N1)/DYT1(IPIVH)
      ITER=0
      X0=Y0(IPIVH)/YW(IPIVH)
      X1=Y1(IPIVH)/YW(IPIVH)
      H=X1-X0
      DEL=DABS(H)
C-----------------------------------------------------------------------
C  ITERATION LOOP
2000  CONTINUE
C  HERMITE INTERPOLATION
      F0=Y0(N1)/YW(N1)
      F1=Y1(N1)/YW(N1)
      A=6.D0*(F0-F1)+3.D0*H*(F0P+F1P)
      B=-6.D0*(F0-F1)-2.D0*H*(2.D0*F0P+F1P)
      C=H*F0P
      Z=B**2-4.D0*A*C
      Z=DSIGN(DSQRT(Z),-B)
      Z=(-B+Z)/(2.D0*A)
      IF (Z.LT.0.D0 .OR. Z.GT.1.D0) Z=(C/A)/Z
      YH(IPIVH)=Y0(IPIVH)+H*Z*YW(IPIVH)
      DO 2020 I=1,N1
      IF (I.EQ.IPIVH) GOTO 2010
      F0H=Y0(I)/YW(I)
      F1H=Y1(I)/YW(I)
      F0PH=DYT0(I)/DYT0(IPIVH)
      F1PH=DYT1(I)/DYT1(IPIVH)
      A=2.D0*(F0H-F1H)+H*(F0PH+F1PH)
      B=-3.D0*(F0H-F1H)-H*(2.D0*F0PH+F1PH)
      C=H*F0PH
      R=F0H+Z*(C+Z*(B+Z*A))
      YH(I)=R*YW(I)
2010  CONTINUE
2020  CONTINUE
C  ITERATION BACK TO CONTINUATION PATH
      ITER=ITER+1
      IF (ITER.GT.ITMAX) GOTO 9910
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      INFO(8)=2
      CALL A2GNHM(N,N1,YH,FY,JFY,EPS,INFO,SUMYH,SUMYQH,
     &  TAUMIN,TAUMAX,IPIV,V,D,YW,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U,COND)
      IF (INFO(9).NE.0) GOTO 9920
      CALL A2TAND(N,N1,JFY,YH,YW,DFY,AQ,COND,DETF,IPIV,AH,V,
     &   D,IRANK)
      IF (IRANK.LT.N) GOTO 9920
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      IF (DABS(YH(IPIVH)-Y0(IPIVH)).GT.DEL*YW(IPIVH) .OR.
     &    DABS(YH(IPIVH)-Y1(IPIVH)).GT.DEL*YW(IPIVH)) GOTO 9930
      DO 2030 I=1,N
      W(IPIV(I))=-V(I)/D(N1)
2030  CONTINUE
      W(IPIV(N1))=1.D0/D(N1)
      FQP=W(N1)/W(IPIVH)
C  NEW INCLUSION INTERVAL
      IF (F0P*FQP.GT.0.D0) GOTO 2200
      DQ=DABS(YH(IPIVH)-Y1(IPIVH))
      F1P=FQP
      DO 2110 I=1,N1
      Y1(I)=YH(I)
      DYT1(I)=W(I)
2110  CONTINUE
      GOTO 2300
2200  CONTINUE
      DQ=DABS(YH(IPIVH)-Y0(IPIVH))
      F0P=FQP
      DO 2210 I=1,N1
      Y0(I)=YH(I)
      DYT0(I)=W(I)
2210  CONTINUE
2300  CONTINUE
      X0=Y0(IPIVH)/YW(IPIVH)
      X1=Y1(IPIVH)/YW(IPIVH)
      HRR=X1-X0
      DEL=DABS(HRR)
      H = HRR
      IF (KPRINT.GT.0 .AND. ITER.EQ.1) WRITE(UPR,60101)
      IF (KPRINT.GT.0) WRITE(UPR,60102) ITER,YH(IPIVH),DEL,FQP,INFO(7)
      IF (FQP.EQ.0.D0) GOTO 9000
      IF (DEL.LT.EPS) GOTO 9000
      IF (DQ/YW(IPIVH).LT.EPS) GOTO 9000
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9000  CONTINUE
      INFO(9)=0
      IF (KPRINT.GT.0) WRITE(UPR,60103) ITER
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(N1),(YH(I),I=1,N)
      GOTO 9990
C  FAIL EXIT
C  MORE THAN ITMAX ITERATIONS
9910  CONTINUE
      INFO(9)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      GOTO 9990
C  GAUSS-NEWTON METHOD FAILED
9920  CONTINUE
      INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      GOTO 9990
C  ITERATION STRATEGY FAILED
9930  CONTINUE
      INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60903)
C  RETURN TO CALLING ROUTINE
9990  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      INFO(1)=KPRINT
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(' ',132('-'))
60002 FORMAT(' POSSIBLY TURNING POINT DETECTED')
60003 FORMAT(/)
60101 FORMAT(/4X,'IGNC',7X,'YHIT',8X,'DELTAX',8X,'DERIV',6X,'ITER')
60102 FORMAT(5X,I2,4X,D11.4,3X,D10.4,3X,D11.4,4X,I2)
60103 FORMAT(/' TURN REQUIRED',I3,' GAUSS-NEWTON-CALLS'/)
60104 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60201 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60901 FORMAT(/' TURN TERMINATED AFTER ITMAX ITERATIONS')
60902 FORMAT(/' TURN TERMINATED SINCE GAUSS-NEWTON METHOD FAILED')
60903 FORMAT(/' TURN TERMINATED SINCE ITERATION STRATEGY FAILED')
      END
C
C
      SUBROUTINE A2BIFU (N,N1,NH,NLT,IPIVQ,EPS,INFO,Y0,Y1,YH,
     &   V,YW,IPIV,LT,FY,JFY,TAUMIN,TAUMAX,D,ETA,DFY,AQ,AH,YA,
     &   DY,DYQ,W,F,FH,U,ZQ,ZA,ZH,DFYH,BQ,DETF,DETFQ,W1,W2)
C
C  SUBROUTINE  A2BIFU  TO BE USED WITH ROUTINE ALCON2
C
C  DETERMINATION OF (SIMPLE) BIFURCATION POINTS
C  NEWTON METHOD FOR AUGMENTED SYSTEM OF MOORE
C  DRIVER ROUTINE FOR SUBROUTINE A2NEWA
C
C-
C- CHANGED: EXIT LABEL 9940 NO LONGER REFERENCED DUE TO CHANGES OF CALLS
C-          OF EQUATION SOLVER ROUTINES IN SUBROUTINE A2NEWA.
C- ADAPTED BY: Lutz Weimann    11.12.86
C-
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y0(N1),Y1(N1),YH(N1)
      DOUBLE PRECISION V(N1),YW(N1)
      DOUBLE PRECISION D(N1),ETA(N1),DFY(N,N1),AQ(1),AH(N1,N1)
      DOUBLE PRECISION YA(N1),DY(NH),DYQ(NH),W(1)
      DOUBLE PRECISION F(NH),FH(N),U(NH)
      DOUBLE PRECISION ZQ(N),ZA(N),ZH(N),DFYH(N,N1),BQ(1)
      DOUBLE PRECISION W1(N1),W2(N1)
      INTEGER INFO(9),IPIV(N1),LT(NLT)
      INTEGER UPR,UDIAG
      COMMON /A2UNIT/ UPR,UDIAG
      EXTERNAL FY,JFY
C  PREPARATIONS
      KPRINT=INFO(1)
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INFO(3).NE.0 .AND. INFO(4).LT.2) GOTO 9990
      X0=Y0(IPIVQ)/YW(IPIVQ)
      X1=Y1(IPIVQ)/YW(IPIVQ)
      H=X1-X0
      DEL=DABS(H)
C-----------------------------------------------------------------------
C  INITIAL ESTIMATE BY LINEAR INTERPOLATION
      Z=-DETFQ/(DETF-DETFQ)
      DO 1000 I=1,N1
      F0H=Y0(I)/YW(I)
      F1H=Y1(I)/YW(I)
      R=F0H+Z*(F1H-F0H)
      YH(I)=R*YW(I)
1000  CONTINUE
C  NEWTON METHOD FOR AUGMENTED SYSTEM OF MOORE
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      CALL A2NEWA(N,N1,NH,NLT,YH,FY,JFY,EPS,INFO,TAUMIN,TAUMAX,
     &     IPIV,LT,V,D,YW,ETA,DFY,AQ,AH,YA,DY,DYQ,W,F,FH,U,
     &     ZQ,ZA,ZH,DFYH,BQ,W1,W2)
C-    IF (INFO(9).EQ.3) GOTO 9940
      IF (INFO(9).NE.0) GOTO 9920
      IF (INFO(8).EQ.2) GOTO 9000
      IF (DABS(YH(IPIVQ)-Y0(IPIVQ)).GT.DEL*YW(IPIVQ) .OR.
     &    DABS(YH(IPIVQ)-Y1(IPIVQ)).GT.DEL*YW(IPIVQ)) GOTO 9930
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9000  CONTINUE
      INFO(9)=0
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(N1),(YH(I),I=1,N)
      GOTO 9990
C  FAIL EXIT
C  NEWTON METHOD FOR AUGMENTED SYSTEM FAILED
9920  CONTINUE
      INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      GOTO 9990
C  SOLUTION OUT OF INTERVAL
9930  CONTINUE
      INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60903)
      GOTO 9990
C  RANK DEFICIENT AUGMENTED SYSTEM FOR SOLUTION
Crr 9940  CONTINUE
      INFO(9)=4
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(N1),(YH(I),I=1,N)
C  RETURN TO CALLING ROUTINE
9990  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(' ',132('-'))
60002 FORMAT(' POSSIBLY BIFURCATION POINT DETECTED')
60003 FORMAT(/)
60104 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60201 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60902 FORMAT(/' BIFURC TERMINATED SINCE NEWTON METHOD FAILED')
60903 FORMAT(/' BIFURC TERMINATED SINCE SOLUTION POINT OUT OF ',
     &   'PARAMETRIZATION INTERVAL')
      END
C
C
      SUBROUTINE A2DERI (JFY,N,N1,Y,DFYH,Z,V,VH,YW,ETA,C)
C
C  SUBROUTINE  A2DERI  TO BE USED WITH ROUTINE ALCON2
C
C  COMPUTATION OF (F'TRANSPOSED*Z)' BY FINITE DIFFERENCES
C  INCLUDING FEED BACK DEVICE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),V(N1),VH(N1),YW(N1),ETA(N1),C(N1,N1)
      DOUBLE PRECISION DFYH(N,N1),Z(N)
      EXTERNAL JFY
      COMMON /COUNT/ IFCTEV,IJACEV,ITER,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /A2VALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      EPDIFF=1.D3*EPMACH
      ETAD=DSQRT(1.D1*EPMACH)
      DO 1 K=1,N1
      K1=0
11    T=Y(K)
      ETAH=ETA(K)
      IF (T.LT.0.D0) ETAH=-ETAH
      S=ETAH*YW(K)
      Y(K)=T+S
      CALL JFY(Y,Y(N1),DFYH,DFYH(1,N1))
      IJACEV=IJACEV+1
      DO 12 I=1,N1
      SQ=0.D0
      DO 121 J=1,N
121   SQ=SQ+DFYH(J,I)*Z(J)
12    VH(I)=SQ
      Y(K)=T
      SU=0.D0
      DO 13 I=1,N1
      HG=DMAX1(DABS(V(I)),DABS(VH(I)))
      T=VH(I)-V(I)
      IF (HG.NE.0.D0) SU=SU+(T/HG)**2
13    C(I,K)=T/S
      SU=DSQRT(SU/DBLE(N1))
      IF ((SU.EQ.0.D0).OR.(K1.GT.0)) GOTO 1
      ETA(K)=DMAX1(EPDIFF,DSQRT(ETAD/SU)*ETA(K))
      ETA(K)=DMIN1(SMALL1,ETA(K))
      K1=1
      IF (SU.LT.EPDIFF) GOTO 11
1     CONTINUE
      RETURN
      END
C
C
      SUBROUTINE A2DETH (N,N1,A,AH,D,IPIVA,IPIV,DET)
C
C  SUBROUTINE  A2DETH  TO BE USED WITH ROUTINE ALCON2
C
C  COMPUTATION OF DETERMINANT OF A HESSENBERG MATRIX
C  BY GIVENS ROTATIONS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(N,N1),AH(N,N1),D(N1),IPIV(N1)
C  SEARCH OLD PIVOT
      DO 1010 I=1,N1
      IQ=I
      IF (IPIV(I).EQ.IPIVA) GOTO 1020
1010  CONTINUE
1020  CONTINUE
C  PRODUCT OF FIRST IQ-1 DIAGONAL ELEMENTS
      DET=1.D0
      IF (IQ.EQ.1) GOTO 2011
      IEND=IQ-1
      DO 2010 I=1,IEND
      DET=DET*D(I)
2010  CONTINUE
2011  CONTINUE
C  N-IQ GIVENS TRANSFORMATIONS
      IF (IQ.EQ.N1) GOTO 3050
      IBEG=IQ+1
      X1=A(IBEG-1,IBEG)
      IF (IBEG.GT.N) GOTO 3040
      DO 3010 I=IBEG,N1
      AH(IQ,I)=A(IQ,I)
3010  CONTINUE
      DO 3030 I=IBEG,N
      Y1=D(I)
      R=DSQRT(X1**2+Y1**2)
      C=X1/R
      S=Y1/R
      DET=DET*(C*X1+S*Y1)
      I1=I+1
      DO 3020 J=I1,N1
      X2=AH(I-1,J)
      Y2=A(I,J)
      AH(I,J)=-S*X2+C*Y2
3020  CONTINUE
      X1=AH(I,I1)
3030  CONTINUE
3040  CONTINUE
      DET=DET*X1
3050  CONTINUE
C  SIGN OF PIVOT VECTOR
      DO 4020 I=1,N
      IF (I.EQ.IQ) GOTO 4020
      I1=I+1
      DO 4010 J=I1,N1
      IF (J.EQ.IQ) GOTO 4010
      IF (IPIV(J).LT.IPIV(I)) DET=-DET
4010  CONTINUE
4020  CONTINUE
      RETURN
      END
C
C
      SUBROUTINE A2TAND (N,N1,JFY,Y,YW,DFY,A,COND,DET,PIVOT,AH,V,
     &   D,IRANK)
C
C  SUBROUTINE  A2TAND  TO BE USED WITH ROUTINE ALCON2
C
C  COMPUTATION OF TANGENT AND DETERMINANT AT GIVEN POINT OF
C  CONTINUATION PATH BY QR DECOMPOSITION
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N1),YW(N1),DFY(N,N1),A(N,N1)
      DOUBLE PRECISION AH(N1,N1),V(N1),D(N1)
      INTEGER PIVOT(N1)
      EXTERNAL JFY
      COMMON /COUNT/ IFCTEV,IJACEV,ITER,IDECS,ISOLS
      COMMON /MACHIN/ EPMACH,SMALL
C  COMPUTE NEW TANGENT AND DETERMINANT
      CALL JFY(Y,Y(N1),DFY,DFY(1,N1))
      IJACEV=IJACEV+1
      DO 1010 J=1,N1
      S=YW(J)
      DO 1010 I=1,N
      A(I,J)=DFY(I,J)*S
1010  CONTINUE
      IRANK=N
      COND=1.D0/EPMACH
      CALL A1DECC(A,N,N1,0,N,N1,IRANK,COND,D,PIVOT,0,AH,V)
      IDECS=IDECS+1
      IF (IRANK.LT.N) RETURN
      DET=1.D0
      DO 1020 I=1,N
      DET=DET*D(I)
1020  CONTINUE
      NM1=N-1
      IF (NM1.EQ.0) GOTO 1031
      DO 1030 I=1,NM1
      I1=I+1
      DO 1030 J=I1,N
      IF (PIVOT(J).LT.PIVOT(I)) DET=-DET
1030  CONTINUE
1031  CONTINUE
      RETURN
      END
C
C
      SUBROUTINE A2BIFT (N,N1,INFO,V,W,PIVOT,AH,BS,W1,W2,EPS)
C
C  SUBROUTINE  A2BIFT  TO BE USED WITH ROUTINE ALCON2
C
C  COMPUTATION OF TANGENTS IN BIFURCATION POINT
C
C- CHANGED: LOC-FUNCTION FOR ARRAY BS TO REFLECT CHANGE OF LOC-FUNCTION
C-          IN SUBROUTINE A2NEWA
C- Lutz Weimann  8.1.87
C-
C-
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION V(N1),W(N1)
      DOUBLE PRECISION AH(N1,N1),BS(1)
      DOUBLE PRECISION W1(N1),W2(N1)
      INTEGER INFO(9),PIVOT(N1)
C-    LOC(I,J)=I+(4*N-J)*(J-1)/2
      LOC(I,J)=I*(I-1)/2+J
C  COMPUTATION OF KERNEL OF F'
      I1 = 1
      NM1=N-1
      IF (NM1.EQ.0) GOTO 1031
      DO 1030 J=N,N1
      DO 1020 II=1,NM1
      I=N-II
      S=BS(LOC(N1+I,J))
      IF (II.EQ.1) GOTO 1011
      DO 1010 JJ=I1,NM1
      S=S-BS(LOC(N1+I,JJ))*V(JJ)
1010  CONTINUE
1011  CONTINUE
      I1=I
      V(I)=S/BS(LOC(N1+I,I))
      AH(I,J)=V(I)
1020  CONTINUE
1030  CONTINUE
1031  CONTINUE
C  DETERMINATION OF QUADRATIC EQUATION A*COS**2+B*COS*SIN+C*SIN**2
      A=BS(LOC(N,N))
      IF (NM1.EQ.0) GOTO 2041
      DO 2010 I=1,NM1
      A=A-2.D0*BS(LOC(N,I))*AH(I,N)
2010  CONTINUE
      DO 2030 I=1,NM1
      S=0.D0
      DO 2020 J=1,NM1
      IF (I.GE.J) S=S+BS(LOC(I,J))*AH(J,N)
      IF (I.LT.J) S=S+BS(LOC(J,I))*AH(J,N)
2020  CONTINUE
      W(I)=S
2030  CONTINUE
      DO 2040 I=1,NM1
      A=A+W(I)*AH(I,N)
2040  CONTINUE
2041  CONTINUE
      B=BS(LOC(N1,N))
      IF (NM1.EQ.0) GOTO 2051
      DO 2050 I=1,NM1
      B=B-BS(LOC(N,I))*AH(I,N1)-BS(LOC(N1,I))*AH(I,N)+W(I)*AH(I,N1)
2050  CONTINUE
2051  CONTINUE
      C=BS(LOC(N1,N1))
      IF (NM1.EQ.0) GOTO 2091
      DO 2060 I=1,NM1
      C=C-2.D0*BS(LOC(N1,I))*AH(I,N1)
2060  CONTINUE
      DO 2080 I=1,NM1
      S=0.D0
      DO 2070 J=1,NM1
      IF (I.GE.J) S=S+BS(LOC(I,J))*AH(J,N1)
      IF (I.LT.J) S=S+BS(LOC(J,I))*AH(J,N1)
2070  CONTINUE
      W(I)=S
2080  CONTINUE
      DO 2090 I=1,NM1
      C=C+W(I)*AH(I,N1)
2090  CONTINUE
2091  CONTINUE
C  COMPUTATION OF TANGENTS
      VSIN=0.D0
      VCOS=1.D0
      IF (A.EQ.0.D0 .AND. C.EQ.0.D0) GOTO 3120
      S=B**2-A*C
      IF (S.LT.(B*EPS)**2) GOTO 9900
      S=DSIGN(DSQRT(S),-B)
      S=-B+S
      IF (DABS(A).GT.DABS(C)) GOTO 3110
      R1=S/C
      R2=(A/C)/R1
      VN=DSQRT(1.D0+R1**2)
      VSIN=R1/VN
      VCOS=1.D0/VN
      GOTO 3120
3110  CONTINUE
      R1=S/A
      R2=(C/A)/R1
      VN=DSQRT(1.D0+R1**2)
      VSIN=1.D0/VN
      VCOS=R1/VN
3120  CONTINUE
      VN=1.D0
      IF (NM1.EQ.0) GOTO 3131
      DO 3130 I=1,NM1
      S=-AH(I,N)*VCOS-AH(I,N1)*VSIN
      VN=VN+S**2
      W(PIVOT(I))=S
3130  CONTINUE
3131  CONTINUE
      W(PIVOT(N))=VCOS
      W(PIVOT(N1))=VSIN
      VN=DSQRT(VN)
      DO 3140 I=1,N1
      W1(I)=W(I)/VN
3140  CONTINUE
      VSIN=1.D0
      VCOS=0.D0
      IF (A.EQ.0.D0 .AND. C.EQ.0.D0) GOTO 3220
      VN=DSQRT(1.D0+R2**2)
      IF (DABS(A).GT.DABS(C)) GOTO 3210
      VSIN=R2/VN
      VCOS=1.D0/VN
      GOTO 3220
3210  CONTINUE
      VSIN=1.D0/VN
      VCOS=R2/VN
3220  CONTINUE
      VN=1.D0
      IF (NM1.EQ.0) GOTO 3231
      DO 3230 I=1,NM1
      S=-AH(I,N)*VCOS-AH(I,N1)*VSIN
      VN=VN+S**2
      W(PIVOT(I))=S
3230  CONTINUE
3231  CONTINUE
      W(PIVOT(N))=VCOS
      W(PIVOT(N1))=VSIN
      VN=DSQRT(VN)
      DO 3240 I=1,N1
      W2(I)=W(I)/VN
3240  CONTINUE
C  CHECK ANGLE BETWEEN COMPUTED TANGENTS
      S=0.D0
      DO 4010 I=1,N1
      S=S+W1(I)*W2(I)
4010  CONTINUE
      IF (DABS(S).GT.0.95D0) GOTO 9900
C  SOLUTION EXIT
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  DISCRIMINANTE CONDITION VIOLATED OR DOUBLE TANGENT DETECTED
C  POSSIBLY NON-SIMPLE BIFURCATION POINT DETECTED
9900  CONTINUE
      INFO(9)=1
      RETURN
      END
C
C
      BLOCK DATA A2BLOK
C
C  BLOCK DATA  TO BE USED WITH ROUTINE ALCON2
C
C    COMMON /A2UNIT/ UPR,UDIAG,UBIF
C      INPUT/OUTPUT UNITS
C      - UPR     PRINT UNIT                            OUTPUT
C      - UDIAG   UNIT FOR PLOT INFORMATION             OUTPUT
C                (IN CONNECTION WITH ROUTINE PLTHM)
C      - UBIF    UNIT FOR BIFURCATION INFORMATION      INPUT/OUTPUT
C
C    COMMON /A2VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C      LIMIT VALUES IN CONTINUATION PROCESS
C      - SMALL1  USED AS INITIAL STEPLENGTH ESTIMATE AND AS MAXIMUM
C                PERMITTED STEPSIZE IN NUMERICAL DIFFERENTIATION
C      - THMAX   MAXIMUM PERMITTED RATIO IN RESTRICTED MONOTONICITY TEST
C      - THR     THRESHOLD IN NEARLY LINEAR CASE
C      - ITMAX   MAXIMUM PERMITTED NUMBER OF ITERATIONS
C                IN INTERVAL METHODS
C      - IGNMAX  MAXIMUM PERMITTED NUMBER OF GAUSS-NEWTON ITERATIONS
C                PER CORRECTOR CALL
C
C    COMMON /A2METH/ FACTOR,IEXTR,IJACM
C      METHOD DESCRIBING PARAMETERS
C      - FACTOR  FACTOR FOR EXTRAPOLATED STEPLENGTH BOUND
C      - IEXTR   METHOD FOR EXTRAPOLATED STEPLENGTH BOUND
C                  0  NO EXTRAPOLATED STEPLENGTH BOUND
C                  1  WITH EXTRAPOLATED STEPLENGTH BOUND
C      - IJACM   METHOD FOR UPDATING JACOBIAN IN GAUSS-NEWTON METHOD
C                  0  KEEP FIRST JACOBIAN
C                  1  BROYDEN UPDATES OF JACOBIAN
C                  2  NEW JACOBIAN IN EACH ITERATION
C                  3  AS 1  BUT IN A COMPUTING TIME SAVING WAY
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER UPR,UDIAG
      COMMON /A2UNIT/ UPR,UDIAG
      COMMON /A2VALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
      COMMON /A2METH/ FACTOR,IEXTR,IJACM
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      DATA UPR/6/,UDIAG/2/
      DATA SMALL1/0.01D0/,THMAX/0.25D0/,THR/0.1D0/,ITMAX/20/,IGNMAX/10/
      DATA FACTOR/1.1D0/,IEXTR/1/,IJACM/3/
      DATA IP/1/,IPI/1/,IPD/1/,IPMAX/100/,IIMAX/1000/,IDMAX/10000/
      END
Crr
      SUBROUTINE A2RRD(N,X,Y,Z)
      INTEGER N
      DOUBLE PRECISION X(N),Y(N),Z(N)
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      INTEGER K
C
      IF (ICHECK(IP).NE.1) THEN
        WRITE(6,*) 'Calling missmatch (1) ICHECK(IP)=',ICHECK(IP)
        STOP
      ENDIF
      IP = IP+1
      IF (ISTACK(IPI).NE.N) THEN
        WRITE(6,*) 'Size missmatch N=',N,' should be ',ISTACK(IPI)
        STOP
      ENDIF
      IPI = IPI+1
      DO 100 K=1,N 
        X(K) = DSTACK(IPD)
        IPD = IPD+1
        Y(K) = DSTACK(IPD)
        IPD = IPD+1
        Z(K) = DSTACK(IPD)
        IPD = IPD+1
100   CONTINUE
      RETURN
      END
Crr
      SUBROUTINE A2RRD4(N,X1,X2,X3)
      INTEGER N
      DOUBLE PRECISION X1,X2,X3
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      IF (ICHECK(IP).NE.2) THEN
        WRITE(6,*) 'Calling missmatch (2) ICHECK(IP)=',ICHECK(IP)
        STOP
      ENDIF
      IP = IP+1
      N = ISTACK(IPI)
      IPI = IPI+1
      X1 = DSTACK(IPD)
      IPD = IPD+1
      X2 = DSTACK(IPD)
      IPD = IPD+1
      X3 = DSTACK(IPD)
      IPD = IPD+1
      RETURN
      END
Crr
      SUBROUTINE A2RRD3(N,X1,X2)
      INTEGER N
      DOUBLE PRECISION X1,X2
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      IF (ICHECK(IP).NE.3) THEN
        WRITE(6,*) 'Calling missmatch (3) ICHECK(IP)=',ICHECK(IP)
        STOP
      ENDIF
      IP = IP+1
      N = ISTACK(IPI)
      IPI = IPI+1
      X1 = DSTACK(IPD)
      IPD = IPD+1
      X2 = DSTACK(IPD)
      IPD = IPD+1
      RETURN
      END
Crr
      SUBROUTINE A2RWR(N,X,Y,Z)
      INTEGER N
      DOUBLE PRECISION X(N),Y(N),Z(N)
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      INTEGER K
C
      IF (IP.GE.IPMAX) THEN
        WRITE(6,*) 'Stack overflow IP=',IP
        STOP
      ENDIF
      ICHECK(IP) = 1
      IP = IP+1
      IF (IPI.GE.IIMAX) THEN
        WRITE(6,*) 'Stack overflow IPI=',IPI
        STOP
      ENDIF
      ISTACK(IPI) = N
      IPI = IPI+1
      IF (IPD.GE.(IDMAX-3*N)) THEN
        WRITE(6,*) 'Stack overflow IPD=',IPD
        STOP
      ENDIF
      DO 100 K=1,N 
        DSTACK(IPD) = X(K)
        IPD = IPD+1
        DSTACK(IPD) = Y(K)
        IPD = IPD+1
        DSTACK(IPD) = Z(K)
        IPD = IPD+1
100   CONTINUE
      RETURN
      END
Crr
      SUBROUTINE A2RWR4(N,X1,X2,X3)
      INTEGER N
      DOUBLE PRECISION X1,X2,X3
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      IF (IP.GE.IPMAX) THEN
        WRITE(6,*) 'Stack overflow IP=',IP
        STOP
      ENDIF
      ICHECK(IP) = 2
      IP = IP+1
      IF (IPI.GE.IIMAX) THEN
        WRITE(6,*) 'Stack overflow IPI=',IPI
        STOP
      ENDIF
      ISTACK(IPI) = N
      IPI = IPI+1
      IF (IPD.GE.(IDMAX-3)) THEN
        WRITE(6,*) 'Stack overflow IPD=',IPD
        STOP
      ENDIF
      DSTACK(IPD) = X1
      IPD = IPD+1
      DSTACK(IPD) = X2
      IPD = IPD+1
      DSTACK(IPD) = X3
      IPD = IPD+1
      RETURN
      END
Crr
      SUBROUTINE A2RWR3(N,X1,X2)
      INTEGER N
      DOUBLE PRECISION X1,X2
C
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
C
      IF (IP.GE.IPMAX) THEN
        WRITE(6,*) 'Stack overflow IP=',IP
        STOP
      ENDIF
      ICHECK(IP) = 3
      IP = IP+1
      IF (IPI.GE.IIMAX) THEN
        WRITE(6,*) 'Stack overflow IPI=',IPI
        STOP
      ENDIF
      ISTACK(IPI) = N
      IPI = IPI+1
      IF (IPD.GE.(IDMAX-2)) THEN
        WRITE(6,*) 'Stack overflow IPD=',IPD
        STOP
      ENDIF
      DSTACK(IPD) = X1
      IPD = IPD+1
      DSTACK(IPD) = X2
      IPD = IPD+1
      RETURN
      END
Crr
      SUBROUTINE A2RREW
      INTEGER IP,IPI,IPD,IPMAX,IDMAX,ISTACK,ICHECK
      DOUBLE PRECISION DSTACK
      COMMON /IO/DSTACK(10000),IP,IPI,IPD,IPMAX,IIMAX,IDMAX,
     *           ICHECK(100),ISTACK(1000)
      IP = 1
      IPI = 1
      IPD = 1
      RETURN
      END
