C ---------------------------------------------------------- 
C     NUMERICAL SOLUTION OF A STIFF DIFFERENTIAL  
C     (OR DIFFERENTIAL ALGEBRAIC) SYSTEM OF FIRST 0RDER  
C     DELAY DIFFERENTIAL EQUATIONS   
C                     M*Y'(X)=F(X,Y(X),Y(X-A),...). 
C     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I) 
C     OR EXPLICIT (M = I). 
C     NOTE: THIS FORM ALSO ALLOWS TO SOLVE NEUTRAL DIFFERENTIAL PROBLEMS 
C 
C     NOTE: THIS VERSION ALLOWS ARBITRARILY LARGE STEP SIZES 
C     (POSSIBLY LARGER THAN THE DELAYS) 
C 
C     THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (3 STAGE  
C     RADAU IIA) OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS  
C     EXTENSION OF ORDER 3  (C.F. SECTION IV.8 OF (HW)) 
C 
C     AUTHORS: N. GUGLIELMI(*) AND E. HAIRER($)  
C          (*) DIVISION OF MATHEMATICS, GRAN SASSO SCIENCE INSTITUTE
C              VIA CRISPI 7, 67100 L'AQUILA, ITALY 
C          ($) UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES 
C              CH-1211 GENEVE 24, SWITZERLAND  
C              E-MAIL ADRESSES:   
C                        nicola.guglielmi@gssi.it 
C                            ernst.hairer@unige.ch 
C      
C     THIS PROGRAM EXTENDS THE CODE RADAU5 (BY E. HAIRER AND G. WANNER) 
C     TO THE CASE OF DELAY DIFFERENTIAL EQUATIONS.  
C     DETAILS ABOUT RADAU5 CAN BE FOUND IN THE BOOK: 
C     (HW)  E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL 
C           EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. 
C           SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, 
C           SPRINGER-VERLAG 1991, SECOND EDITION 1996. 
C     DETAILS ABOUT RADAR5 CAN BE FOUND IN THE PAPERS:
C     (GH1) N. GUGLIELMI AND E. HAIRER, Implementing Radau IIA methods for
C           stiff delay differential equations , Computing 67, 1-12 (2001). 
C     (GH2) N. GUGLIELMI AND E. HAIRER, Computing breaking points in
C           implicit delay differential equations, 
C           Adv. Comput. Math. 29 (2008) 229-247    
C
C ===============================================================================================      
C     VERSION 2.2 OF FEBRUARY 8, 2024 
C     FIXES A FEW BUGS, AND MAKES OTHER LITTLE CHANGES WITH RESPECT TO PREVIOUS VERSION 2.1, 
C     ALL CHANGES HAVE MINIMAL IMPACT, IF ANY, ON BACKWARD COMPATIBILITY WITH VERSION 2.1. 
C     BASIC DIFFERENCES FROM VERSION 2.1: CONCERN THE CALL TO RADAR5, WHICH NOW INCLUDES A FEW
C     FURTHER PARAMETERS CONCERNING VECTOR DIMENSIONS.
C     THE SUBROUTINE CONTR5 - WHICH WAS GIVEN AS AN EXTERNAL FILES - IS NOW INCLUDED IN THIS FILE.    
C     THE ALGORITHM IMPLEMENTING THE NUMERICAL INTEGRATOR HAS NOT CHANGED WITH RESPECT TO VERSION 
C     2.1. 
C ===============================================================================================
C      
C     AN ADDITIONAL SEPARATED ROUTINE (DC_SUMEXPDEL) MIGHT BE USED TO TREAT DISTRIUTED DELAYS BY 
C     A FAST LINEAR SOLVER. IT SHOULD REPLACE THE SUBROUTINE DC_DECDEL
C     A NEW DRIVER PROGRAM FOR AN EXAMPLE WITH IS AVAILABLE IN THE SAME WEBPAGE OF THE CODE. 
C      
C ---------------------------------------------------------- 
      SUBROUTINE RADAR5(N,FCN,PHI,ARGLAG,X,Y,XEND,H, 
     &                  RTOL,ATOL,ITOL, 
     &                  JAC,IJAC,MLJAC,MUJAC, 
     &                  JACLAG,NLAGS,NJACL, 
     &                  MAS,IMAS,MLMAS,MUMAS,SOLOUT,IOUT, 
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID, 
     &                  GRID,LGRID,IPAST,NRDENS) 
C ---------------------------------------------------------- 
C     INPUT PARAMETERS   
C --------------------   
C     N           DIMENSION OF THE SYSTEM  
C 
C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE RIGHT- 
C                 HAND-SIDE OF THE DELAY EQUATION, E.G., 
C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR,...) 
C                    DOUBLE PRECISION X,Y(N),F(N) 
C                    EXTERNAL PHI 
C                    F(1)=G1(X,Y(*),YLAGR5(*,X-TAU(X,Y(*))),PHI,...)) 
C                    F(2)=G2(X,Y(*),YLAGR5(*,X-TAU(X,Y(*))),PHI,...)) 
C                    ETC. 
C                    (*) MEANS ALL POSSIBLE COMPONENTS 
C                 FOR AN EXPLICATION OF YLAGR5 SEE BELOW. 
C                 DO NOT USE YLAGR5(I,X-0.D0,PHI,RPAR,IPAR,...) ! 
C                 Note: 
C                 THE INITIAL FUNCTION HAS TO BE SUPPLIED BY: 
C                    FUNCTION PHI(I,X,RPAR,IPAR) 
C                    DOUBLE PRECISION PHI,X 
C                 WHERE I IS THE COMPONENT AND X THE ARGUMENT 
C                 RPAR, IPAR (SEE BELOW) 
C 
C     PHI         NAME (EXTERNAL) OF FUNCTION COMPUTING THE INITIAL
C                 FUNCTIONS FOR COMONENTS I (WITH RETARDED ARGUMENT),
C                    FUNCTION PHI(I,X,RPAR,IPAR)
C                    DOUBLE PRECISION PHI,X
C
C     ARGLAG      NAME (EXTERNAL) OF FUNCTION COMPUTING THE IL-TH LAG TERM,
C                    FUNCTION ARGLAG(IL,X,N,Y,RPAR,IPAR,PHI,...)
C                    DOUBLE PRECISION PHI,X,Y(N)
C     X           INITIAL X-VALUE 
C 
C     Y(N)        INITIAL VALUES FOR Y (MAY BE DIFFERENT FROM PHI (I,X), 
C                 IN THIS CASE IT IS HIGHLY RECOMMENDED TO SET IWORK(13) 
C                 AND GRID(1),..., (SEE BELOW) 
C 
C     XEND        FINAL X-VALUE (XEND-X HAS TO BE POSITIVE) 
C 
C     H           INITIAL STEP SIZE GUESS; 
C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,  
C                 H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD. 
C                 THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS 
C                 QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6). 
C 
C     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY 
C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. 
C 
C     ITOL        SWITCH FOR RTOL AND ATOL: 
C                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. 
C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF 
C                     Y(I) OVER RTOL*ABS(Y(I))+ATOL 
C                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. 
C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) OVER 
C                     RTOL(I)*ABS(Y(I))+ATOL(I). 
C 
C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES 
C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y 
C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1;  
C                 THE USER HAS TO SUPPLY A DUMMY SUBROUTINE  
C                 IN THE CASE IJAC=0). 
C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM 
C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR,...) 
C                    DOUBLE PRECISION X,Y(N),DFY(LDFY,N) 
C                    DFY(1,1)= ... 
C                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS 
C                 FURNISHED BY THE CALLING PROGRAM. 
C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO 
C                    BE FULL AND THE PARTIAL DERIVATIVES ARE 
C                    STORED IN DFY AS 
C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) 
C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND 
C                    THE PARTIAL DERIVATIVES ARE STORED 
C                    DIAGONAL-WISE AS 
C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). 
C 
C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN: 
C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE 
C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. 
C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC. 
C 
C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: 
C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR 
C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. 
C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN  
C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW 
C                       THE MAIN DIAGONAL). 
C 
C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON- 
C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). 
C                 DOES NOT NEED TO BE DEFINED IF MLJAC=N. 
C 
C     JACLAG      NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES 
C                 THE PARTIAL DERIVATIVES OF F(X,Y,YLAG) WITH RESPECT TO  
C                 YLAG(*) (YLAG DENOTE THE DELAYED VARIABLES) 
C 
C     NLAGS       DENOTES THE NUMBER OF DELAY ARGUMENTS.  
C                 THIS PARAMETER IS OF INTEREST FOR THE COMPUTATION OF THE 
C                 JACOBIAN. 
C                 TO BE SET = 0 IF ONE DOES WANT TO COMPUTE THE TRADITIONAL 
C                 JACOBIAN;  
C                 TO BE SET = NUMBER OF DISTINCT DELAY ARGUMENTS 
C                 IF ONE WANTS TO CORRECT THE STANDARD JACOBIAN (THROUGH 
C                 THE SUBROUTINE JACLAG) WHEN ADVANCED ARGUMENTS ARE USED. 
C 
C     NJACL       NUMBER OF TERMS IN THE JACOBIAN W.R.T. 
C                 RETARDED COMPONENTS (WHICH IS THOUGHT AS A SPARSE MATRIX). 
C 
C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      ----- 
C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - 
C 
C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- 
C                 MATRIX M. 
C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY 
C                 MATRIX AND NEEDS NOT TO BE DEFINED; 
C                 THE USER HAS TO SUPPLY A DUMMY SUBROUTINE IN THIS CASE. 
C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM 
C                    SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR) 
C                    DOUBLE PRECISION AM(LMAS,N) 
C                    AM(1,1)= .... 
C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED 
C                    AS FULL MATRIX LIKE 
C                         AM(I,J) = M(I,J) 
C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED 
C                    DIAGONAL-WISE AS 
C                         AM(I-J+MUMAS+1,J) = M(I,J). 
C 
C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX: 
C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY 
C                       MATRIX, MAS IS NEVER CALLED. 
C                    IMAS=1: MASS-MATRIX  IS SUPPLIED. 
C 
C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: 
C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR 
C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. 
C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE 
C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW 
C                       THE MAIN DIAGONAL). 
C                 MLMAS IS SUPPOSED TO BE <= MLJAC. 
C 
C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- 
C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). 
C                 DOES NOT NEED TO BE DEFINED IF MLMAS=N. 
C                 MUMAS IS SUPPOSED TO BE <= MUJAC. 
C 
C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE 
C                 NUMERICAL SOLUTION DURING INTEGRATION.  
C                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. 
C                 THE USER HAS TO SUPPLY A DUMMY SUBROUTINE IF IOUT=0.  
C                 IT MUST HAVE THE FORM 
C                    SUBROUTINE SOLOUT (NR,XOLD,X,HSOL,Y,CONT,LRC,N, 
C                                       RPAR,IPAR,IRTRN) 
C                    DOUBLE PRECISION X,Y(N),CONT(LRC) 
C                    ....   
C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH 
C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS 
C                    THE FIRST GRID-POINT). 
C                 "XOLD" IS THE PRECEEDING GRID-POINT. 
C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN 
C                    IS SET <0, RADAR5 RETURNS TO THE CALLING PROGRAM. 
C            
C          -----  CONTINUOUS OUTPUT: ----- 
C                 DURING CALLS TO "SOLOUT" AS WELL AS TO "FCN", A 
C                 CONTINUOUS SOLUTION IS AVAILABLE THROUGH HTHE FUNCTION 
C                        >>>   YLAGR5(I,S,PHI,RPAR,IPAR,...)   <<< 
C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH 
C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE S 
C                 HAS TO LIE IN AN INTERVAL WHERE THE NUMERICAL SOLUTION 
C                 IS ALREADY COMPUTED. IT DEPENDS ON THE SIZE OF LRPAST 
C                 (SEE BELOW) HOW FAR BACK THE SOLUTION IS AVAILABLE. 
C 
C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT: 
C                    IOUT=0: SUBROUTINE IS NEVER CALLED 
C                    IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT. 
C 
C     WORK        ARRAY OF STATE VARIABLES OF REAL TYPE FOR EXECUTION. 
C                 WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS 
C                 FOR THE CODE. FOR STANDARD USE OF THE CODE 
C                 WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE 
C 
C     LWORK       DECLARED LENGTH OF ARRAY WORK
C
C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK". 
C                 IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS 
C                 FOR THE CODE. FOR STANDARD USE, SET IWORK(1),.., 
C                 IWORK(20) TO ZERO BEFORE CALLING. 
C 
C     LIWORK       DECLARED LENGTH OF ARRAY IWORK
C 
C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH   
C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING 
C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES.  
C 
C     GRID        CONTAINS PRESCRIBED GRID POINTS, WHICH THE 
C                 INTEGRATION METHOD HAS TO TAKE AS GRID-POINTS 
C                 NORMALLY, IF GRID(1) > X, THEN ONE HAS 
C                 X < GRID(1) < GRID(2) < ... < GRID(NGRID) <= XEND 
C                 IN SOME CASES IF THERE ARE DISCONTINUITIES IN THE  
C                 INITIAL FUNCTIONS THEY ARE ALSO SET IN THE GRID  
C                 VECTOR; THEN X < GRID(J) < GRID(J+1) ... < XEND 
C                 WHERE J ADDRESSES THE FIRST ASCISSA IN GRID > X 
C 
C     LGRID       DECLARED LENGTH OF GRID VECTOR, 
C                 GRID(LGRID), 
C                 WHICH MUST BE DECLARED IN THE CALLING PROGRAM. 
C                 "LGRID" MUST BE AT LEAST 
C                              NGRID + 1 
C 
C     IPAST       FOR   0 < NRDENS < N   THE COMPONENTS (FOR WHICH DENSE 
C                 OUTPUT IS REQUIRED) HAVE TO BE SPECIFIED IN 
C                 IPAST(1),...,IPAST(NRDENS),
C                 FOR  NRDENS=N  THIS IS DONE BY THE CODE. 
C
C     NRDENS =    NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT 
C                 IS REQUIRED (EITHER BY "SOLOUT" OR BY "FCN"); 
C ---------------------------------------------------------------------- 
C  
C     SOPHISTICATED SETTING OF PARAMETERS 
C     ----------------------------------- 
C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK  
C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),... 
C              AS WELL AS IWORK(1),... DIFFERENT FROM ZERO. 
C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: 
C 
C    IWORK(1)  IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN 
C              MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY 
C              ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN. 
C              IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N) 
C              AND NOT FOR IMPLICIT SYSTEMS (IMAS=1). 
C 
C    IWORK(2)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. 
C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000. 
C 
C    IWORK(3)  THE MAXIMUM NUMBER OF NEWTON ITERATIONS FOR THE 
C              SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP. 
C              THE DEFAULT VALUE (FOR IWORK(3)=0) IS 7. 
C 
C    IWORK(4)  IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION 
C              IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD. 
C              IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED. 
C              THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS 
C              DIFFICULTIES WITH CONVERGENCE (THIS IS SEEN IN THE CASE WHEN 
C              NSTEP IS LARGER THAN NACCPT + NREJCT; SEE OUTPUT PARAM.). 
C              DEFAULT IS IWORK(4)=0. 
C 
C       THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR 
C       DELAY DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1. 
C       THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT 
C       THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER.  
C       IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE 
C       MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2. 
C 
C    IWORK(5)  DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR  
C              DDE'S THIS EQUALS THE DIMENSION OF THE SYSTEM. 
C              DEFAULT IWORK(5)=N. 
C 
C    IWORK(6)  DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0. 
C 
C    IWORK(7)  DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0. 
C 
C    IWORK(8)  SWITCH FOR STEP SIZE STRATEGY 
C              IF IWORK(8).EQ.1  MODIFIED PREDICTIVE CONTROLLER  
C              (GUSTAFSSON) 
C              IF IWORK(8).EQ.2  CLASSICAL STEP SIZE CONTROL 
C              THE DEFAULT VALUE (FOR IWORK(8)=0) IS IWORK(8)=1. 
C              THE CHOICE IWORK(8).EQ.1 SEEMS TO PRODUCE SAFER RESULTS; 
C              FOR SIMPLE PROBLEMS, THE CHOICE IWORK(8).EQ.2 PRODUCES 
C              OFTEN SLIGHTLY FASTER RUNS 
C 
C       IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT 
C            Y(I)' = Y(I+M2)   FOR  I=1,...,M1, 
C       WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME 
C       CAN BE ACHIEVED BY SETTING THE PARAMETERS IWORK(9) AND IWORK(10). 
C       E.G., FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE  
C       VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2. 
C       FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS: 
C       - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE 
C              JACOBIAN HAVE TO BE STORED 
C              IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL 
C                 DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J) 
C                FOR I=1,N-M1 AND J=1,N. 
C              ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM ) 
C                 DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2) 
C                FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM. 
C       - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL 
C                0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM) 
C                     PARTIAL F(I+M1) / PARTIAL Y(J+K*M2),  I,J=1,M2 
C                    ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH 
C                    OF THESE MM+1 SUBMATRICES 
C       - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES; 
C                DOES NOT NEED TO BE DEFINED IF MLJAC=N-M1 
C       - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND 
C              DOES NOT NEED TO BE DEFINED.  
C              THE USER HAS TO SUPPLY A DUMMY SUBROUTINE IN THIS CASE. 
C              IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF 
C              DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX. 
C              IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL 
C                 AM(I,J) = M(I+M1,J+M1)     FOR I=1,N-M1 AND J=1,N-M1. 
C              ELSE, THE MASS MATRIX IS BANDED 
C                 AM(I-J+MUMAS+1,J) = M(I+M1,J+M1) 
C       - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL 
C                0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX 
C       - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX 
C                DOES NOT NEED TO BE DEFINED IF MLMAS=N-M1 
C 
C    IWORK(9)  THE VALUE OF M1.  DEFAULT M1=0. 
C 
C    IWORK(10) THE VALUE OF M2.  DEFAULT M2=M1. 
C 
C 
C 
C    IWORK(11) SELECT THE TYPE OF ERROR CONTROL:  
C              -1: FOR A PURE CONTROL OF THE DENSE OUTPUT  
C                  (MAKES USE OF A QUADRATIC AND A LINEAR INTERPOLATING 
C                   POLYNOMIALS);  
C               1: FOR A MIXED CONTROL OF DENSE OUTPUT AND DISCRETE OUTPUT 
C               2: FOR A PURE CONTROL OF THE DISCRETE OUTPUT 
C                  (ERROR CONTROL PROVIDED BY THE SUBROUTINE ESTRAD). 
C               3: FOR A SIMPLER MIXED CONTROL OF DENSE OUTPUT AND   
C                  DISCRETE OUTPUT 
C               DEFAULT VALUE IWORK(11)=2. 
C 
C    IWORK(12) = MXST = (ON INPUT)  
C              DECLARED NUMBER OF STEPS STORED IN THE ``PAST VECTOR'',  
C              PAST(LRPAST), 
C              WHICH MUST BE DECLARED IN THE CALLING PROGRAM. 
C              "MXST" MUST BE SUFFICIENTLY LARGE. IF THE DENSE 
C              OUTPUT OF MXST BACK STEPS HAS TO BE STORED,  
C              THE DIMENSION OF PAST MUST BE  
C                       LRPAST=MXST*(4*NRDENS+2) 
C              WHERE NRDENS=IWORK(15) (SEE BELOW). 
C 
C    IWORK(13) = NGRID = (ON INPUT) 
C              NUMBER OF PRESCRIBED POINTS IN THE 
C              INTEGRATION INTERVAL WHICH HAVE TO BE GRID-POINTS 
C              IN THE INTEGRATION. USUALLY, AT THESE POINTS THE 
C              SOLUTION OR ONE OF ITS DERIVATIVE HAS A DISCONTINUITY. 
C              DEFINE THESE POINTS IN GRID(1),...,GRID(NGRID) 
C              DEFAULT VALUE:  IWORK(13)=0 
C 
C    IWORK(14) = SELECTOR FOR FULL ITERATION (2) OR SIMPLIFIED  
C              ITERATION (1) (TAKING INTO ACCOUNT POSSIBLE  
C              ADVANCED ARGUMENTS BUT PRESERVING TENSOR STRUCTURE  
C              OF THE JACOBIAN. 
C              DEFAULT VALUE:  IWORK(14)=1 
C 
C    IWORK(16) = NDIMN = (ON INPUT)
C              OPTION VALID FOR NEUTRAL PROBLEMS
C              NUMBER OF DERIVATIVE COMPONENTS (Z) OF THE NEUTRAL PROBLEM 
C              EXCLUDED TRUE SOLUTION COMPONENTS
C 
C ---------- 
C 
C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. 
C 
C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION, 
C              DEFAULT 0.9D0. 
C 
C    WORK(3)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; 
C              INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS 
C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER  
C              (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE TO 
C              COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP.      
C              DEFAULT 0.001D0. 
C 
C    WORK(4)   STOPPING CRITERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1. 
C              SMALLER VALUES OF WORK(4) MAKE THE CODE SLOWER, BUT SAFER. 
C              DEFAULT MIN(0.03D0,RTOL(1)**0.5D0) 
C 
C    WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE 
C              STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A 
C              LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR 
C              LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE 
C              WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS 
C              WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD. 
C              DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 . 
C 
C    WORK(7)   MAXIMAL STEP SIZE, DEFAULT XEND-X. 
C 
C    WORK(8), WORK(9)   PARAMETERS FOR STEP SIZE SELECTION 
C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION 
C                 WORK(8) <= HNEW/HOLD <= WORK(9) 
C              DEFAULT VALUES: WORK(8)=0.2D0, WORK(9)=8.D0 
C 
C    WORK(10)  PARAMETER FOR CONTROLLING THE ERROR CONTROL OF DENSE 
C              OUTPUT (0 <= WORK(10) <= 1). (0: STRONG CONTROL, 1: WEAKER) 
C              SUGGESTED VALUES: 
C              FOR PROBLEMS WITH `ALMOST DISCONTINUOUS' SOLUTIONS  
C              (LIKE SHOCKS):  WORK(10)=0.D0 
C              FOR PROBLEMS WITH FAIRLY SMOOTH SOLUTION:  WORK(10)=1.D0 
C              FOR INTERMEDIATE PROBLEMS:  WORK(10)=1.D-M (M=1,2,3,..,) 
C              DEFAULT VALUE: WORK(10)=0.D0  
C    WORK(11)  PARAMETER FOR CONTROLLING THE SEARCH OF BREAKING POINTS:
C              IF THE ERROR INCREASES OF A FACTOR LARGER THAN WORK(11)
C              FROM A STEP TO THE SUBSEQUENT, THE ROUTINE SEARCHING BREAKING
C              POINTS ACTIVATES.
C              DEFAULT VALUE: WORK(11)=5.D0
C----------------------------------------------------------------------- 
C 
C     OUTPUT PARAMETERS  
C     -----------------  
C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED 
C                 (AFTER SUCCESSFUL RETURN X=XEND). 
C 
C     Y(N)        NUMERICAL SOLUTION AT X 
C  
C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP 
C 
C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN: 
C                   IDID= 1  COMPUTATION SUCCESSFUL, 
C                   IDID= 2  COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT) 
C                   IDID=-1  INPUT IS NOT CONSISTENT, 
C                   IDID=-2  LARGER NMAX IS NEEDED, 
C                   IDID=-3  STEP SIZE BECOMES TOO SMALL, 
C                   IDID=-4  MATRIX IS REPEATEDLY SINGULAR. 
C                   IDID=-5  COMPUTATION INTERRUPTED BY YLAGR5.    
C                   IDID=-6  THE EQUATION USES ADVANCED ARGUMENTS 
C 
C   IWORK(13)  NFULL   NUMBER OF FULL NEWTON ITERATIONS                    
C   IWORK(14)  NFCN    NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL 
C                      EVALUATION OF THE JACOBIAN ARE NOT COUNTED)   
C   IWORK(15)  NJAC    NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY 
C                      OR NUMERICALLY) 
C   IWORK(16)  NSTEP   NUMBER OF COMPUTED STEPS 
C   IWORK(17)  NACCPT  NUMBER OF ACCEPTED STEPS 
C   IWORK(18)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST), 
C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED) 
C   IWORK(19)  NDEC    NUMBER OF LU-DECOMPOSITIONS OF BOTH MATRICES 
C   IWORK(20)  NSOL    NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF BOTH 
C                      SYSTEMS; THE NSTEP FORWARD-BACKWARD SUBSTITUTIONS, 
C                      NEEDED FOR STEP SIZE SELECTION, ARE NOT COUNTED 
C----------------------------------------------------------------------- 
C *** *** *** *** *** *** *** *** *** *** *** *** *** 
C          DECLARATIONS  
C *** *** *** *** *** *** *** *** *** *** *** *** *** 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      INTEGER, PARAMETER :: DP=kind(1D0) 
      REAL(kind=DP), dimension(N), intent(inout) :: Y 
      REAL(kind=DP), dimension(LWORK), intent(inout) :: WORK 
      REAL(kind=DP), dimension(1), intent(inout) :: ATOL,RTOL 
      INTEGER, dimension(LIWORK), intent(inout) :: IWORK 
      REAL(kind=DP), dimension(LGRID), intent(inout) :: GRID 
      INTEGER, dimension(NRDENS+2*N), intent(inout) :: IPAST 
      REAL(kind=DP), dimension(1), intent(in) :: RPAR 
      INTEGER, dimension(1), intent(in) :: IPAR 
      LOGICAL IMPLCT,NEUTRAL,JBAND,ARRET,STARTN,PRED 
      LOGICAL FLAGS,FLAGN 
      EXTERNAL FCN,PHI,ARGLAG,JAC,JACLAG,MAS,SOLOUT 
C ----> COMMON BLOCKS <---- 
      COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
C *** *** *** *** *** *** *** 
C        SETTING THE PARAMETERS  
C *** *** *** *** *** *** *** 
       NN=N 
       NFCN=0 
       NJAC=0 
       NSTEP=0 
       NACCPT=0 
       NREJCT=0 
       NDEC=0 
       NSOL=0 
       ARRET=.FALSE. 
       FLAGS=.FALSE. 
       FLAGN=.FALSE. 
 
       IF (IOUT.EQ.1) WRITE (6,*) 'STARTING INTEGRATION...' 
C        
C ------> OPERATIONS RELEVANT TO THE DELAY DEPENDENCE <------ 
C 
C -------- ERROR CONTROL 
      IF (IWORK(11).EQ.0) THEN 
       IEFLAG=2 
      ELSE 
       IEFLAG=IWORK(11) 
      END IF 
      IF (IEFLAG.EQ.2) WORK(10)=1.D0 
C -------- NGRID   NUMBER OF PRESCRIBED GRID-POINTS 
      NGRID=IWORK(13) 
      IF (NGRID.LT.0) NGRID=0 
      IF (IOUT.EQ.1) WRITE(6,*)  
     &           'NUMBER OF PRESCRIBED GRID POINTS: ',NGRID 
C ------- NDIMN   NUMBER OF COMPONENTS OF A NEUTRAL PROBLEM
      IF (IMAS.EQ.2) THEN
	    IF (IWORK(16).EQ.0) THEN
	      WRITE(6,*) 'NUMBER OF Y COMPONENTS HAS TO BE SPECIFIED'
	      ARRET=.TRUE.
	    END IF
	    NDIMN=IWORK(16)
	  ELSE 
        NDIMN=0
	  END IF  	 
C ------- LIPAST   DIMENSION OF VECTOR IPAST 
      LIPAST=NRDENS+1  
      IF(NRDENS.LT.0.OR.NRDENS.GT.N) THEN 
         IF (IOUT.GT.0) WRITE(6,*) 
     &           ' CURIOUS INPUT IWORK(15)=',IWORK(15) 
         ARRET=.TRUE. 
      ELSE 
         NRDS=NRDENS 
      END IF 
      IF (NRDS.EQ.N) THEN 
            DO I=1,NRDS 
              IPAST(I)=I 
            END DO
      END IF 
      IF (IOUT.EQ.1) WRITE(6,*) 'NUMBER OF DELAYED COMPONENTS: ',NRDS 
C ------- LRPAST   DIMENSION OF VECTOR PAST 
      MXST=IWORK(12) 
C ------- CONTROL OF LENGTH OF PAST  ------- 
      IF(MXST.LT.1)THEN 
         IF (IOUT.GT.0) WRITE(6,*) 
     & ' INSUFFICIENT STORAGE FOR PAST, MIN. LRPAST=',1 
         ARRET=.TRUE. 
      END IF 
C ------- DIM. of PAST  -------- 
      IDIF=4*NRDS+2 
      LRPAST=MXST*IDIF 
C -------------------------------------------------  
C ------- CONTROL OF SIMPLE NEWTON ITERATION  ------- 
      ISWJL=IWORK(14) 
      IF (ISWJL.EQ.0) ISWJL=1 
 
C -------- UROUND : SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0   
      IF (WORK(1).EQ.0.0D0) THEN 
         UROUND=1.0D-16 
      ELSE 
         UROUND=WORK(1) 
         IF (UROUND.LE.1.0D-19.OR.UROUND.GE.1.0D0) THEN 
            WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1) 
            ARRET=.TRUE. 
         END IF 
      END IF 
 
C -------> CHECK AND CHANGE THE TOLERANCES <------ 
      EXPM=2.0D0/3.0D0 
      IF (ITOL.EQ.0) THEN 
          IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN 
              WRITE (6,*) ' TOLERANCES ARE TOO SMALL' 
              ARRET=.TRUE. 
          ELSE 
              QUOT=ATOL(1)/RTOL(1) 
              RTOL(1)=0.1D0*RTOL(1)**EXPM 
              ATOL(1)=RTOL(1)*QUOT 
          END IF 
      ELSE 
          DO I=1,N 
          IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN 
              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' 
              ARRET=.TRUE. 
          ELSE 
              QUOT=ATOL(I)/RTOL(I) 
              RTOL(I)=0.1D0*RTOL(I)**EXPM 
              ATOL(I)=RTOL(I)*QUOT 
          END IF 
          END DO 
      END IF 
 
C -------> NMAX : THE MAXIMAL NUMBER OF STEPS <------- 
      IF (IWORK(2).EQ.0) THEN 
         NMAX=100000 
      ELSE 
         NMAX=IWORK(2) 
         IF (NMAX.LE.0) THEN 
            WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2) 
            ARRET=.TRUE. 
         END IF 
      END IF 
 
C -------> NIT :  MAXIMAL NUMBER OF NEWTON ITERATIONS <------- 
      IF (IWORK(3).EQ.0) THEN 
         NIT=7 
      ELSE 
         NIT=IWORK(3) 
         IF (NIT.LE.0) THEN 
            WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3) 
            ARRET=.TRUE. 
         END IF 
      END IF 
C -------- STARTN : SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS 
      IF(IWORK(4).EQ.0)THEN 
         STARTN=.FALSE. 
      ELSE 
         STARTN=.TRUE. 
      END IF 
 
C -------> PARAMETERS (IF ANY) FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS <------- 
      NIND1=IWORK(5) 
      NIND2=IWORK(6) 
      NIND3=IWORK(7) 
      IF (NIND1.EQ.0) NIND1=N 
      IF (NIND1+NIND2+NIND3.NE.N) THEN 
       WRITE(6,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',NIND1,NIND2,NIND3 
       ARRET=.TRUE. 
      END IF 
 
C -------> PRED   STEP SIZE CONTROL <------- 
      IF(IWORK(8).LE.1)THEN 
         PRED=.TRUE. 
      ELSE 
         PRED=.FALSE. 
      END IF 
 
C -------> PARAMETER FOR SECOND ORDER EQUATIONS <------- 
      M1=IWORK(9) 
      M2=IWORK(10) 
      NM1=N-M1 
      IF (M1.EQ.0) M2=N 
      IF (M2.EQ.0) M2=M1 
      IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN 
       WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2 
       ARRET=.TRUE. 
      END IF 
 
C -------> SAFE  :  SAFETY FACTOR IN STEP SIZE PREDICTION <------- 
      IF (WORK(2).EQ.0.0D0) THEN 
         SAFE=0.9D0 
      ELSE 
         SAFE=WORK(2) 
         IF (SAFE.LE.0.001D0.OR.SAFE.GE.1.0D0) THEN 
            WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2) 
            ARRET=.TRUE. 
         END IF 
      END IF 
 
C ------> THET : DETERMINES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; 
      IF (WORK(3).EQ.0.D0) THEN 
         THET=0.001D0 
      ELSE 
         THET=WORK(3) 
         IF (THET.GE.1.0D0) THEN 
            WRITE(6,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3) 
            ARRET=.TRUE. 
         END IF 
      END IF 
 
C ---> FNEWT : STOPPING CRITERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1. <--- 
      TOLST=RTOL(1) 
      IF (WORK(4).EQ.0.D0) THEN 
         FNEWT=MAX(10*UROUND/TOLST,MIN(0.03D0,TOLST**0.5D0)) 
      ELSE 
         FNEWT=WORK(4) 
         IF (FNEWT.LE.UROUND/TOLST) THEN 
            WRITE(6,*)' CURIOUS INPUT FOR WORK(4)=',WORK(4) 
            ARRET=.TRUE. 
         END IF 
      END IF 
 
C ---> QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST. <--- 
      IF (WORK(5).EQ.0.D0) THEN 
         QUOT1=1.D0 
      ELSE 
         QUOT1=WORK(5) 
      END IF 
      IF (WORK(6).EQ.0.D0) THEN 
         QUOT2=1.2D0 
      ELSE 
         QUOT2=WORK(6) 
      END IF 
      IF (QUOT1.GT.1.0D0.OR.QUOT2.LT.1.0D0) THEN 
         WRITE(6,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2 
         ARRET=.TRUE. 
      END IF 
C -------------------------------------------------------  
 
C ---->    GRID WITH DISCONTINUITIES  <---- 
      XURO=100*UROUND*ABS(XEND) 
      IF (NGRID.GT.0) THEN 
         IF (GRID(NGRID)-XEND.GE.XURO) THEN 
            IF(IOUT.GT.0) WRITE(6,*) 
     &          ' GRID(NGRID) HAS TO BE <= XEND' 
            ARRET=.TRUE. 
         ENDIF 
         IF (ABS(GRID(NGRID)-XEND).GE.XURO) NGRID=NGRID+1 
      ELSE 
         NGRID=NGRID+1 
      END IF 
      GRID(NGRID)=XEND 
C -------------------------------------------------------  
 
C -------> MAXIMAL STEP SIZE <------- 
      IF (WORK(7).EQ.0.D0) THEN 
         DO I=1,NGRID 
          IF (GRID(I).GT.X) THEN  
           IGRID=I 
           GO TO 2 
          END IF 
         END DO 
 2       CONTINUE 
         HMAX=GRID(IGRID)-X 
      ELSE 
         HMAX=WORK(7) 
      END IF  
 
C ------->  FACL,FACR     PARAMETERS FOR STEP SIZE SELECTION <------- 
      IF(WORK(8).EQ.0.D0)THEN 
         FACL=5.D0 
      ELSE 
         FACL=1.D0/WORK(8) 
      END IF 
      IF(WORK(9).EQ.0.D0)THEN 
         FACR=1.D0/8.0D0 
      ELSE 
         FACR=1.D0/WORK(9) 
      END IF 
      IF (FACL.LT.1.0D0.OR.FACR.GT.1.0D0) THEN 
            WRITE(6,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9) 
            ARRET=.TRUE. 
      END IF 
C ------->  PARAMETER FOR THE CONTROL OF DENSE OUTPUT <------- 
      ALPHA=WORK(10) 
      IF (ALPHA.LT.0.D0.OR.ALPHA.GT.1.D0) THEN 
            WRITE(6,*)' CURIOUS INPUT WORK(10)=',WORK(10) 
            ARRET=.TRUE. 
      END IF 
C ------->   PARAMETER FOR CONTROLLING THE SEARCH OF BP <-------
      TCKBP=WORK(11)
	  IF (TCKBP.LE.0.D0) THEN
	      TCKBP=5.D0
	  END IF	       
C *** *** *** *** *** *** *** *** *** *** *** *** *** 
C         COMPUTATION OF ARRAY ENTRIES 
C *** *** *** *** *** *** *** *** *** *** *** *** *** 
C ---- IMPLICIT, BANDED OR NOT ? 
      IMPLCT=IMAS.NE.0 
       NEUTRAL=IMAS.EQ.2
      JBAND=MLJAC.LT.NM1 
C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- 
C -- JACOBIAN  AND  MATRICES E1, E2 
      IF (JBAND) THEN 
         LDJAC=MLJAC+MUJAC+1 
         LDE1=MLJAC+LDJAC 
      ELSE 
         MLJAC=NM1 
         MUJAC=NM1 
         LDJAC=NM1 
         LDE1=NM1 
      END IF 
C -- MASS MATRIX 
      IF (IMPLCT) THEN 
          IF (MLMAS.NE.NM1) THEN 
              LDMAS=MLMAS+MUMAS+1 
              IF (JBAND) THEN 
                 IJOB=4 
              ELSE 
                 IJOB=3 
              END IF 
          ELSE 
              LDMAS=NM1 
              IJOB=5 
          END IF 
C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC" 
          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN 
             WRITE (6,*) 'BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF 
     & "JAC"' 
            ARRET=.TRUE. 
          END IF 
      ELSE 
          LDMAS=0 
          IF (JBAND) THEN 
             IJOB=2 
          ELSE 
             IJOB=1 
             IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7 
          END IF 
      END IF 
      LDMAS2=MAX(1,LDMAS) 
C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN 
      IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN 
         WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH  
     &FULL JACOBIAN' 
         ARRET=.TRUE. 
      END IF 
 
C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 
      IF (ARRET) THEN 
         IDID=-1 
         RETURN 
      END IF 
 
C     NUMERICAL KERNEL
      WRITE(6,*) 'INTEGRATION...'
C -------- CALL TO CORE INTEGRATOR ------------ 
      CALL RADCOR(N,X,Y,XEND,H,FCN,PHI,ARGLAG,RTOL,ATOL,ITOL, 
     &   JAC,IJAC,MLJAC,MUJAC,JACLAG,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID, 
     &   NMAX,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN, 
     &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1, 
     &   IMPLCT,NEUTRAL,NDIMN,JBAND,LDJAC,LDE1,LDMAS2, 
     &   NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,NFULL,RPAR,IPAR, 
     &   IPAST,GRID,LGRID,NRDS,NLAGS,NJACL, 
     &   NGRID,IEFLAG,WORK(7),TCKBP,ALPHA,ISWJL) 
      IWORK(13)=NFULL 
      IWORK(14)=NFCN 
      IWORK(15)=NJAC 
      IWORK(16)=NSTEP 
      IWORK(17)=NACCPT 
      IWORK(18)=NREJCT 
      IWORK(19)=NDEC 
      IWORK(20)=NSOL 
C -------- RESTORE TOLERANCES 
      EXPM=1.0D0/EXPM 
      IF (ITOL.EQ.0) THEN 
              QUOT=ATOL(1)/RTOL(1) 
              RTOL(1)=(10.0D0*RTOL(1))**EXPM 
              ATOL(1)=RTOL(1)*QUOT 
      ELSE 
          DO I=1,N 
              QUOT=ATOL(I)/RTOL(I) 
              RTOL(I)=(10.0D0*RTOL(I))**EXPM 
              ATOL(I)=RTOL(I)*QUOT 
          END DO 
      END IF 
C ----------- RETURN ----------- 
      RETURN 
      END 
C 
C     END OF SUBROUTINE RADAR5 
C 
C *********************************************************** 
C 
 
      SUBROUTINE RADCOR(N,X,Y,XEND,H,FCN,PHI,ARGLAG,RTOL,ATOL,ITOL, 
     &   JAC,IJAC,MLJAC,MUJAC,JACLAG,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID, 
     &   NMAX,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN, 
     &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1, 
     &   IMPLCT,NEUTRAL,NDIMN,BANDED,LDJAC,LDE1,LDMAS, 
     &   NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,NFULL,RPAR,IPAR, 
     &   IPAST,GRID,LGRID,NRDS,NLAGS,NJACL, 
     &   NGRID,IEFLAG,WORK7,TCKBP,ALPHA,ISWJL) 
C ---------------------------------------------------------- 
C     CORE INTEGRATOR FOR RADAR5 
C     PARAMETERS SAME AS IN RADAR5 WITH WORKSPACE ADDED  
C ----------------------------------------------------------  
C         DECLARATIONS  
C ----------------------------------------------------------  
C     use definitions 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      INTEGER, PARAMETER :: DP=kind(1D0) 
      REAL(kind=DP), dimension(1), intent(inout) :: Y 
      REAL(kind=DP), dimension(:), allocatable ::  
     &                              Z1,Z2,Z3,Y0,SCAL,F1,F2,F3,CONT,PAST 
      REAL(kind=DP), dimension(:), allocatable :: BPV,UCONT 
      REAL(kind=DP), dimension(:,:), allocatable ::  
     &                              FJAC,FJACS,FMAS,E1,E2R,E2I 
      REAL(kind=DP), dimension(1), intent(inout) :: ATOL,RTOL 
      REAL(kind=DP), dimension(1), intent(in) :: RPAR 
      INTEGER, dimension(1), intent(in) :: IPAR
      INTEGER, dimension(NRDS+2*N), intent(inout) :: IPAST 
      REAL(kind=DP), dimension(LGRID), intent(inout) :: GRID 
      REAL(kind=DP), dimension(:,:), allocatable :: FJACL,XLAG 
      REAL(kind=DP), dimension(:), allocatable :: FJACLAG,ZL  
      INTEGER, dimension(:,:), allocatable :: ICOUN 
      INTEGER, dimension(:), allocatable :: IVL,IVE,IVC,ILS 
      INTEGER, dimension(:), allocatable :: IP1,IP2,IPHES,IPJ 
      LOGICAL FLAGS,FLAGN,FLAGUS 
      LOGICAL QUADR 
      LOGICAL BPC,BPD,BPDMEM,LEFT 
      LOGICAL REPEAT
C ----> COMMON BLOCKS <---- 
      COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
      COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2 
      COMMON /BPLOG/FIRST,LAST,REJECT,BPD 
      COMMON /BPCOM/BPP,ILBP,LEFT 
      COMMON /LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG 
C ---- 
      LOGICAL REJECT,FIRST,IMPLCT,NEUTRAL,PROJECT,BANDED,CALJAC 
      LOGICAL STARTN,CALHES,CALJACL,CALLAG 
      LOGICAL INDEX1,INDEX2,INDEX3,LAST,PRED 
      EXTERNAL FCN,PHI,ARGLAG 
C *** *** *** *** *** *** *** 
C  INITIALISATIONS 
C *** *** *** *** *** *** *** 
      ALLOCATE (Z1(N),Z2(N),Z3(N),Y0(N),SCAL(N),F1(N),F2(N),F3(N)) 
      ALLOCATE(BPV(10000)) 
      ALLOCATE (FJAC(LDJAC,N),ZL(3*N)) 
      IF (IMPLCT) ALLOCATE(FMAS(LDMAS,NM1)) 
      ALLOCATE (IP1(NM1),IP2(NM1),IPHES(NM1)) 
      ALLOCATE (E1(LDE1,NM1),E2R(LDE1,NM1),E2I(LDE1,NM1)) 
      ALLOCATE (PAST(MXST*IDIF)) 
      IF (NLAGS.GT.0) THEN 
       ALLOCATE (FJACS(LDJAC,N),FJACLAG(NJACL)) 
       ALLOCATE (IVL(NJACL),IVE(NJACL),IVC(NJACL), 
     &           ILS(2*NLAGS+NJACL),ICOUN(3,NLAGS)) 
       IF (ISWJL.NE.1) THEN 
        ALLOCATE (IPJ(3*N),FJACL(3*N,3*N)) 
       END IF 
       ALLOCATE (XLAG(3,NLAGS)) 
      END IF 
  
C     AMPLITUDE OF CONT 
      LRC=4*N 
      ALLOCATE (CONT(LRC)) 
      ALLOCATE (UCONT(LRC+2)) 
C --- 
       OPEN(8,FILE='radar5.log')
       REWIND 8

C -------------------------------------------------  
      BPC=.FALSE. 
      BPD=.FALSE. 
      BPDMEM=.FALSE. 
      QUADR=.FALSE. 
      ILBP=0
       
C --- INITIAL PREPARATIONS 
      CALLAG=.FALSE. 
      IACT=1 
      IPOS=1 
      DO I=1,NLAGS 
       IPAST(I+2*N)=1 
      END DO 
      X0B=X 
      DO I=1,NGRID 
       IF (GRID(I).GT.X0B) THEN  
        IGRID=I 
        GO TO 2 
       END IF 
      END DO 
 2    XEND=GRID(IGRID) 
      IBP=1 
      BPV(1)=X0B 
      BTOL=1.D1*UROUND 
      KMAX=10
      IMANT=0
       
C --- GUSTAFFSON TECHNIQUE AFTER BREAKING POINTS IS NOT APPLIED        
      FLAGUS=.FALSE. 
      ERRACC=1.D0
 
      IRTRN=2 
      CALL FCN(N,X,Y,Y0,ARGLAG,PHI,RPAR,IPAR,PAST,IPAST,NRDS) 
      IRTRN=1 
 
C     TOLMIN
      IF (ITOL.EQ.0) THEN
	 RTOLM=RTOL(1)
	ELSE
       RTOLM=RTOL(1)
       DO I=2,N
        IF (RTOL(I).LT.RTOLM) RTOLM=RTOL(I)
	 END DO
	END IF
 
C -------- CHECK THE INDEX OF THE PROBLEM -----  
      INDEX1=NIND1.NE.0 
      INDEX2=NIND2.NE.0 
      INDEX3=NIND3.NE.0 
C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- 
      IF (IMPLCT) CALL MAS(NM1,FMAS,LDMAS,RPAR,IPAR) 
C ---------> REQUIRED CONSTANTS <--------- 
      SQ6=DSQRT(6.D0) 
      C1=(4.D0-SQ6)/10.D0 
      C2=(4.D0+SQ6)/10.D0 
      C1M1=C1-1.D0 
      C2M1=C2-1.D0 
      C1MC2=C1-C2 
      CQ1=(2.D0+3.D0*SQ6)/6.D0 
      CQ2=(2.D0-3.D0*SQ6)/6.D0 
      CQ3=1.D0/3.D0 
      CL1=10.D0/(6.D0+SQ6)                   
      CL2=0.D0                 
      CL3=(-4.D0+SQ6)/(6.D0+SQ6)  
      CERS=5.D-1 
      CERC=5.D-1 
      CERLQ=1.D-2 
      THRS=100.D0 
      DD1=-(13.D0+7.D0*SQ6)/3.D0 
      DD2=(-13.D0+7.D0*SQ6)/3.D0 
      DD3=-1.D0/3.D0 
      U1=(6.D0+81.D0**(1.D0/3.D0)-9.D0**(1.D0/3.D0))/30.D0 
      ALPH=(12.D0-81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))/60.D0 
      BETA=(81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))*DSQRT(3.D0)/60.D0 
      CNO=ALPH**2+BETA**2 
      U1=1.0D0/U1 
      ALPH=ALPH/CNO 
      BETA=BETA/CNO 
      T11=9.1232394870892942792D-02 
      T12=-0.14125529502095420843D0 
      T13=-3.0029194105147424492D-02 
      T21=0.24171793270710701896D0 
      T22=0.20412935229379993199D0 
      T23=0.38294211275726193779D0 
      T31=0.96604818261509293619D0 
      TI11=4.3255798900631553510D0 
      TI12=0.33919925181580986954D0 
      TI13=0.54177053993587487119D0 
      TI21=-4.1787185915519047273D0 
      TI22=-0.32768282076106238708D0 
      TI23=0.47662355450055045196D0 
      TI31=-0.50287263494578687595D0 
      TI32=2.5719269498556054292D0 
      TI33=-0.59603920482822492497D0 
C 
C     INVERSE OF A 
      IF (NLAGS.GT.0.AND.ISWJL.NE.1) THEN 
       AI11= 3.22474487139158904909864D0 
       AI12= 1.16784008469040549492404D0 
       AI13=-0.25319726474218082618594D0 
       AI21=-3.56784008469040549492404D0 
       AI22= 0.77525512860841095090136D0 
       AI23= 1.05319726474218082618594D0 
       AI31= 5.53197264742180826185942D0 
       AI32=-7.53197264742180826185942D0 
       AI33= 5.00000000000000000000000D0 
      END IF 
C 
      IF (M1.GT.0) IJOB=IJOB+10 
      HMAXN=MIN(HMAX,XEND-X)  
      IF (H.LE.10.D0*UROUND) H=1.0D-6  
      H=MIN(H,HMAXN) 
      HOLD=H 
      REJECT=.FALSE. 
      FIRST=.TRUE. 
      LAST=.FALSE. 
	  NITER=0
      IF ((X+H*1.0001D0-XEND).GE.0.D0) THEN 
         H=XEND-X 
         LAST=.TRUE. 
      END IF 
C ---  INITIALIZATION FOR THE ARRAY PAST    
       DO 3 I=0,MXST-1 
          PAST(1+IDIF*I)=X  
   3   CONTINUE 
       IPA=(MXST-1)*IDIF+1 
       DO J=1,NRDS 
          K=IPAST(J) 
          PAST(J+IPA)=Y(K) 
          PAST(J+1*NRDS+IPA)=0.D0 
          PAST(J+2*NRDS+IPA)=0.D0 
          PAST(J+3*NRDS+IPA)=0.D0 
       END DO 
          PAST(IPA+IDIF-1)=H  
C ---  END OF THE INITIALIZATION      
      FACCON=1.D0 
      CFAC=SAFE*(1+2*NIT) 
      NSING=0 
      XOLD=X 
      IF (IOUT.NE.0) THEN 
          IRTRN=1 
          NRSOL=1 
          XOSOL=XOLD 
          XSOL=X 
          CONT(1:N)=Y(1:N) 
          NSOLU=N 
          HSOL=HOLD
          CALL SOLOUT(NRSOL,XOSOL,XSOL,HSOL,Y,CONT,LRC,NSOLU, 
     &                RPAR,IPAR,IRTRN) 
          IF (IRTRN.LT.0) GOTO 179 
      END IF 
      MLE=MLJAC 
      MUE=MUJAC 
      MBJAC=MLJAC+MUJAC+1 
      MBB=MLMAS+MUMAS+1 
      MDIAG=MLE+MUE+1 
      MDIFF=MLE+MUE-MUMAS 
      MBDIAG=MUMAS+1 
      N2=2*N 
      N3=3*N 
      IF (ITOL.EQ.0) THEN 
          DO I=1,N 
             SCAL(I)=ATOL(1)+RTOL(1)*ABS(Y(I)) 
          END DO 
      ELSE 
          DO I=1,N 
             SCAL(I)=ATOL(I)+RTOL(I)*ABS(Y(I)) 
          END DO 
      END IF 
      HHFAC=H 
      NFCN=NFCN+1 
 
      NFULL=0 
C -------------------------- 
C --- BASIC INTEGRATION STEP   
C -------------------------- 
  10  CONTINUE 
C *** *** *** *** *** *** *** 
C  COMPUTATION OF THE JACOBIAN 
C *** *** *** *** *** *** *** 
C ----------------------- 
      FLAGS=.FALSE. 
      FLAGN=.FALSE. 
C ----------------------- 
      ALOPT=0.D0 
      NJAC=NJAC+1 
	IF (BPD) THEN
	 BPDMEM=.TRUE.
	 BPD=.FALSE.
      END IF
      IF (IJAC.EQ.0) THEN 
C --- COMPUTE JACOBIAN MATRIX NUMERICALLY 
         IF (BANDED) THEN 
C --- JACOBIAN IS BANDED 
            MUJACP=MUJAC+1 
            MD=MIN(MBJAC,M2) 
            DO MM=1,M1/M2+1 
               DO K=1,MD 
                  J=K+(MM-1)*M2 
 12               F1(J)=Y(J) 
                  F2(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J)))) 
                  Y(J)=Y(J)+F2(J) 
                  J=J+MD 
                  IF (J.LE.MM*M2) GOTO 12  
                  CALL FCN(N,X,Y,CONT,ARGLAG,PHI,RPAR,IPAR,PAST,IPAST,
     &                     NRDS) 
                  J=K+(MM-1)*M2 
                  J1=K 
                  LBEG=MAX(1,J1-MUJAC)+M1 
 14               LEND=MIN(M2,J1+MLJAC)+M1 
                  Y(J)=F1(J) 
                  MUJACJ=MUJACP-J1-M1 
                  DO L=LBEG,LEND 
                     FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/F2(J)  
                  END DO 
                  J=J+MD 
                  J1=J1+MD 
                  LBEG=LEND+1 
                  IF (J.LE.MM*M2) GOTO 14 
               END DO 
            END DO 
         ELSE 
C --- JACOBIAN IS FULL 
            DO I=1,N 
               YSAFE=Y(I) 
	         DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE))) 
			 Y(I)=YSAFE+DELT 
                 CALL FCN(N,X,Y,CONT,ARGLAG,PHI, 
     1	              RPAR,IPAR,PAST,IPAST,NRDS)
	         DO J=M1+1,N 
                 FJAC(J-M1,I)=(CONT(J)-Y0(J))/DELT 
               END DO 
               Y(I)=YSAFE 
            END DO
         END IF 
      ELSE 
C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
         CALL JAC(N,X,Y,FJAC,LDJAC,ARGLAG,PHI,RPAR,IPAR,PAST,IPAST,NRDS) 
      END IF
	IF (BPDMEM) THEN
	 BPDMEM=.FALSE.
	 BPD=.TRUE.
      END IF 
      CALJAC=.TRUE. 
      CALHES=.TRUE. 
      JLFLAG=0 
C --- SAVE FJAC 
      IF (NLAGS.GT.0) FJACS=FJAC 
 
C ------------------------------------------------ 
C --- GLOBAL ITERATION BEGINS HERE
C ------------------------------------------------ 
  20  CONTINUE 
 
C *** *** *** *** *** *** *** 
C  STARTING VALUES FOR NEWTON ITERATION 
C *** *** *** *** *** *** *** 
    
      IF (FIRST.OR.STARTN) THEN 
            DO I=1,N 
             Z1(I)=0.D0 
             Z2(I)=0.D0 
             Z3(I)=0.D0 
             F1(I)=0.D0 
             F2(I)=0.D0 
             F3(I)=0.D0 
C 
             A1=Y(I) 
             ZL(I)=A1 
             ZL(I+N)=A1 
             ZL(I+N2)=A1 
            END DO 
      ELSE 
         C3Q=H/HOLD 
         C1Q=C1*C3Q 
         C2Q=C2*C3Q 
         DO I=1,N 
            A1=Y(I) 
            AK1=CONT(I+N) 
            AK2=CONT(I+N2) 
            AK3=CONT(I+N3) 
            Z1I=C1Q*(AK1+(C1Q-C2M1)*(AK2+(C1Q-C1M1)*AK3)) 
            Z2I=C2Q*(AK1+(C2Q-C2M1)*(AK2+(C2Q-C1M1)*AK3)) 
            Z3I=C3Q*(AK1+(C3Q-C2M1)*(AK2+(C3Q-C1M1)*AK3)) 
            Z1(I)=Z1I 
            Z2(I)=Z2I 
            Z3(I)=Z3I 
C 
            ZL(I)=A1+Z1I 
            F1(I)=TI11*Z1I+TI12*Z2I+TI13*Z3I 
            F2(I)=TI21*Z1I+TI22*Z2I+TI23*Z3I 
            F3(I)=TI31*Z1I+TI32*Z2I+TI33*Z3I 
         END DO 
            IF (NLAGS.GT.0) THEN 
             DO I=1,N 
              A1=Y(I) 
              ZL(I+N)=A1+Z2(I) 
              ZL(I+N2)=A1+Z3(I) 
             END DO 
            END IF 
      END IF 
C --- 
      IF (JLFLAG.EQ.0) THEN 
       IF (NLAGS.EQ.0) THEN 
        IF (.NOT.(REJECT.OR.LAST)) THEN  
         IF (IMANT.EQ.1) GOTO 30 
        END IF 
        GOTO 22  
       END IF 
C --------------------- 
C      INITIALIZATION   
C 
C ---  LAGS CYCLE --- C 
C --- 
       CALJACL=.FALSE. 
       X1=X+C1*H 
       X2=X+C2*H 
       X3=X+H 
       DO IL=1,NLAGS 
        XLAG(1,IL)=0.D0 
        XLAG(2,IL)=0.D0 
        XLAG(3,IL)=0.D0 
        ICOUN(1,IL)=0 
        ICOUN(2,IL)=0 
        ICOUN(3,IL)=0 
       END DO 
C ---  LOOP ON LAG TERMS 
       DO IL=1,NLAGS 
C ---   DELAYED ARGUMENTS ARE COMPUTED 
        XLAG(1,IL)=ARGLAG(IL,X1,N,ZL,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
           IF (XLAG(1,IL).GT.X) ICOUN(1,IL)=1 
        XLAG(2,IL)=ARGLAG(IL,X2,N,ZL(N+1),RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
           IF (XLAG(2,IL).GT.X) ICOUN(2,IL)=1 
        XLAG(3,IL)=ARGLAG(IL,X3,N,ZL(N2+1),RPAR,IPAR,PHI,
     &             PAST,IPAST,NRDS)
           IF (XLAG(3,IL).GT.X) ICOUN(3,IL)=1 
        IF (ICOUN(1,IL)+ICOUN(2,IL)+ICOUN(3,IL).GE.1) CALJACL=.TRUE.
       END DO 
 
       IF (CALJACL) THEN 
        CALL JACLAG(N,X,Y,FJACLAG,ARGLAG,PHI,IVE,IVC,IVL, 
     &                  RPAR,IPAR,PAST,IPAST,NRDS) 
        IF (.NOT.CALLAG) THEN 
         CALLAG=.TRUE. 
C --     ORDERING STEP 
         LL=2*NLAGS+1 
         DO L=1,NLAGS 
          NL=0 
          DO I=1,NJACL 
           IF (IVL(I).EQ.L) THEN 
            ILS(LL)=I 
            NL=NL+1 
            LL=LL+1 
           END IF 
          END DO 
          ILS(2*L-1)=NL 
          ILS(2*L)=LL-1 
         END DO 
        END IF 
C 
        DO IL=1,NLAGS 
          DO IS=1,3 
           SELECT CASE (IS) 
            CASE (1) 
             XACT=X1 
            CASE (2) 
             XACT=X2 
            CASE (3) 
             XACT=X3 
           END SELECT  
           IF (XLAG(IS,IL).GT.XACT) THEN 
             IF (IOUT.EQ.2) WRITE (6,*) 
     &          ' WARNING!: ADVANCED ARGUMENTS ARE USED AT X= ',XACT 
             XLAG(IS,IL)=XACT 
           END IF 
          END DO 
 
C ------ JACOBIAN MAINTAINS THE TENSOR STRUCTURE
C ------ UPDATING CONDITION
          ALOPT=0.D0 
          IF (ICOUN(1,IL).EQ.1) THEN 
             S1=DIM(XLAG(1,IL),X)/H 
             ALOPT=(-1.D0+S1)*S1*(-13.D0-7.D0*Sqrt(6.D0) + 
     &        5.D0*(2.D0+3.D0*Sqrt(6.D0))*S1) 
          END IF 
          IF (ICOUN(2,IL).EQ.1) THEN 
             S2=DIM(XLAG(2,IL),X)/H 
             ALOPT=ALOPT-(-1+S2)*S2* 
     &       (13.D0-7.D0*Sqrt(6.D0)+5.D0*(-2.D0+3.D0*Sqrt(6.D0))*S2) 
          END IF 
          IF (ICOUN(3,IL).EQ.1) THEN 
             S3=DIM(XLAG(3,IL),X)/H 
             ALOPT=ALOPT+S3*(1.D0 - 8.D0*S3 + 10.D0*S3**2) 
          END IF 
          ALOPT=ALOPT/9.D0         
          
C         OPTIMAL COEFFICIENT (W.R.T. FROBENIUS NORM)
C         JACLAG ~= ALOPT*I 
C 
C         ACTIVATES IF ALOPT DIFFERENT FROM ZERO
          IF (ABS(ALOPT).GE.1.D-8) THEN 
            NL =ILS(2*IL-1) 
            ILE=ILS(2*IL) 
            DO K=1,NL 
              KK=ILS(ILE-K+1) 
              IK=IVE(KK) 
              JK=IVC(KK) 
              FJAC(IK,JK)=FJAC(IK,JK)+ALOPT*FJACLAG(KK)
            END DO 
          ELSE 
            CALJACL=.FALSE. 
          END IF 
        END DO 
C 
        ELSE 
C ---   FACTORIZATION IS CONSERVED  
        END IF 
        IF (.NOT.(REJECT.OR.LAST.OR.CALJACL)) THEN  
          IF (IMANT.EQ.1) GOTO 30 
        END IF 
        GOTO 22  
C --- 
      ELSE IF (JLFLAG.EQ.1) THEN 
C ---  THE FULL ITERATION MAKES USE OF THE SAME STEPSIZE OF THE SIMPLIFIED
        GOTO 23 
C --- 
      END IF 
 
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
C - - THE DIFFERENCE IN THE SOLUTION OF THE SYSTEM STARTS HERE 
C ------------------------------- - - - - - - - - - - - - - - - 
C --- SIMPLE NEWTON ITERATION 
C ------------------------------- 
  22  CONTINUE 
C 
C --- COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS 
      FAC1=U1/H 
      ALPHN=ALPH/H 
      BETAN=BETA/H 
C --- RK JACOBIAN FACTORIZATION 
      CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, 
     &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES) 
      IF (IER.NE.0) THEN 
          GOTO 78 
      END IF 
      CALL DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, 
     &            M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB) 
      IF (IER.NE.0) THEN 
          GOTO 78 
      END IF 
      NDEC=NDEC+1 
C --- 
  30  CONTINUE 
C --- UPDATE NSTEP 
      NSTEP=NSTEP+1 
      IF (NSTEP.GT.NMAX) GOTO 178 
      IF (0.1D0*H.LE.ABS(X)*UROUND) GOTO 177 
          IF (INDEX2) THEN 
             DO I=NIND1+1,NIND1+NIND2 
                SCAL(I)=SCAL(I)/HHFAC 
             END DO 
          END IF 
          IF (INDEX3) THEN 
             DO I=NIND1+NIND2+1,NIND1+NIND2+NIND3 
                SCAL(I)=SCAL(I)/(HHFAC*HHFAC) 
             END DO 
          END IF 
      XPH=X+H 
 
C *** *** *** *** *** *** *** 
C  LOOP FOR THE SIMPLE NEWTON ITERATION 
C *** *** *** *** *** *** *** 
C              -----------      
            NEWT=0 
            IMANT=0 
C ----------------------- 
            FLAGS=.FALSE. 
            FLAGN=.FALSE. 
C ----------------------- 
            FACCON=MAX(FACCON,UROUND)**0.8D0 
            THETA=ABS(THET) 
C ------------------------------------------------------- 
C           REFERENCE POINT FOR THE SIMPLE AZIONE SIMPLE 
C ------------------------------------------------------- 
  40        CONTINUE  
C ------------------------- 
            IF (FLAGS) THEN 
             FLAGN=.TRUE. 
C ---        CALLS SUBROUTINE LAGR5 
C ------------------------------------------------------------- 
C ---------- DYNAMIC UPDATE OF INTERPOLANT (in PAST).  
C ------------------------------------------------------------- 
             DO J=1,NRDS 
                I=IPAST(J) 
                PAST(J+IACT)=Y(I)+Z3(I) 
                 Z2I=Z2(I) 
                 Z1I=Z1(I) 
                A1=(Z2I-Z3(I))/C2M1 
                PAST(J+1*NRDS+IACT)=A1 
                 AK=(Z1I-Z2I)/C1MC2 
                 ACONT3=Z1I/C1 
                 ACONT3=(AK-ACONT3)/C2 
                A2=(AK-A1)/C1M1  
                PAST(J+2*NRDS+IACT)=A2 
                PAST(J+3*NRDS+IACT)=A2-ACONT3 
             ENDDO 
C ---        UPDATE DI PAST 
             PAST(IACT)=X 
C            LEFT ENDPOINT 
             PAST(IACT+IDIF-1)=H 
C            STEPSIZE 
            END IF 
C ---------------- 
            IF (NEWT.GE.NIT) THEN 
             INREJ=2 
             GOTO 421 
            END IF 
C ----------------------------------- 
C ---     COMPUTE THE RIGHT-HAND SIDE 
            DO I=1,N 
             A1=Y(I) 
             ZL(I)=A1+Z1(I) 
             ZL(I+N)=A1+Z2(I) 
             ZL(I+N2)=A1+Z3(I) 
            END DO 
C           COMPUTATION OF STAGE VALUES
            CALL FCN(N,X+C1*H,ZL,Z1,ARGLAG,PHI,RPAR,IPAR,PAST,
     &               IPAST,NRDS) 
            CALL FCN(N,X+C2*H,ZL(N+1),Z2,ARGLAG,PHI,RPAR,IPAR,PAST,
     &               IPAST,NRDS)
            IF (BPD) THEN
C ---------------------------------------------------------------------- 
C ---         A BREAKING POINT HAS BEEN DETECTED
C ---------------------------------------------------------------- 
              LEFT=.FALSE. 
              DBP=ARGLAG(ILBP,X+C2*H,N,ZL(N+1),RPAR,IPAR,PHI, 
     &                   PAST,IPAST,NRDS) 
              IF (DBP.LT.BPP) LEFT=.TRUE. 
            ELSE IF (LAST) THEN 
C ---         DISC. FLAG 
              XX=(X+H)*(1.D0-BTOL) 
              DO I=1,N 
                 ZL(N2+I)=(1.D0-BTOL)*ZL(N2+I)+BTOL*ZL(N+I) 
              END DO 
              CALL FCN(N,XX,ZL(N2+1),Z3,ARGLAG,PHI,RPAR,IPAR,PAST, 
     &                 IPAST,NRDS) 
              GO TO 42 
            END IF 
            CALL FCN(N,X+H,ZL(N2+1),Z3,ARGLAG,PHI,RPAR,IPAR,PAST, 
     &               IPAST,NRDS) 
 42         CONTINUE
    
            NFCN=NFCN+3 
C ----------------------------------- 
C ---     SOLVE THE LINEAR SYSTEMS 
           DO I=1,N 
              A1=Z1(I) 
              A2=Z2(I) 
              A3=Z3(I) 
              Z1(I)=TI11*A1+TI12*A2+TI13*A3 
              Z2(I)=TI21*A1+TI22*A2+TI23*A3 
              Z3(I)=TI31*A1+TI32*A2+TI33*A3 
           END DO 
            
        CALL SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, 
     &          M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3, 
     &          F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB) 
            NSOL=NSOL+1 
            NEWT=NEWT+1 
C ---       NORM OF DY 
		  DYNO=0.D0 
	      DO I=1,N
            END DO
		  DO I=1,N 
               DENOM=SCAL(I) 
               DYNO=DYNO+(Z1(I)/DENOM)**2+(Z2(I)/DENOM)**2 
     &          +(Z3(I)/DENOM)**2 
            END DO 
            DYNO=DSQRT(DYNO/N3) 
C -------------------------------------------------------- 
C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE 
C --------------------------------------------------------
            IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN 
                THQ=DYNO/DYNOLD 
                IF (NEWT.EQ.2) THEN 
                   THETA=THQ 
                ELSE 
                   THETA=SQRT(THQ*THQOLD) 
                END IF 
                THQOLD=THQ 
                INREJ=0 
C --- 1 
                IF (THETA.LT.0.99D0) THEN 
                    FACCON=THETA/(1.0D0-THETA) 
                    DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT 
                    IF (DYTH.GE.1.0D0) INREJ=1 
C ----------------  SLOW CONVERGENCE --- 
                ELSE 
                    INREJ=2 
C ----------------  DIVERGENCE --- 
                END IF 
C --- 1 
            ELSE 
C ---           NEWT=1 
                INREJ=0 
            END IF 
C ---------------------------------------- 
C ------------------ THE STEP IS REPEATED   
C ---------------------------------------- 
C ---- 2 
 421        IF (INREJ.GT.0) THEN 
C 
C ----- 3 
             IF (.NOT.(BPC.OR.FIRST)) THEN    
              IF (BPD) THEN
C ---          BP IS WRONG                      
               IBP=IBP-1
               BPD=.FALSE.
              END IF
              HP=H*0.99D0
              CALL BPDTCT(N,X,HP,Y,ARGLAG,RPAR,IPAR,UCONT,GRID,NLAGS, 
     &                    FIRST,LAST,XEND,IGRID,BPV,IBP,ILBP,BPP,BPD, 
     &                    KMAX,PHI,PAST,IPAST,NRDS) 
            
              BPC=.TRUE. 
C ------ 4 
              IF (BPD) THEN 
               H=HP 
               LAST=.TRUE.
              ELSE IF (INREJ.EQ.1) THEN 
               QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH)) 
               HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT)) 
               H=HHFAC*H 
               LAST=.FALSE. 
C
              ELSE IF (INREJ.EQ.2) THEN 
               H=H*0.55D0 
               HHFAC=0.55D0 
               LAST=.FALSE. 
C
              END IF 
C ------ 4 
             ELSE   
C ----- 3 
              IF (BPD) THEN  
C              BP WRONG
               IBP=IBP-1 
               BPD=.FALSE. 
               LAST=.FALSE. 
              END IF 
C ------------------------------------------- 
              IF (.NOT.CALJACL.OR.ISWJL.EQ.1) THEN  
C ---          THERE ARE NOT SMALL DELAYS 
               IF (REJECT.AND.CALJAC) THEN  
                H=H*0.12D0 
                HHFAC=0.12D0 
               ELSE IF (INREJ.EQ.1) THEN 
                QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH)) 
                HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT)) 
                H=HHFAC*H 
	     	LAST=.FALSE. 
               ELSE IF (INREJ.EQ.2) THEN 
                H=H*0.55D0 
                HHFAC=0.55D0 
                LAST=.FALSE. 
               END IF 
              ELSE 
C ---          CALJACL IS TRUE: FULL ITEARATION IS DONE          
C ---          THE STEPSIZE DOES NOT CHANGET
               JLFLAG=1 
              END IF 
             END IF 
C ----- 3 
             REJECT=.TRUE. 
             IF (CALJAC) GOTO 20 
             GOTO 10 
            END IF 
C ---- 2 
C -------------------------------------------------------- 
            DYNOLD=MAX(DYNO,UROUND)
            DO I=1,N 
               F1I=F1(I)+Z1(I) 
               F2I=F2(I)+Z2(I) 
               F3I=F3(I)+Z3(I) 
               F1(I)=F1I 
               F2(I)=F2I 
               F3(I)=F3I 
               Z1(I)=T11*F1I+T12*F2I+T13*F3I 
               Z2(I)=T21*F1I+T22*F2I+T23*F3I 
               Z3I=T31*F1I+    F2I 
               Z3(I)=Z3I 
C ---          APPROX DELLA SOLUZIONE 
               ZL(I+N2)=Y(I)+Z3I 
            END DO 
C -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
            IF (NEWT.EQ.1.OR.FACCON*DYNO.GT.FNEWT) THEN 
C ---        NEWTON PROCEDE
             GOTO 40 
            END IF 

C ----------------------------------------------------------------
C ---       ITERATIVE CORRECTION OF THE BREAKING POINT
C ----------------------------------------------------------------
            IF (BPD) THEN
	     HNEWT=H
	     NITER = NITER+1
	     IF (NITER.GT.NIT) THEN  
C              BP WRONG                 
               IBP=IBP-1 
               BPD=.FALSE. 
               LAST=.FALSE. 
               H=H*0.55D0 
               HHFAC=0.55D0
	       REJECT=.TRUE.
	       NITER=0 
               IF (CALJAC) GOTO 20 
               GOTO 10   
             END IF  
		    
             CALL BPACC(N,X,H,Y,ARGLAG,RPAR,IPAR,Z1,Z2,Z3,FIRST,
     &                  BPV,IBP,ILBP,BPP,KMAX,PHI,PAST,IPAST,NRDS) 
             IF (ABS(H-HNEWT)/HNEWT.GE.MAX(BTOL,RTOLM*1.D-2)) THEN
              GOTO 20
C             REF POINT
             ELSE
	        H=HNEWT
	        NITER=0
             END IF
            END IF  
C ----------------------------------------------------------------


C *** *** *** *** *** *** *** *** *** *** *** 
C END LOOP 
C *** *** *** *** *** *** *** *** *** *** *** 
      GOTO 55 
C *** *** *** *** *** 
C 
 
C --- FULL NEWTON ITERATION 
  23  CONTINUE 
      NFULL =NFULL +1 
C --- 
 
C ---------------------------------------------------- 
C --- ALTERNATIVE FULL NEWTON JACOB. INTEGRATION STEP   
C ---------------------------------------------------- 
C 
      DO I=1,N 
       DO J=1,N 
        FJACL(I+N,J)=0.D0 
        FJACL(I+2*N,J)=0.D0 
        FJACL(I+2*N,J+N)=0.D0 
        FJACL(I,J+N)=0.D0 
        FJACL(I,J+2*N)=0.D0 
        FJACL(I+N,J+2*N)=0.D0 
        FIJ=FJACS(I,J) 
        FJACL(I,J)=FIJ 
        FJACL(I+N,J+N)=FIJ 
        FJACL(I+2*N,J+2*N)=FIJ 
       END DO 
      END DO 
 
      QUADR=FIRST 
      IF (.NOT.QUADR) THEN 
       DL1=C1*(C1-C2)*(C1-1.D0) 
       DL2=C2*(C2-C1)*(C2-1.D0) 
       DL3=(1.D0-C1)*(1.D0-C2) 
      ELSE 
       DL1=(C1-C2)*(C1-1.D0) 
       DL2=(C2-C1)*(C2-1.D0) 
       DL3=(1.D0-C1)*(1.D0-C2) 
      ENDIF 
      DO IL=1,NLAGS 
       DO I=1,3 
            XL=XLAG(I,IL) 
             IF (XL.GT.X) THEN 
CCC 
CCC            DERIVATIVES OF THE COLLOCATION POLYNOMIAL WRT STAGES 
CCC            d u/d Y_k = L_k(xlag_i) 
CCC                                    
               IF (.NOT.QUADR) THEN 
                DCOLI1=((XL-X)/H)*((XL-X2)/H)*((XL-X3)/H)/DL1 
                DCOLI2=((XL-X)/H)*((XL-X1)/H)*((XL-X3)/H)/DL2 
                DCOLI3=((XL-X)/H)*((XL-X1)/H)*((XL-X2)/H)/DL3 
               ELSE 
                DCOLI1=((XL-X2)/H)*((XL-X3)/H)/DL1 
                DCOLI2=((XL-X1)/H)*((XL-X3)/H)/DL2 
                DCOLI3=((XL-X1)/H)*((XL-X2)/H)/DL3 
               END IF 
 
C -----------> JACOBIAN IS UPDATED 
               NL=ILS(2*IL-1) 
               ILE=ILS(2*IL) 
 
C -----------> FULL JACOBIAN MATRIX FJACL 
               DO K=1,NL 
                 KK=ILS(ILE-K+1) 
                 IK=IVE(KK) 
                 JK=IVC(KK) 
                 FJLK=FJACLAG(KK) 
                 IKI=IK+(I-1)*N 
C 
                  FJACL(IKI,JK)=FJACL(IKI,JK)+FJLK*DCOLI1 
                  FJACL(IKI,JK+N)=FJACL(IKI,JK+N)+FJLK*DCOLI2    
                  FJACL(IKI,JK+2*N)=FJACL(IKI,JK+2*N)+FJLK*DCOLI3 
               END DO 
             END IF 
       END DO            
CCC 
CCC --> NLAGS 
      END DO 
CCC <-- 
 
      AI11H   =-AI11/H 
      AI12H   =-AI12/H 
      AI13H   =-AI13/H 
      AI21H   =-AI21/H 
      AI22H   =-AI22/H 
      AI23H   =-AI23/H 
      AI31H   =-AI31/H 
      AI32H   =-AI32/H 
      AI33H   =-AI33/H 
 
C --- FJACL 
      IF (IMPLCT) THEN 
       DO I1=1,N 
         DO J1=1,N 
           FJACL(I1,J1)=FJACL(I1,J1)+AI11H*FMAS(I1,J1) 
           FJACL(I1,J1+N)=FJACL(I1,J1+N)+AI12H*FMAS(I1,J1) 
           FJACL(I1,J1+2*N)=FJACL(I1,J1+2*N)+AI13H*FMAS(I1,J1) 
 
           FJACL(I1+N,J1)=FJACL(I1+N,J1)+AI21H*FMAS(I1,J1) 
           FJACL(I1+N,J1+N)=FJACL(I1+N,J1+N)+AI22H*FMAS(I1,J1) 
           FJACL(I1+N,J1+2*N)=FJACL(I1+N,J1+2*N)+AI23H*FMAS(I1,J1) 
 
           FJACL(I1+2*N,J1)=FJACL(I1+2*N,J1)+AI31H*FMAS(I1,J1) 
           FJACL(I1+2*N,J1+N)=FJACL(I1+2*N,J1+N)+AI32H*FMAS(I1,J1) 
           FJACL(I1+2*N,J1+2*N)=FJACL(I1+2*N,J1+2*N)+AI33H*FMAS(I1,J1) 
         END DO 
       END DO 
      ELSE 
C ---  EXPLICIT CASE 
       DO I1=1,N 
         FJACL(I1,I1)=FJACL(I1,I1)+AI11H 
         FJACL(I1,I1+N)=FJACL(I1,I1+N)+AI12H 
         FJACL(I1,I1+2*N)=FJACL(I1,I1+2*N)+AI13H 
 
         FJACL(I1+N,I1)=FJACL(I1+N,I1)+AI21H 
         FJACL(I1+N,I1+N)=FJACL(I1+N,I1+N)+AI22H 
         FJACL(I1+N,I1+2*N)=FJACL(I1+N,I1+2*N)+AI23H 
 
         FJACL(I1+2*N,I1)=FJACL(I1+2*N,I1)+AI31H 
         FJACL(I1+2*N,I1+N)=FJACL(I1+2*N,I1+N)+AI32H 
         FJACL(I1+2*N,I1+2*N)=FJACL(I1+2*N,I1+2*N)+AI33H 
       END DO 
      END IF 
C ----- FACTORIZATION OF THE FULL JACOBIAN 
      CALL DEC(3*N,3*LDJAC,FJACL,IPJ,IER)  
      IF (IER.NE.0) THEN 
       GOTO 78 
      END IF 
CCC --->                                
      NDEC=NDEC+1 
  33  CONTINUE 
C --- 
C --- R: EVERY STEP STARTS WITH A SIMPLE ITERATION
      NSTEP=NSTEP+1 
      IF (NSTEP.GT.NMAX) GOTO 178 
      IF (0.1D0*H.LE.ABS(X)*UROUND) GOTO 177 
C --- 
      XPH=X+H 
 
C *** *** *** *** *** *** *** 
C  LOOP FOR NEWTON ITERATION 
C *** *** *** *** *** *** *** 
C              -----------      
            NEWT=0 
C ----------------------- 
            FLAGS=.FALSE. 
            FLAGN=.FALSE. 
C ----------------------- 
            FACCON=MAX(FACCON,UROUND)**0.8D0 
            THETA=ABS(THET) 
CCC --- --- --- --- --- --- --- --- --- --- --- --- --- 
CCC         REFERENCE POINT FOR FULL ITERATION          
CCC --- --- --- --- --- --- --- --- --- --- --- --- --- 
  43        CONTINUE 
C --- --- --- --- --- --- --- --- --- --- --- --- 
            IF (FLAGS) THEN 
             FLAGN=.TRUE. 
C ***************** 
CCC ---      DYNAMIC UPDATE OF THE CURRENT INTETPOLANT (in PAST)        
             DO J=1,NRDS 
                I=IPAST(J) 
                PAST(J+IACT)=Y(I)+Z3(I) 
                 Z2I=Z2(I) 
                 Z1I=Z1(I) 
                PAST(J+1*NRDS+IACT)=(Z2I-Z3(I))/C2M1 
                 AK=(Z1I-Z2I)/C1MC2 
                 ACONT3=Z1I/C1 
                 ACONT3=(AK-ACONT3)/C2 
                PAST(J+2*NRDS+IACT)=(AK-PAST(J+1*NRDS+IACT))/C1M1 
                PAST(J+3*NRDS+IACT)=PAST(J+2*NRDS+IACT)-ACONT3 
             ENDDO 
CCC          UPDATE 
             PAST(IACT)=X 
CCC          LEFT ENDPOINT OF CURRENT INTERVAL 
             PAST(IACT+IDIF-1)=H 
CCC          USED STEPSIZE                                                 
            END IF 
C --- --- --- --- --- --- --- --- --- --- --- --- 
            IF (NEWT.GE.NIT) THEN  
             INREJ=2 
             GOTO 431 
C ---------> UNEXPECTED STEP-REJECTION 
            END IF 
C ---     COMPUTE THE RIGHT-HAND SIDE 
            CONT(1:N)=Y(1:N)+Z1(1:N) 
            CALL FCN(N,X+C1*H,CONT,F1,ARGLAG,PHI,RPAR,IPAR,PAST,
     &               IPAST,NRDS) 
            CONT(1:N)=Y(1:N)+Z2(1:N) 
            CALL FCN(N,X+C2*H,CONT,F2,ARGLAG,PHI,RPAR,IPAR,PAST,
     &               IPAST,NRDS) 
            CONT(1:N)=Y(1:N)+Z3(1:N) 
            CALL FCN(N,XPH,CONT,F3,ARGLAG,PHI,RPAR,IPAR,PAST,
     &               IPAST,NRDS) 
            NFCN=NFCN+3 
C 
C --->    RHS COMPUTATION 
            IF (IMPLCT) THEN 
             ZL=0.D0 
             DO I1=1,N 
               DO J1=1,N 
                 ZL(I1)=ZL(I1)+AI11H*FMAS(I1,J1)*Z1(J1) 
                 ZL(I1)=ZL(I1)+AI12H*FMAS(I1,J1)*Z2(J1) 
                 ZL(I1)=ZL(I1)+AI13H*FMAS(I1,J1)*Z3(J1) 
 
                 ZL(I1+N)=ZL(I1+N)+AI21H*FMAS(I1,J1)*Z1(J1) 
                 ZL(I1+N)=ZL(I1+N)+AI22H*FMAS(I1,J1)*Z2(J1) 
                 ZL(I1+N)=ZL(I1+N)+AI23H*FMAS(I1,J1)*Z3(J1) 
 
                 ZL(I1+2*N)=ZL(I1+2*N)+AI31H*FMAS(I1,J1)*Z1(J1) 
                 ZL(I1+2*N)=ZL(I1+2*N)+AI32H*FMAS(I1,J1)*Z2(J1) 
                 ZL(I1+2*N)=ZL(I1+2*N)+AI33H*FMAS(I1,J1)*Z3(J1) 
               END DO 
             END DO 
            ELSE 
             ZL(1:N)=AI11H*Z1(1:N)+AI12H*Z2(1:N)+AI13H*Z3(1:N) 
             ZL(N+1:2*N)=AI21H*Z1(1:N)+AI22H*Z2(1:N)+AI23H*Z3(1:N) 
             ZL(2*N+1:3*N)=AI31H*Z1(1:N)+AI32H*Z2(1:N)+AI33H*Z3(1:N) 
            END IF 
C 
            ZL(1:N)=-ZL(1:N)-F1(1:N) 
            ZL(N+1:2*N)=-ZL(N+1:2*N)-F2(1:N) 
            ZL(2*N+1:3*N)=-ZL(2*N+1:3*N)-F3(1:N)       
C --------> SOLVE THE LINEAR SYSTEMS 
            CALL SOL(3*N,3*LDJAC,FJACL,ZL,IPJ) 
            NSOL=NSOL+1 
            NEWT=NEWT+1 
            DYNO=0.D0 
            DO I=1,N 
               DENOM=SCAL(I) 
               DYNO=DYNO+(ZL(I)/DENOM)**2+(ZL(I+N)/DENOM)**2 
     &          +(ZL(I+2*N)/DENOM)**2 
            END DO 
            DYNO=DSQRT(DYNO/N3) 
C -------------------------------------------------------- 
C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE 
C -------------------------------------------------------- 
              IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN 
                THQ=DYNO/DYNOLD 
                IF (NEWT.EQ.2) THEN 
                   THETA=THQ 
                ELSE 
                   THETA=SQRT(THQ*THQOLD) 
                END IF 
                THQOLD=THQ 
                INREJ=0 
C --- 1 
                IF (THETA.LT.0.99D0) THEN 
                    FACCON=THETA/(1.0D0-THETA) 
                    DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT 
                    IF (DYTH.GE.1.0D0) INREJ=1 
C ----------------  CONVERGENZA LENTA --- 
                ELSE 
                    INREJ=2 
C ----------------  DIVERGENZA --- 
                END IF 
C --- 1 
            ELSE 
                INREJ=0 
            END IF 
C ---------------------------------------- 
C ------------------ THE STEP IS REPEATED  
C ---------------------------------------- 
 431        IF (INREJ.GT.0) THEN 
             REJECT=.TRUE. 
             LAST=.FALSE.
             JLFLAG=0 
C -------------------------------------------- 
C ---        TURN TO A SIMPLE ITERATION         
             IF (INREJ.EQ.1) THEN 
                QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH)) 
                HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT)) 
                H=HHFAC*H 
             ELSE IF (INREJ.EQ.2) THEN 
                H=H*0.55D0 
                HHFAC=0.55D0 
             END IF 
             IF (CALJAC) GOTO 20 
             GOTO 10 
            END IF 
C -------------------------------------------------------- 
            DYNOLD=MAX(DYNO,UROUND) 
C --        UPDATE OF Z VALUES  
            Z1(1:N)=Z1(1:N)+ZL(1:N) 
            Z2(1:N)=Z2(1:N)+ZL(N+1:2*N) 
            Z3(1:N)=Z3(1:N)+ZL(2*N+1:3*N) 
C -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
            IF (FACCON*DYNO.GT.FNEWT) THEN 
             GOTO 43 
            END IF 
C - - - - - END FULL NEWTON ITERATION
C ----------------------------------------------------------------- 
C          RK EQUATIONS SUCCESFULLY SOLVED 
C ----------------------------------------------------------------- 
C *** *** *** *** *** *** *** *** *** *** *** 
C END LOOP 
C *** *** *** *** *** *** *** *** *** *** *** 
C -----------------------------------------------------------------
  55  CONTINUE   
C ******************** 
C --- ERROR ESTIMATION   
C ******************** 
C 
C    ERROR ESTIMATES                                        
       CALL ESTRAD (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, 
     &              H,U1,DD1,DD2,DD3,CL1,CL3,CQ1,CQ2,CQ3,CERLQ, 
     &              FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,E1,LDE1,ALPHA, 
     &              Z1,Z2,Z3,CONT,F1,F2,F3,IP1,IPHES,SCAL,SERR,CERR, 
     &              FIRST,REJECT,FAC1,ARGLAG,PHI,RPAR,IPAR, 
     &              IOUT,PAST,IPAST,NRDS,JLFLAG,IEFLAG) 
      FAC=MIN(SAFE,CFAC/(NEWT+2*NIT)) 

      IF (FIRST) THEN 
C ------------------------------------------------------------ 
C ---  AFTER A GRID OR BREAKING POINT            
C --------------------------------------------- 
       ERR=SERR
C ---  WE REQUIRE .2<=HNEW/H<=8. 
       QUOT=MAX(FACR,MIN(FACL,ERR**0.25D0/FAC)) 
       HNEW=H/QUOT 
      ELSE 
C   ---------------------------------------------------------- 
       IF (IEFLAG.EQ.-1) THEN 
        CERR2=H*CERR 
C ---    
        ERR=CERR2 
       ELSE IF (IEFLAG.EQ.1) THEN 
C -----    
        CERR2=H*CERR 
        ERR=CERS*SERR+CERC*CERR2 
       ELSE IF (IEFLAG.EQ.2) THEN 
C ----- STANDARD DISCRETE ERROR 
        ERR=SERR 
       ELSE IF (IEFLAG.EQ.3) THEN 
C ----- NOT AVAILABLE AT THE MOMENT 
        CERR2=H*CERR 
        ERR=CERS*SERR+CERC*CERR2 
       END IF 
C ------------------------------------ 
C ---  COMPUTATION OF HNEW 
C ---  AS PREVIOUSLY COMPUTED: 
C ------------------------------------ 
C ---  WE REQUIRE .2<=HNEW/H<=8. 
C ---  LINEAR COMBINATION OF BOTH ERROR COMPONENTS
       QUOT=MAX(FACR,MIN(FACL,ERR**.25D0/FAC)) 
 
       HNEW=H/QUOT 
      END IF
C *** *** *** *** *** *** *** 
C  DOES THE ERROR INCREASE MUCH? 
C *** *** *** *** *** *** *** 
      REPEAT=.FALSE.
      IF (ERR.GE.1.D0.OR.((ERR/ERRACC.GE.TCKBP).AND.(.NOT.BPD))) THEN 
C ---   KIND OF STEPSIZE
        IF (BPD) THEN
C ---     BP IS WRONG! REPEAT
          IBP=IBP-1
          BPD=.FALSE.
          REPEAT=.TRUE.
          IF (ERR.LT.1.D0) HNEW=H*0.55D0
        ELSE
C ---     LOOK FOR A BP        
          HP=H*0.99D0
          CALL BPDTCT(N,X,HP,Y,ARGLAG,RPAR,IPAR,UCONT,GRID,NLAGS, 
     &                FIRST,LAST,XEND,IGRID,BPV,IBP,ILBP,BPP,BPD,
     &                KMAX,PHI,PAST,IPAST,NRDS) 
          IF (BPD) REPEAT=.TRUE.
        END IF
      END IF 
      BPC=.FALSE. 
C *** *** *** *** *** *** *** 
C  IS THE ERROR SMALL ENOUGH ? 
C *** *** *** *** *** *** *** 
C     IF (ERR.LT.1.D0) THEN 
 56   IF (.NOT.REPEAT.AND.ERR.LT.1.D0) THEN
C -------------------- 
C --- STEP IS ACCEPTED   
C -------------------- 
         FLAGN=.FALSE. 
         CALJACL=.FALSE. 
C ---    UPDATE NUMBER OF ACCEPTED STEPS              
         NACCPT=NACCPT+1 
         IF (PRED) THEN 
C --------> PREDICTIVE CONTROLLER OF GUSTAFSSON 
            IF (NACCPT.GT.1) THEN 
             IF (FLAGUS) THEN 
               FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE 
               FACGUS=MAX(FACR,MIN(FACL,FACGUS)) 
               QUOT=MAX(QUOT,FACGUS) 
               HNEW=H/QUOT 
             ELSE 
               FLAGUS=.TRUE. 
             END IF 
            END IF 
            HACC=H 
         END IF 
         ERRACC=MAX(1.0D-2,ERR)
C        ERRACC=ERR 
         XOLD=X 
         HOLD=H 
         X=XPH  
C ---    AGGIORNAMENTO DELLA SOLUZIONE 
         DO I=1,N 
            Z3I=Z3(I) 
            YI=Y(I)+Z3I   
            Y(I)=YI 
C ------------------------------------------------------------------------- 
C --------- UPDATE THE STAGES AND SOLUTION 
C ------------------------------------------------------------------------- 
C --- 
              CONT(I)=YI 
              Z2I=Z2(I) 
              Z1I=Z1(I) 
C ------------------------------------------------------------------- 
C ---         INVERSE DEVIDED DIFFERENCES
C ------------------------------------------------------------------- 
              A1=(Z2I-Z3I)/C2M1 
              CONT(I+N)=A1 
              AK=(Z1I-Z2I)/C1MC2 
              ACONT3=Z1I/C1 
              ACONT3=(AK-ACONT3)/C2 
              A2=(AK-CONT(I+N))/C1M1 
              CONT(I+N2)=A2 
            IF (.NOT.FIRST) THEN 
              CONT(I+N3)=A2-ACONT3 
            ELSE 
C ---         QUADRATIC APPROXIMATION
              CONT(I+N3)=0.D0			
C ---         INVECE DI:  
C ---         CONT(I+N3)=CONT(I+N2)-ACONT3 
            END IF 
         END DO 
C ---------------------------------------- 
C ---    SAVE LAST ACCEPTED STEP INFORMATION   
         DO I=1,LRC 
          UCONT(I)=CONT(I) 
         END DO 
         UCONT(LRC+1)=X 
         UCONT(LRC+2)=H 
C ---    FOR POSSIBLE SEARCH OF BREAKING POINTS     
C ------------------------------------------------------------------ 
C ---    STEP IS ACCEPTED> DENSE OUTPUT IS STORED IN PAST        
C ---    
         DO J=1,NRDS 
            I=IPAST(J) 
            PAST(J+IACT)=CONT(I) 
            PAST(J+1*NRDS+IACT)=CONT(I+N) 
            PAST(J+2*NRDS+IACT)=CONT(I+N2) 
            IF (.NOT.FIRST) THEN 
              PAST(J+3*NRDS+IACT)=CONT(I+N3) 
            ELSE 
C ---       QUADRATIC APPROXIMATION                        
              PAST(J+3*NRDS+IACT)=0.D0 
            END IF 
         ENDDO 
C ---------> AGGIORNAMENTO DI PAST <--------- 
         PAST(IACT)=XOLD 
C ---    
         IACT=IACT+IDIF 
C ---    POINTER TO NEXT STEP 
         PAST(IACT-1)=H 
C ---    
         IF (IACT+IDIF-1.GT.MXST*IDIF) IACT=1  
C ---    CONTROL ON THE MEMORY DIMENSION
C ----------------------------------------------------------------- 
C ---    COMPUTATION AT BP FOR NEXT STEP
         IF (LAST) THEN 
C ---     LAST HAS TO BE RE/DEFINED
          IF (BPD) THEN 
           IGRID=IGRID-1
	  END IF 
C --- 
          IF (.NOT.IMPLCT.OR.NEUTRAL) THEN  ! EXPLICIT PROBLEM
            HE=DMAX1(H/1.D4,10.D0*UROUND) 
C -------------------- 
            LEFT=.TRUE. 
C -------------------- 
C --- 
C ---       EULER STEP
            CALL FCN(N,X,Y,F2,ARGLAG,PHI,RPAR,IPAR,PAST,IPAST,NRDS) 
            IF (NEUTRAL) THEN
	        DO I=1,N-NDIMN 
               Z2(I)=Y(I)+HE*F2(I) 
              END DO
		    DO I=1,NDIMN
	         Z2(N-NDIMN+I)=F2(IPAST(NRDS+I))+HE*F2(N-NDIMN+I)
	        END DO
            ELSE
		    DO I=1,N 
               Z2(I)=Y(I)+HE*F2(I) 
              END DO
		  END IF  
C -------------------- 
            LEFT=.FALSE. 
C -------------------- 
            CALL FCN(N,X,Y,F3,ARGLAG,PHI,RPAR,IPAR,PAST,IPAST,NRDS) 
            IF (NEUTRAL) THEN
	        DO I=1,N-NDIMN 
               Z3(I)=Y(I)+HE*F3(I) 
              END DO
		    DO I=1,NDIMN
	         Z3(N-NDIMN+I)=F3(IPAST(NRDS+I))+HE*F3(N-NDIMN+I)
	        END DO
            ELSE
		    DO I=1,N 
               Z3(I)=Y(I)+HE*F3(I) 
              END DO 
            END IF
            XL =ARGLAG(ILBP,X+HE,N,Z2,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
            XLR=ARGLAG(ILBP,X+HE,N,Z3,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
            IF (XL.GE.BPP.AND.XLR.GE.BPP) THEN 
             LEFT=.FALSE. 
            ELSE IF (XL.LT.BPP.AND.XLR.LT.BPP) THEN 
             LEFT=.TRUE. 
            ELSE 
             IF (IOUT.EQ.1) THEN 
		    IF (XL.GT.BPP) THEN
			 WRITE (6,*) 
     &        ' WARNING!: SOLUTION DOES NOT EXIST AT X= ',X
              ELSE
	         WRITE (6,*) 
     &        ' WARNING!: SOLUTION IS  NOT UNIQUE AT X= ',X
	        END IF
              GO TO 980
C             RETURN 
             END IF
            END IF
C ---       PROJECTION FOR DERIVATIVE COMPONENTS OF NEUTRAL EXPLICIT PROBLEMS
            PROJECT=.TRUE.
	    IF (NEUTRAL.AND.PROJECT) THEN
	       IF (LEFT) THEN
		    DO J=1,NDIMN
	         Y(N-NDIMN+J)=F2(IPAST(NRDS+J))
	        END DO
            ELSE
	        DO J=1,NDIMN
	         Y(N-NDIMN+J)=F3(IPAST(NRDS+J))
	        END DO
	       END IF
            END IF
	  ELSE ! GENERAL IMPLICIT
	      LEFT=.TRUE.
	  END IF
C --- 
          BPD=.FALSE. 
         END IF  
C ----------------------------------------------------------------- 
         IF (ITOL.EQ.0) THEN 
             DO I=1,N 
                SCAL(I)=ATOL(1)+RTOL(1)*ABS(Y(I)) 
             END DO 
         ELSE 
             DO I=1,N 
                SCAL(I)=ATOL(I)+RTOL(I)*ABS(Y(I)) 
             END DO 
         END IF 
         IF (IOUT.NE.0) THEN 
C --- 
             NRSOL=NACCPT+1 
             XSOL=X 
             XOSOL=XOLD 
             NSOLU=N 
             HSOL=HOLD 
C --- 
             CALL SOLOUT(NRSOL,XOSOL,XSOL,HSOL,Y,CONT,LRC,NSOLU, 
     &                   RPAR,IPAR,IRTRN)
             IF (IRTRN.LT.0) GOTO 179 
         END IF 
         CALJAC=.FALSE. 
C ----------------------------- 
C ---    FIRST IS RESTORED 
         FIRST=.FALSE. 
         IF (LAST) FIRST=.TRUE. 
C ---    COMPUTATION OF Y0 
          POSNEG=1.0D0
          IF (LEFT) POSNEG=-POSNEG
          XEPS = POSNEG*(ABS(X)+1.0D0)*100*UROUND
          CALL FCN(N,X+XEPS,Y,Y0,ARGLAG,PHI,RPAR,IPAR,PAST,IPAST,NRDS)
          IF (FIRST) THEN
	        DO I=1,NDIMN
	          Y0(N-NDIMN+I)=0.0D0
	        END DO
	      END IF
         NFCN=NFCN+1 
         FIRST=.FALSE.  
C -----------------------------  
C ------ FINAL POINT                               
         IF (LAST) THEN 
 45         CONTINUE
            IF (IGRID.GE.NGRID) THEN 
               H=HOPT 
               IDID=1 
C ---          END OF COMPUTATION 
               GOTO 980 
            ELSE 
               IGRID=IGRID+1 
               LAST=.FALSE.            
C              LEFT=.FALSE. 
               FIRST=.TRUE. 
               XEND=GRID(IGRID)
	       IF (ABS(XEND-X).LE.(H*1.D-2)) THEN
		  IGRID=IGRID+1
		  GO TO 45 
	       END IF 
               FLAGUS=.FALSE. 
               IF (WORK7.EQ.0.D0) THEN 
                  HMAXN=XEND-X  
                  HMAX=HMAXN 
               END IF 
               HNEW=MIN(HNEW,H)
            END IF 
         END IF 
         HNEW=MIN(HNEW,HMAXN) 
         HOPT=MIN(H,HNEW) 
         IF (REJECT) HNEW=MIN(HNEW,H)  
         REJECT=.FALSE. 
         IF ((X+HNEW/QUOT1-XEND).GE.0.D0) THEN 
            H=XEND-X 
            IF (H.LT.0.D0) THEN 
             WRITE(6,*) 'ERROR!: NEGATIVE STEPSIZE! AT ' 
             WRITE(6,*) 'X > XEND = ',X,XEND 
             STOP 
            END IF 
            LAST=.TRUE. 
         ELSE 
C ----------------------------------------------------- 
C --------  IN ORDER TO AVOID VERY SMALL END-STEPSIZES:    
            IF ((X+1.8D0*HNEW-XEND).GT.0.D0) THEN  
              H=(XEND-X)*0.55D0 
            ELSE 
              QT=HNEW/H  
              HHFAC=H 
              IMANT=0 
C ---         STEP IS MAINTAINED     
              IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) THEN 
               IF ((JLFLAG.EQ.0).AND.(.NOT.FIRST)) IMANT=1 
              ELSE  
               H=HNEW  
              END IF 
            END IF 
         END IF 
         HHFAC=H 
C ------------------------------------------- 
C ---    SIMPLE ITERATION FIRST                 
         JLFLAG=0 
C ---    
         IF ((.NOT.FIRST).AND.(THETA.LE.THET)) GOTO 20 
         GOTO 10 
C -------------------- 
C --- STEP IS ACCEPTED   
C -------------------- 
      ELSE 
C -------------------- 
C --- STEP IS REJECTED   
C -------------------- 
         IF (BPD) THEN  
          IF (.NOT.FIRST) THEN 
           HNEW=HP 
           LAST=.TRUE. 
          END IF 
         ELSE 
          LAST=.FALSE. 
         END IF 
         FLAGN=.FALSE. 
C ------------------------------- 
         IF (IRTRN.LT.0) GOTO 179 
         REJECT=.TRUE. 
C --- 
         IF (FIRST) THEN 
             H=H*0.12D0 
             HHFAC=0.12D0 
         ELSE  
             HHFAC=HNEW/H 
             H=HNEW  
         END IF 
         IF (NACCPT.GE.1) NREJCT=NREJCT+1 
C ---> 
         JLFLAG=0  
         IF (CALJAC) GOTO 20 
         GOTO 10 
C -------------------- 
C --- STEP IS REJECTED   
C -------------------- 
      END IF 
C ---------------------------------- 
C --- END OF ERROR CONTROL PROCEDURE 
C ---------------------------------- 
C --- UNEXPECTED STEP-REJECTION 
  78  CONTINUE 
      REJECT=.TRUE. 
      FLAGN=.FALSE. 
C --- 
      LAST=.FALSE. 
C --- 
      IF (IER.NE.0) THEN 
          NSING=NSING+1 
          IF (NSING.GE.5) GOTO 176 
      END IF 
C --- 
      H=H*0.55D0  
      HHFAC=0.55D0 
C --- 
      IF (CALJAC) GOTO 20 
      IF (IRTRN.LT.0) GOTO 175 
      GOTO 10 
C --- FAIL EXIT 
 175  CONTINUE 
      IDID=-5 
      GOTO 980 
 176  CONTINUE 
      WRITE(6,979)X    
      WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER 
      IDID=-4 
      GOTO 980 
 177  CONTINUE 
      WRITE(6,979)X    
      WRITE(6,*) ' STEP SIZE TOO SMALL, H=',H 
      IDID=-3 
      GOTO 980 
 178  CONTINUE 
      WRITE(6,979)X    
      WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'  
      IDID=-2 
      GOTO 980 
C --- EXIT CAUSED BY SOLOUT 
 179  CONTINUE 
      WRITE(6,979)X 
 979  FORMAT(' EXIT OF RADAR5 AT X=',E18.4)  
      IDID=2 
C --- RETURN LABEL 
 980  CONTINUE 
      WRITE(6,*) IBP,' COMPUTED BREAKING POINTS: '
      WRITE(8,*) 'BREAKING POINTS: '
      DO I=1,IBP
	    WRITE(8,*) BPV(I)
      END DO
      WRITE(8,*) ' -------------- '
      CLOSE(8)
C --- DEALLOCATION OF THE MEMORY 
      DEALLOCATE (Z1,Z2,Z3,Y0,SCAL,F1,F2,F3) 
      DEALLOCATE(BPV) 
      DEALLOCATE(FJAC,ZL) 
      IF (IMPLCT) DEALLOCATE(FMAS) 
      DEALLOCATE (IP1,IP2,IPHES) 
      DEALLOCATE (E1,E2R,E2I) 
      DEALLOCATE (PAST) 
      IF (NLAGS.GT.0) THEN 
       DEALLOCATE (FJACS,FJACLAG) 
       DEALLOCATE (IVL,IVE,IVC,ILS,ICOUN) 
       IF (ISWJL.NE.1) DEALLOCATE (IPJ,FJACL,XLAG) 
      END IF 
      DEALLOCATE (CONT,UCONT) 
      RETURN 
      END 
C 
C     END OF SUBROUTINE RADCOR 
C 
C *********************************************************** 
C 
      SUBROUTINE LAGR5(IL,X,N,Y,ARGLAG,PAST,THETA,IPOS,RPAR,IPAR,
     &                 PHI,IPAST,NRDS)
C ---------------------------------------------------------- 
C     THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT IN CONNECTION 
C     WITH THE OUTPUT-SUBROUTINE FOR RADAR5. IT PROVIDES THE 
C     POSITION OF THE DENSE OUTPUT AT THE IL-TH DELAY. 
C ---------------------------------------------------------- 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      INTEGER, PARAMETER :: DP=kind(1D0) 
      REAL(kind=DP), dimension(1), intent(in) :: Y,PAST 
      REAL(kind=DP), dimension(1), intent(in) :: RPAR 
      INTEGER, dimension(1), intent(in) :: IPAR 
      INTEGER, dimension(NRDS+2*N), intent(inout) :: IPAST 
      EXTERNAL PHI,ARGLAG
      LOGICAL FLAGS,FLAGN,FIRST,LAST,REJECT,BPD,LEFT 
C --- COMMON BLOCKS 
      COMMON /BPLOG/FIRST,LAST,REJECT,BPD 
      COMMON /BPCOM/BPP,ILBP,LEFT 
      COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
C --- COMPUTE DEVIATED ARGUMENT FOR IL-TH DELAY 
      XLAG=ARGLAG(IL,X,N,Y,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
C --- INITIAL PHASE 
      THETA=XLAG 
      IPOS=-1 
C     EPSILON GIVES THE OVERLAPPING
C     MIN VALUE FOR THE SUPER-POSITION NEIGHBOURHOOD OF A BP
      COMPAR=UROUND*MAX(ABS(XLAG),ABS(X0B))
      EPSACT=10.D0*COMPAR
	  IF (IACT.GT.1) EPSACT=DMAX1(PAST(IACT-1)*1.D-2,EPSACT)
      IF (XLAG.LE.X0B) THEN 
C ---     DEVIATING ARGUMENT ON THE INITIAL SEGMENT         
        IF (.NOT.((IL.EQ.ILBP).AND.(BPD.OR.FIRST))) THEN 
	      IF (XLAG-X0B.LT.0.D0) THEN
	        RETURN
	      ELSE
	        IPOS=1
	        THETA=-1.D0 
	      END IF
	    ELSE
          IF (ABS(XLAG-X0B).LE.EPSACT) THEN
	        IF (LEFT) THEN
                IPOS=-1 
                THETA=XLAG
            ELSE
                IPOS=1
	            THETA=(XLAG-(PAST(IPOS)+PAST(IPOS+IDIF-1)))/
     1     	    PAST(IPOS+IDIF-1) 
            END IF
          ELSE IF (ABS(XLAG-BPP).LE.EPSACT) THEN 
	        IPOS=-1 
            IF (LEFT) THEN 
              IF (XLAG.GT.BPP) THEN 
                IF (BPP.GT.0.D0) THEN 
                  THETA=BPP*(1.D0-100*UROUND) 
                ELSE 
                  THETA=BPP*(1.D0+100*UROUND) 
                END IF 
              END IF 
            ELSE 
              IF (XLAG.LT.BPP) THEN 
                IF (BPP.GT.0) THEN 
                  THETA=BPP*(1.D0+100*UROUND) 
                ELSE 
                  THETA=BPP*(1.D0-100*UROUND) 
                END IF 
              END IF 
            END IF 
          END IF
	    END IF    
	    RETURN
	  END IF 
C --- COMPUTE THE POSITION OF XLAG 
      IPA = IACT+IDIF 
      IF (IPA.GT.(MXST-1)*IDIF+1) IPA=1 
      IF (XLAG-PAST(IPA).LT.0.D0) THEN
         WRITE (6,*) ' MEMORY FULL, MXST = ',MXST 
         IRTRN=-1
	     STOP 
      END IF 
      INEXT=IACT-IDIF 
      IF (INEXT.LT.1) INEXT=(MXST-1)*IDIF+1 
      XRIGHT=PAST(INEXT)+PAST(INEXT+IDIF-1) 
C -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
C --- INITIALIZE POSITION INSIDE THE MEMORY 
      IPOS=IPAST(IL+2*N) 
C --- STEP AND DELAYS 
      IF (XLAG-XRIGHT.GT.0.D0) THEN 
        IF (.NOT.FLAGN) THEN 
           IPOS=IACT-IDIF 
           IF (IPOS.LT.1) IPOS=(MXST-1)*IDIF+1 
        ELSE 
           IPOS=IACT 
        END IF 
        FLAGS=.TRUE. 
C ------------------------ 
C ----- COMPUTE THETA (>0)     
C ------------------------ 
        H=PAST(IPOS+IDIF-1) 
        THETA=(XLAG-(PAST(IPOS)+H))/H 
C -----    EXTRAPOLATION USE OF THE COLLOCATION POLYNOMIAL            
C ----- 
        RETURN                       
      ELSE 
   1    CONTINUE 
        IF (XLAG-PAST(IPOS).LE.0.D0) THEN 
           IPOS=IPOS-IDIF 
           IF (IPOS.LT.1) IPOS=(MXST-1)*IDIF+1 
           GOTO 1 
        END IF 
   2    CONTINUE 
        INEXT=IPOS+IDIF 
        IF (INEXT.GT.(MXST-1)*IDIF+1) INEXT=1 
        IF (XLAG.GT.PAST(INEXT).AND.INEXT.NE.IACT) THEN 
           IPOS=INEXT 
           GOTO 2 
        END IF
	  IF (IPOS.EQ.1) THEN 
            IPREV=(MXST-1)*IDIF+1 
        ELSE 
            IPREV=IPOS-IDIF 
        END IF
C --- 
C ---   IN CORRESPONDENCE OF BREAKING POINTS                
C --- 
 3      CONTINUE 
        IF (.NOT.((IL.EQ.ILBP).AND.(BPD.OR.FIRST))) GOTO 10
        IF (BPP.EQ.X0B) THEN
	      IF (LEFT) THEN
	        IPOS=-1
	        THETA=XLAG
	        RETURN
          ELSE
	        IF (IPOS.EQ.-1) IPOS=1
	        GO TO 10
	      END IF
        END IF
	    IPOSB=0  
	    IF (ABS(BPP-PAST(IPOS)).LE.10.D0*UROUND) THEN
            IPOSB=IPOS
	    ELSE IF (ABS(BPP-PAST(INEXT)).LE.10.D0*UROUND) THEN
	        IPOSB=INEXT
	      END IF
          IF (IPOSB.EQ.0) THEN
	        GO TO 10
	      END IF
           IF (IPOSB.EQ.1) THEN
	         EPSILON=(PAST(IPOSB+IDIF)-PAST(IPOSB))
	       ELSE IF (IPOSB.EQ.(MXST-1)*IDIF+1) THEN
	           EPSILON=(PAST(IPOSB)-PAST(IPOSB-IDIF))
	         ELSE
	           EPSILON=DMIN1(PAST(IPOSB+IDIF)-PAST(IPOSB),
     1                  PAST(IPOSB)-PAST(IPOSB-IDIF))
	       END IF
	       EPSILON=DMAX1(EPSILON*1.D-2,EPSACT)
           IF (ABS(XLAG-BPP).GT.EPSILON) GOTO 10 
           IF (IPOSB.EQ.1) THEN
	       IF (LEFT) THEN 
                IPOS=-1 
                THETA=XLAG
                RETURN
	       ELSE 
                IPOS=1 
	            GO TO 10
	       END IF
        END IF
        IF (LEFT) THEN           
C ---     PREVIOUS INTERVAL HAS TO BE SELECTED
          IPOS=IPOSB-IDIF
        ELSE 
C ---     NEXT INTERVAL HAS TO BE SELECTED
          IPOS=IPOSB
        END IF    
C ----- COMPUTE THETA (<0): SITUAZIONE PIU' TIPICA     
  10    THETA=(XLAG-(PAST(IPOS)+PAST(IPOS+IDIF-1)))/PAST(IPOS+IDIF-1) 
C ----- REM: THETA IS NEGATIVE                                  
	  END IF 
C --- UPDATE POSITION INSIDE THE MEMORY 
      IPAST(IL+2*N)=IPOS 
      RETURN 
      END 
C 
C     END OF FUNCTION LAGR5 
C 
C -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 

C *********************************************************** 
C 
      FUNCTION YLAGR5(IC,THETA,IPOS,PHI,RPAR,IPAR,PAST,IPAST,NRDS) 
C ---------------------------------------------------------- 
C     THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT IN CONNECTION 
C     WITH THE SUBROUTINE LAGR5. IT PROVIDES AN APPROXIMATION  
C     TO THE IC-TH COMPONENT OF THE SOLUTION AT XLAG. 
C ---------------------------------------------------------- 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      INTEGER, PARAMETER :: DP=kind(1D0) 
      REAL(kind=DP), dimension(1), intent(in) :: RPAR 
      INTEGER, dimension(1), intent(in) :: IPAR 
      REAL(kind=DP), dimension(1), intent(in) :: PAST 
      INTEGER, dimension(1), intent(in) :: IPAST 
      LOGICAL FLAGS,FLAGN 
      EXTERNAL PHI
C---- COMMON BLOCKS 
      COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2 
      COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
C 
C --- INITIAL PHASE 
      IF (IPOS.EQ.-1) THEN 
            YLAGR5=PHI(IC,THETA,RPAR,IPAR)   
C ---       CALL PHI(IC,THETA,YLAGR5,RPAR,IPAR)
            RETURN 
      END IF 
C --- 
C --- COMPUTE PLACE OF IC-TH COMPONENT  
      I=0  
      DO 5 J=1,NRDS  
      IF (IPAST(J).EQ.IC) I=J 
   5  CONTINUE 
      IF (I.EQ.0) THEN 
         WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',IC 
         RETURN 
      END IF   
C ----- COMPUTE DESIRED APPROXIMATION 
      I=I+IPOS 
      YLAGR5=PAST(I)+THETA*(PAST(NRDS+I)+(THETA-C2M1)*(PAST(2*NRDS+I) 
     &            +(THETA-C1M1)*(PAST(3*NRDS+I)))) 
      RETURN 
      END 
C 
C     END OF FUNCTION YLAGR5 
C 
 
C *********************************************************** 
C 
      FUNCTION DLAGR5(IC,THETA,IPOS,PHI,RPAR,IPAR,PAST,IPAST,NRDS) 
C ---------------------------------------------------------- 
C     THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT IN CONNECTION 
C     WITH THE SUBROUTINE LAGR5. IT PROVIDES AN APPROXIMATION  
C     TO THE IC-TH COMPONENT OF THE SOLUTION DERIVATIVE AT XLAG. 
C ---------------------------------------------------------- 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      INTEGER, PARAMETER :: DP=kind(1D0) 
      REAL(kind=DP), dimension(1), intent(in) :: RPAR 
      INTEGER, dimension(1), intent(in) :: IPAR 
      REAL(kind=DP), dimension(1), intent(in) :: PAST 
      INTEGER, dimension(1), intent(in) :: IPAST 
      LOGICAL FLAGS,FLAGN 
      EXTERNAL PHI
C---- COMMON BLOCKS 
      COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2 
      COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
C 
C --- INITIAL PHASE 
      IF (IPOS.EQ.-1) THEN 
C           DLAGR5=DPHI(IC,THETA,RPAR,IPAR)   
            DLAGR5=0.D0 
            RETURN 
      END IF 
C --- 
C --- COMPUTE PLACE OF IC-TH COMPONENT  
      I=0  
      DO 5 J=1,NRDS  
      IF (IPAST(J).EQ.IC) I=J 
   5  CONTINUE 
      IF (I.EQ.0) THEN 
         WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',IC 
         RETURN 
      END IF   
C ----- COMPUTE DESIRED APPROXIMATION 
      H=PAST(IPOS+IDIF-1) 
      I=I+IPOS 
      DLAGR5=PAST(NRDS+I)+(THETA-C2M1)*(PAST(2*NRDS+I) 
     &     +(THETA-C1M1)*PAST(3*NRDS+I)) 
     &     +THETA*(PAST(2*NRDS+I)+(2.D0*THETA-C2M1-C1M1)*PAST(3*NRDS+I)) 
      DLAGR5=DLAGR5/H 
      RETURN 
      END 
C 
C     END OF FUNCTION DLAGR5 
C 
 
C 
      SUBROUTINE BPDTCT(N,X,H,Y,ARGLAG,RPAR,IPAR,UCONT,GRID,NLAGS, 
     &                  FIRST,LAST,XEND,IGRID,BPV,IBP,ILBP,BPP,BPD,
     &                  KMAX,PHI,PAST,IPAST,NRDS)                   
C ---------------------------------------------------------- 
C     THIS SUBROUTINE CAN BE USED FOR DETECTING BREAKING POINTS  
C     WITH THE OUTPUT-SUBROUTINE FOR RADAR5. IT PROVIDES AN 
C     APPROXIMATION TO THE IC-TH COMPONENT OF THE SOLUTION AT X. 
C ---------------------------------------------------------- 
        IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
        INTEGER, PARAMETER :: DP=kind(1D0) 
        REAL(kind=DP), dimension(N) :: Y 
        REAL(kind=DP), dimension(1), intent(in) :: PAST 
        INTEGER, dimension(1), intent(in) :: IPAST 
        REAL(kind=DP), dimension(:), allocatable :: YADV 
        REAL(kind=DP), dimension(1) :: UCONT 
        REAL(kind=DP), dimension(1), intent(inout) :: GRID 
        REAL(kind=DP), dimension(1) :: BPV 
        INTEGER, dimension(1) :: IPAR 
        REAL(kind=DP), dimension(1) :: RPAR 
        EXTERNAL PHI
        LOGICAL FIRST,LAST,BPD,FLAGS,FLAGN 
C----   COMMON BLOCKS 
        COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2 
        COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
C ---     
        IF (FIRST) RETURN
        BPD=.FALSE.  
        LRC=4*N 
        EPSILON=1.D-10
        ALLOCATE(YADV(N)) 
        COMPAR=UROUND*MAX(ABS(X),ABS(X+H)) 
	    XLAST=UCONT(LRC+1) 
        HLAST=UCONT(LRC+2) 
        DO IL=1,NLAGS 
         ALS = ARGLAG(IL,X,N,Y,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
C -----  DEVIATING ARGUMENT AT X 
C -----  EXTRAPOLATION OF THE COLLOCATION POLYNOMIAL       
         DO IC=1,N 
          YADV(IC)=CONTR5(IC,N,X+H,UCONT,XLAST,HLAST) 
         END DO 
         ALD = ARGLAG(IL,X+H,N,YADV,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
C -----  DEVIATING ARGUMENT AT X+H        
         IF (ABS(ALS-ALD).LE.COMPAR) GO TO 33
         DO L=1,IGRID-1 
          BPP=GRID(L) 
          IF ((ALS-BPP)*(ALD-BPP).LT.COMPAR) THEN 
           BPD=.TRUE. 
C          BREAKING POINT! 
           GO TO 33 
          END IF 
         END DO    
         DO L=IBP,1,-1 
          BPP=BPV(L) 
          IF ((ALS-BPP)*(ALD-BPP).LT.COMPAR) THEN 
C          BREAKING POINT! 
           BPD=.TRUE. 
           GO TO 33 
          END IF 
         END DO    
 33      CONTINUE 
         IF (BPD) THEN
C --- 
           THLIM=1.D0 
           THRIGH=THLIM 
C -------------------------------------------------- 
           THLEFT=0.D0
C ---      
           DO K=1,KMAX   
            THNEW = THLEFT - (ALS-BPP)*(THRIGH-THLEFT)/(ALD-ALS)  
C ---       TEST DI CONVERGENZA 
            IF (ABS(THRIGH-THNEW).LE.EPSILON.OR. 
     &          ABS(THLEFT-THNEW).LE.EPSILON) GOTO 36 
            XA=X+THNEW*H 
            DO IC=1,N 
             YADV(IC)=CONTR5(IC,N,XA,UCONT,XLAST,HLAST) 
            END DO 
            ALN = ARGLAG(IL,XA,N,YADV,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
            IF ((ALS-BPP)*(ALN-BPP).LE.0.D0) THEN 
             ALD=ALN 
             THRIGH=THNEW 
            ELSE 
             ALS=ALN 
             THLEFT=THNEW 
            END IF 
C --- 
           END DO 
 36        CONTINUE 
C ---      BP FOUND!      
           IF ((THNEW.GT.COMPAR).AND.(THNEW.LT.THLIM)) THEN 
            HP=THNEW*H 
C ---
            IF (HP.LE.1.D2*COMPAR) THEN 
             BPD=.FALSE. 
             GO TO 37 
            END IF 
C ---
          DO L=1,IGRID-1
            BPP=GRID(L) 
            IF (ABS(BPP-X-HP).LE.COMPAR) THEN
              BPD=.FALSE. 
              GOTO 37
            END IF
          END DO 
            IBP=IBP+1 
            BPV(IBP)=X+HP 
            H=HP 
            ILBP=IL 
            GOTO 37 
           ELSE 
C ---       BP ALREADY PRESENT
            BPD=.FALSE. 
           END IF  
         END IF 
        END DO 
 37     CONTINUE 
        DEALLOCATE(YADV) 
        RETURN 
        END 
 

      SUBROUTINE BPACC(N,X,H,Y,ARGLAG,RPAR,IPAR,Z1,Z2,Z3,FIRST,
     &                 BPV,IBP,ILBP,BPP,KMAX,PHI,PAST,IPAST,NRDS)                     
C ---------------------------------------------------------- 
C     THIS SUBROUTINE CAN BE USED FOR APPROXIMATING BREAKING POINTS  
C     IN TANDEM WITH THE SIMPLIFIED NEWTON ITERATION.. 
C ---------------------------------------------------------- 
        IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
        INTEGER, PARAMETER :: DP=kind(1D0) 
        REAL(kind=DP), dimension(N) :: Y,Z1,Z2,Z3 
        REAL(kind=DP), dimension(1), intent(in) :: PAST
        INTEGER, dimension(1), intent(in) :: IPAST
        REAL(kind=DP), dimension(:), allocatable :: YCONT,YAPP
        REAL(kind=DP), dimension(1) :: BPV 
        INTEGER, dimension(1) :: IPAR 
        REAL(kind=DP), dimension(1) :: RPAR 
        EXTERNAL PHI
        LOGICAL FIRST,FLAGS,FLAGN
C----   COMMON BLOCKS 
        COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2 
        COMMON /POSITS/X0B,UROUND,HMAX,IACT,IRTRN,IDIF,MXST,FLAGS,FLAGN 
C ---     
        ALLOCATE(YCONT(4*N),YAPP(N))
        EPSILON=UROUND*1.D3
C ---   DYNAMIC UPDATE
        DO I=1,N 
            Z3I=Z3(I) 
            YI=Y(I)+Z3I   
            YCONT(I)=YI 
            Z2I=Z2(I) 
            Z1I=Z1(I) 
            A1=(Z2I-Z3I)/C2M1 
            YCONT(I+N)=A1 
            AK=(Z1I-Z2I)/C1MC2 
            ACONT3=Z1I/C1 
            ACONT3=(AK-ACONT3)/C2 
            A2=(AK-YCONT(I+N))/C1M1 
            YCONT(I+2*N)=A2 
            IF (.NOT.FIRST) THEN 
              YCONT(I+3*N)=A2-ACONT3 
            ELSE 
              YCONT(I+3*N)=0.D0 
            END IF 
        END DO 
C ---   INITIAL VALUES FOR THE COMPUTATION
        XSOL=X+H
        HSOL=H
        THLEFT=0.9D0
        THRIGH=1.0D0
        XL=X+THLEFT*H
        DO I=1,N 
            YAPP(I)=CONTR5(I,N,XL,YCONT,XSOL,HSOL) 
        END DO 
        ALS = ARGLAG(ILBP,XL,N,YAPP,RPAR,IPAR,PHI,PAST,IPAST,NRDS)
C ---
        XR=X+THRIGH*H
        ALD = ARGLAG(ILBP,XR,N,YCONT,RPAR,IPAR,PHI,PAST,IPAST,NRDS)
        DO K=1,KMAX 
            IF (ABS(ALD-ALS).LE.EPSILON) THEN
              FAC=0.0D0
            ELSE
              FAC = (ALD-BPP)/(ALD-ALS)
            END IF
            THNEW = THRIGH - FAC*(THRIGH-THLEFT)
            THLEFT= THRIGH
C ---       TEST DI CONVERGENZA 
            IF (ABS(THNEW-THRIGH).LE.EPSILON) GOTO 36 
            IF ((THNEW.LE.0.5D0).OR.(THNEW.GE.1.5D0)) THEN 
              DEALLOCATE(YCONT,YAPP)
              RETURN
            END IF
            THRIGH=THNEW
            XAP=X+THRIGH*H 
            ALS=ALD
            DO I=1,N 
             YAPP(I)=CONTR5(I,N,XAP,YCONT,XSOL,HSOL) 
            END DO 
            ALD = ARGLAG(ILBP,XAP,N,YAPP,RPAR,IPAR,PHI,PAST,IPAST,NRDS) 
            IF (ABS(ALD-ALS).LE.EPSILON) GOTO 36 
        END DO 
 36     CONTINUE
C ---   BP FOUND 
        H=MIN(THRIGH,THLEFT)*H
        BPV(IBP)=X+H 
        DEALLOCATE(YCONT,YAPP)
        RETURN 
        END 
*********************************************************
C
      DOUBLE PRECISION FUNCTION CONTR5(I,N,X,CONT,XSOL,HSOL) 
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
C     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
C     THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAR5).
C ----------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER, PARAMETER :: DP=kind(1D0)
      REAL(kind=DP), dimension(1), intent(in) ::  CONT
C --- REQUIRED CONSTANTS
      COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2
      N2=2*N
      N3=3*N
      S=(X-XSOL)/HSOL
      CONTR5=CONT(I)+S*(CONT(I+N)+(S-C2M1)*(CONT(I+N2)
     &     +(S-C1M1)*CONT(I+N3)))
      RETURN
      END
C
C     END OF FUNCTION CONTR5
C
C ***********************************************************
C
      DOUBLE PRECISION FUNCTION DONTR5(I,N,X,CONT,XSOL,HSOL) 
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION DERIVATIVE AT .
C     X. IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
C     THE LAST SUCCESSFULLY COMPUTED STEP (BY DELAY5).
C ----------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER, PARAMETER :: DP=kind(1D0)
      REAL(kind=DP), dimension(1), intent(in) ::  CONT
C --- REQUIRED CONSTANTS
      COMMON /CONSTN/C1,C2,C1M1,C2M1,C1MC2
      N2=2*N
      N3=3*N
      S=(X-XSOL)/HSOL
	  DONTR5=(CONT(I+N)+
     &   (S-C2M1)*(CONT(I+N2)+CONT(I+N3)*(S-C1M1))+
     &    S*(CONT(I+N2)+CONT(I+N3)*(2.*S-C1M1-C2M1)))/HSOL
      RETURN
      END
C
C     END OF FUNCTION DONTR5
C
C ***********************************************************


