source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/atm_radau5.f @ 3734

Last change on this file since 3734 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 132.4 KB
Line 
1      SUBROUTINE INTEGRATE( TIN, TOUT )
2
3      IMPLICIT KPP_REAL (A-H,O-Z)
4      INCLUDE 'KPP_ROOT_Parameters.h'
5      INCLUDE 'KPP_ROOT_Global.h'
6
7C TIN - Start Time
8      KPP_REAL TIN
9C TOUT - End Time
10      KPP_REAL TOUT, H
11      SAVE H
12      DATA H /0.0d0/
13      INTEGER Nfun, Njac, Nstp, Nacc, Nrej, Ndec, Nsol
14      SAVE Nstp, Nacc, Nrej
15      DATA Nstp /0/
16      DATA Nacc /0/
17      DATA Nrej /0/
18      INTEGER i
19
20      PARAMETER (LWORK=5*NVAR*NVAR+12*NVAR+20,LIWORK=3*NVAR+20)
21      PARAMETER (LRCONT=4*NVAR+4)
22
23      KPP_REAL WORK(LWORK), RPAR(1)
24      INTEGER IWORK(LIWORK), IPAR(1)
25      COMMON /CONT/ ICONT(4),RCONT(LRCONT)
26      EXTERNAL FUNC_CHEM,JAC_CHEM,SOLOUT
27
28
29      ITOL=1     ! --- VECTOR TOLERANCES
30      IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
31      IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
32      IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
33      MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
34      MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
35      MLMAS=NVAR ! --- JACOBIAN IS A FULL MATRIX
36      MUMAS=NVAR ! --- JACOBIAN IS A FULL MATRIX
37
38      DO i=1,20
39        IWORK(i) = 0
40        WORK(i) = 0.D0
41      ENDDO
42
43      IWORK(3) = 8    !  Max no. of Newton iterations
44      IWORK(4) = 0    !  Starting values for Newton are interpolated (0) or zero (1)
45      IWORK(8) = 2    ! Gustaffson (1) or classic(2) controller
46      WORK(2)  = 0.9  ! Safety factor
47
48      CALL RADAU5(NVAR,FUNC_CHEM,TIN,VAR,TOUT,H,
49     &                  RTOL,ATOL,ITOL,
50     &                  JAC_CHEM ,IJAC,MLJAC,MUJAC,
51     &                  FUNC_CHEM ,IMAS,MLMAS,MUMAS,
52     &                  SOLOUT,IOUT,
53     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
54
55
56       Nfun = Nfun + IWORK(14)
57       Njac = Njac + IWORK(15)
58       Nstp = Nstp + IWORK(16)
59       Nacc = Nacc + IWORK(17)
60       Nrej = Nrej + IWORK(18)
61       Ndec = Ndec + IWORK(19)
62       Nsol = Nsol + IWORK(20)
63
64       print("('Nstep=',I5,' Nacc=',I6,' Nrej=',I6)"),Nstp, Nacc, Nrej
65
66      IF (IDID.LT.0) THEN
67        print *,'RADAU: Unsucessfull exit at T=',
68     &          TIN,' (IDID=',IDID,')'
69      ENDIF
70
71      RETURN
72      END
73
74
75      SUBROUTINE RADAU5(N,FCN,X,Y,XEND,H,
76     &                  RelTol,AbsTol,ITOL,
77     &                  JAC ,IJAC,MLJAC,MUJAC,
78     &                  MAS ,IMAS,MLMAS,MUMAS,
79     &                  SOLOUT,IOUT,
80     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
81C ----------------------------------------------------------
82C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
83C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS
84C                     M*Y'=F(X,Y).
85C     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I)
86C     OR EXPLICIT (M=I).
87C     THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA)
88C     OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT.
89C     C.F. SECTION IV.8
90C
91C     AUTHORS: E. HAIRER AND G. WANNER
92C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
93C              CH-1211 GENEVE 24, SWITZERLAND
94C              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
95C
96C     THIS CODE IS PART OF THE BOOK:
97C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
98C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
99C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
100C         SPRINGER-VERLAG (1991)
101C
102C     VERSION OF SEPTEMBER 30, 1995
103C
104C     INPUT PARAMETERS
105C     ----------------
106C     N           DIMENSION OF THE SYSTEM
107C
108C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
109C                 VALUE OF F(X,Y):
110C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR)
111C                    KPP_REAL X,Y(N),F(N)
112C                    F(1)=...   ETC.
113C                 RPAR, IPAR (SEE BELOW)
114C
115C     X           INITIAL X-VALUE
116C
117C     Y(N)        INITIAL VALUES FOR Y
118C
119C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
120C
121C     H           INITIAL STEP SIZE GUESS;
122C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
123C                 H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD.
124C                 THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS
125C                 QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).
126C
127C     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
128C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
129C
130C     ITOL        SWITCH FOR RelTol AND AbsTol:
131C                   ITOL=0: BOTH RelTol AND AbsTol ARE SCALARS.
132C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
133C                     Y(I) BELOW RelTol*ABS(Y(I))+AbsTol
134C                   ITOL=1: BOTH RelTol AND AbsTol ARE VECTORS.
135C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
136C                     RelTol(I)*ABS(Y(I))+AbsTol(I).
137C
138C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
139C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
140C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
141C                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
142C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
143C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR)
144C                    KPP_REAL X,Y(N),DFY(LDFY,N)
145C                    DFY(1,1)= ...
146C                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
147C                 FURNISHED BY THE CALLING PROGRAM.
148C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
149C                    BE FULL AND THE PARTIAL DERIVATIVES ARE
150C                    STORED IN DFY AS
151C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
152C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
153C                    THE PARTIAL DERIVATIVES ARE STORED
154C                    DIAGONAL-WISE AS
155C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
156C
157C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
158C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
159C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
160C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
161C
162C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
163C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
164C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
165C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
166C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
167C                       THE MAIN DIAGONAL).
168C
169C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
170C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
171C                 NEED NOT BE DEFINED IF MLJAC=N.
172C
173C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
174C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
175C
176C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
177C                 MATRIX M.
178C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
179C                 MATRIX AND NEEDS NOT TO BE DEFINED;
180C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
181C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
182C                    SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR)
183C                    KPP_REAL AM(LMAS,N)
184C                    AM(1,1)= ....
185C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
186C                    AS FULL MATRIX LIKE
187C                         AM(I,J) = M(I,J)
188C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
189C                    DIAGONAL-WISE AS
190C                         AM(I-J+MUMAS+1,J) = M(I,J).
191C
192C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
193C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
194C                       MATRIX, MAS IS NEVER CALLED.
195C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
196C
197C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
198C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
199C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
200C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
201C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
202C                       THE MAIN DIAGONAL).
203C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
204C
205C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
206C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
207C                 NEED NOT BE DEFINED IF MLMAS=N.
208C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
209C
210C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
211C                 NUMERICAL SOLUTION DURING INTEGRATION.
212C                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
213C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
214C                 IT MUST HAVE THE FORM
215C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,
216C                                       RPAR,IPAR,IRTRN)
217C                    KPP_REAL X,Y(N),CONT(LRC)
218C                    ....
219C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
220C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
221C                    THE FIRST GRID-POINT).
222C                 "XOLD" IS THE PRECEEDING GRID-POINT.
223C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
224C                    IS SET <0, RADAU5 RETURNS TO THE CALLING PROGRAM.
225C
226C          -----  CONTINUOUS OUTPUT: -----
227C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
228C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
229C                 THE FUNCTION
230C                        >>>   CONTR5(I,S,CONT,LRC)   <<<
231C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
232C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
233C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
234C                 DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE
235C                 DENSE OUTPUT FUNCTION IS USED.
236C
237C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
238C                    IOUT=0: SUBROUTINE IS NEVER CALLED
239C                    IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.
240C
241C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
242C                 WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS
243C                 FOR THE CODE. FOR STANDARD USE OF THE CODE
244C                 WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE
245C                 CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE.
246C                 WORK(21),..,WORK(LWORK) SERVE AS WORKING SPACE
247C                 FOR ALL VECTORS AND MATRICES.
248C                 "LWORK" MUST BE AT LEAST
249C                             N*(LJAC+LMAS+3*LE+12)+20
250C                 WHERE
251C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
252C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
253C                 AND
254C                    LMAS=0              IF IMAS=0
255C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
256C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
257C                 AND
258C                    LE=N               IF MLJAC=N (FULL JACOBIAN)
259C                    LE=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
260C
261C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
262C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
263C                 STORAGE REQUIREMENT IS
264C                             LWORK = 4*N*N+12*N+20.
265C                 IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST
266C                          N*(LJAC+12)+(N-M1)*(LMAS+3*LE)+20
267C                 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE
268C                 NUMBER N CAN BE REPLACED BY N-M1.
269C
270C     LWORK       DECLARED LENGTH OF ARRAY "WORK".
271C
272C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
273C                 IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS
274C                 FOR THE CODE. FOR STANDARD USE, SET IWORK(1),..,
275C                 IWORK(20) TO ZERO BEFORE CALLING.
276C                 IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA.
277C                 "LIWORK" MUST BE AT LEAST 3*N+20.
278C
279C     LIWORK      DECLARED LENGTH OF ARRAY "IWORK".
280C
281C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH
282C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
283C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES.
284C
285C ----------------------------------------------------------------------
286C
287C     SOPHISTICATED SETTING OF PARAMETERS
288C     -----------------------------------
289C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
290C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),...
291C              AS WELL AS IWORK(1),... DIFFERENT FROM ZERO.
292C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
293C
294C    IWORK(1)  IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN
295C              MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY
296C              ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN.
297C              IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N)
298C              AND NOT FOR IMPLICIT SYSTEMS (IMAS=1).
299C
300C    IWORK(2)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
301C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000.
302C
303C    IWORK(3)  THE MAXIMUM NUMBER OF NEWTON ITERATIONS FOR THE
304C              SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP.
305C              THE DEFAULT VALUE (FOR IWORK(3)=0) IS 7.
306C
307C    IWORK(4)  IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION
308C              IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD.
309C              IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED.
310C              THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS
311C              DIFFICULTIES WITH CONVERGENCE (THIS IS THE CASE WHEN
312C              NSTEP IS LARGER THAN NACCPT + NREJCT; SEE OUTPUT PARAM.).
313C              DEFAULT IS IWORK(4)=0.
314C
315C       THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR
316C       DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1.
317C       THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT
318C       THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER.
319C       IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE
320C       MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2.
321C
322C    IWORK(5)  DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR
323C              ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM.
324C              DEFAULT IWORK(5)=N.
325C
326C    IWORK(6)  DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0.
327C
328C    IWORK(7)  DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0.
329C
330C    IWORK(8)  SWITCH FOR STEP SIZE STRATEGY
331C              IF IWORK(8).EQ.1  MOD. PREDICTIVE CONTROLLER (GUSTAFSSON)
332C              IF IWORK(8).EQ.2  CLASSICAL STEP SIZE CONTROL
333C              THE DEFAULT VALUE (FOR IWORK(8)=0) IS IWORK(8)=1.
334C              THE CHOICE IWORK(8).EQ.1 SEEMS TO PRODUCE SAFER RESULTS;
335C              FOR SIMPLE PROBLEMS, THE CHOICE IWORK(8).EQ.2 PRODUCES
336C              OFTEN SLIGHTLY FASTER RUNS
337C
338C       IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT
339C            Y(I)' = Y(I+M2)   FOR  I=1,...,M1,
340C       WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME
341C       CAN BE ACHIEVED BY SETTING THE FOLLOWING TWO PARAMETERS. E.G.,
342C       FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE
343C       VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2.
344C       FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS:
345C       - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE
346C              JACOBIAN HAVE TO BE STORED
347C              IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL
348C                 DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J)
349C                FOR I=1,N-M1 AND J=1,N.
350C              ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM )
351C                 DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2)
352C                FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM.
353C       - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL
354C                0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM)
355C                     PARTIAL F(I+M1) / PARTIAL Y(J+K*M2),  I,J=1,M2
356C                    ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH
357C                    OF THESE MM+1 SUBMATRICES
358C       - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES
359C                NEED NOT BE DEFINED IF MLJAC=N-M1
360C       - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND
361C              NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
362C              IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF
363C              DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX.
364C              IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL
365C                 AM(I,J) = M(I+M1,J+M1)     FOR I=1,N-M1 AND J=1,N-M1.
366C              ELSE, THE MASS MATRIX IS BANDED
367C                 AM(I-J+MUMAS+1,J) = M(I+M1,J+M1)
368C       - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL
369C                0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX
370C       - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX
371C                NEED NOT BE DEFINED IF MLMAS=N-M1
372C
373C    IWORK(9)  THE VALUE OF M1.  DEFAULT M1=0.
374C
375C    IWORK(10) THE VALUE OF M2.  DEFAULT M2=M1.
376C
377C ----------
378C
379C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
380C
381C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION,
382C              DEFAULT 0.9D0.
383C
384C    WORK(3)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
385C              INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS
386C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER
387C              (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE TO
388C              COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP.
389C              DEFAULT 0.001D0.
390C
391C    WORK(4)   STOPPING CRITERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
392C              SMALLER VALUES OF WORK(4) MAKE THE CODE SLOWER, BUT SAFER.
393C              DEFAULT 0.03D0.
394C
395C    WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE
396C              STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A
397C              LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR
398C              LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE
399C              WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS
400C              WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD.
401C              DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 .
402C
403C    WORK(7)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
404C
405C    WORK(8), WORK(9)   PARAMETERS FOR STEP SIZE SELECTION
406C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
407C                 WORK(8) <= HNEW/HOLD <= WORK(9)
408C              DEFAULT VALUES: WORK(8)=0.2D0, WORK(9)=8.D0
409C
410C-----------------------------------------------------------------------
411C
412C     OUTPUT PARAMETERS
413C     -----------------
414C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
415C                 (AFTER SUCCESSFUL RETURN X=XEND).
416C
417C     Y(N)        NUMERICAL SOLUTION AT X
418C
419C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
420C
421C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
422C                   IDID= 1  COMPUTATION SUCCESSFUL,
423C                   IDID= 2  COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
424C                   IDID=-1  INPUT IS NOT CONSISTENT,
425C                   IDID=-2  LARGER NMAX IS NEEDED,
426C                   IDID=-3  STEP SIZE BECOMES TOO SMALL,
427C                   IDID=-4  MATRIX IS REPEATEDLY SINGULAR.
428C
429C   IWORK(14)  NFCN    NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
430C                      EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
431C   IWORK(15)  NJAC    NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
432C                      OR NUMERICALLY)
433C   IWORK(16)  NSTEP   NUMBER OF COMPUTED STEPS
434C   IWORK(17)  NACCPT  NUMBER OF ACCEPTED STEPS
435C   IWORK(18)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
436C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
437C   IWORK(19)  NDEC    NUMBER OF LU-DECOMPOSITIONS OF BOTH MATRICES
438C   IWORK(20)  NSOL    NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF BOTH
439C                      SYSTEMS; THE NSTEP FORWARD-BACKWARD SUBSTITUTIONS,
440C                      NEEDED FOR STEP SIZE SELECTION, ARE NOT COUNTED
441C-----------------------------------------------------------------------
442C *** *** *** *** *** *** *** *** *** *** *** *** ***
443C          DECLARATIONS
444C *** *** *** *** *** *** *** *** *** *** *** *** ***
445      IMPLICIT KPP_REAL (A-H,O-Z)
446      DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK)
447      DIMENSION RPAR(*),IPAR(*)
448      LOGICAL IMPLCT,JBAND,ARRET,STARTN,PRED
449      EXTERNAL FCN,JAC,MAS,SOLOUT
450C *** *** *** *** *** *** ***
451C        SETTING THE PARAMETERS
452C *** *** *** *** *** *** ***
453       NFCN=0
454       NJAC=0
455       NSTEP=0
456       NACCPT=0
457       NREJCT=0
458       NDEC=0
459       NSOL=0
460       ARRET=.FALSE.
461C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
462      IF (IWORK(2).EQ.0) THEN
463         NMAX=100000
464      ELSE
465         NMAX=IWORK(2)
466         IF (NMAX.LE.0) THEN
467            WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
468            ARRET=.TRUE.
469         END IF
470      END IF
471C -------- NIT    MAXIMAL NUMBER OF NEWTON ITERATIONS
472      IF (IWORK(3).EQ.0) THEN
473         NIT=7
474      ELSE
475         NIT=IWORK(3)
476         IF (NIT.LE.0) THEN
477            WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
478            ARRET=.TRUE.
479         END IF
480      END IF
481C -------- STARTN  SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS
482      IF(IWORK(4).EQ.0)THEN
483         STARTN=.FALSE.
484      ELSE
485         STARTN=.TRUE.
486      END IF
487C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS
488      NIND1=IWORK(5)
489      NIND2=IWORK(6)
490      NIND3=IWORK(7)
491      IF (NIND1.EQ.0) NIND1=N
492      IF (NIND1+NIND2+NIND3.NE.N) THEN
493       WRITE(6,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',NIND1,NIND2,NIND3
494       ARRET=.TRUE.
495      END IF
496C -------- PRED   STEP SIZE CONTROL
497      IF(IWORK(8).LE.1)THEN
498         PRED=.TRUE.
499      ELSE
500         PRED=.FALSE.
501      END IF
502C -------- PARAMETER FOR SECOND ORDER EQUATIONS
503      M1=IWORK(9)
504      M2=IWORK(10)
505      NM1=N-M1
506      IF (M1.EQ.0) M2=N
507      IF (M2.EQ.0) M2=M1
508      IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN
509       WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
510       ARRET=.TRUE.
511      END IF
512C -------- UROUND   SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0
513      IF (WORK(1).EQ.0.0D0) THEN
514         UROUND=1.0D-16
515      ELSE
516         UROUND=WORK(1)
517         IF (UROUND.LE.1.0D-19.OR.UROUND.GE.1.0D0) THEN
518            WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
519            ARRET=.TRUE.
520         END IF
521      END IF
522C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
523      IF (WORK(2).EQ.0.0D0) THEN
524         SAFE=0.9D0
525      ELSE
526         SAFE=WORK(2)
527         IF (SAFE.LE.0.001D0.OR.SAFE.GE.1.0D0) THEN
528            WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
529            ARRET=.TRUE.
530         END IF
531      END IF
532C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
533      IF (WORK(3).EQ.0.D0) THEN
534         THET=0.001D0
535      ELSE
536         THET=WORK(3)
537         IF (THET.LE.0.0D0.OR.THET.GE.1.0D0) THEN
538            WRITE(6,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3)
539            ARRET=.TRUE.
540         END IF
541      END IF
542C --- FNEWT   STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
543      IF (WORK(4).EQ.0.D0) THEN
544         FNEWT=0.03D0
545      ELSE
546         FNEWT=WORK(4)
547         IF (FNEWT.LE.UROUND) THEN
548            WRITE(6,*)' CURIOUS INPUT FOR WORK(4)=',WORK(4)
549            ARRET=.TRUE.
550         END IF
551      END IF
552C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
553      IF (WORK(5).EQ.0.D0) THEN
554         QUOT1=1.D0
555      ELSE
556         QUOT1=WORK(5)
557      END IF
558      IF (WORK(6).EQ.0.D0) THEN
559         QUOT2=1.2D0
560      ELSE
561         QUOT2=WORK(6)
562      END IF
563      IF (QUOT1.GT.1.0D0.OR.QUOT2.LT.1.0D0) THEN
564         WRITE(6,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2
565         ARRET=.TRUE.
566      END IF
567C -------- MAXIMAL STEP SIZE
568      IF (WORK(7).EQ.0.D0) THEN
569         HMAX=XEND-X
570      ELSE
571         HMAX=WORK(7)
572      END IF
573C -------  FACL,FACR     PARAMETERS FOR STEP SIZE SELECTION
574      IF(WORK(8).EQ.0.D0)THEN
575         FACL=5.D0
576      ELSE
577         FACL=1.D0/WORK(8)
578      END IF
579      IF(WORK(9).EQ.0.D0)THEN
580         FACR=1.D0/8.0D0
581      ELSE
582         FACR=1.D0/WORK(9)
583      END IF
584      IF (FACL.LT.1.0D0.OR.FACR.GT.1.0D0) THEN
585            WRITE(6,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9)
586            ARRET=.TRUE.
587         END IF
588C --------- CHECK IF TOLERANCES ARE O.K.
589      IF (ITOL.EQ.0) THEN
590          IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
591              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
592              ARRET=.TRUE.
593          END IF
594      ELSE
595          DO I=1,N
596          IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
597              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
598              ARRET=.TRUE.
599          END IF
600          END DO
601      END IF
602C *** *** *** *** *** *** *** *** *** *** *** *** ***
603C         COMPUTATION OF ARRAY ENTRIES
604C *** *** *** *** *** *** *** *** *** *** *** *** ***
605C ---- IMPLICIT, BANDED OR NOT ?
606      IMPLCT=IMAS.NE.0
607      JBAND=MLJAC.LT.NM1
608C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
609C -- JACOBIAN  AND  MATRICES E1, E2
610      IF (JBAND) THEN
611         LDJAC=MLJAC+MUJAC+1
612         LDE1=MLJAC+LDJAC
613      ELSE
614         MLJAC=NM1
615         MUJAC=NM1
616         LDJAC=NM1
617         LDE1=NM1
618      END IF
619C -- MASS MATRIX
620      IF (IMPLCT) THEN
621          IF (MLMAS.NE.NM1) THEN
622              LDMAS=MLMAS+MUMAS+1
623              IF (JBAND) THEN
624                 IJOB=4
625              ELSE
626                 IJOB=3
627              END IF
628          ELSE
629              LDMAS=NM1
630              IJOB=5
631          END IF
632C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC"
633          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
634             WRITE (6,*) 'BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF
635     & "JAC"'
636            ARRET=.TRUE.
637          END IF
638      ELSE
639          LDMAS=0
640          IF (JBAND) THEN
641             IJOB=2
642          ELSE
643             IJOB=1
644             IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7
645          END IF
646      END IF
647      LDMAS2=MAX(1,LDMAS)
648C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
649      IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN
650         WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
651     &FULL JACOBIAN'
652         ARRET=.TRUE.
653      END IF
654C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
655      IEZ1=21
656      IEZ2=IEZ1+N
657      IEZ3=IEZ2+N
658      IEY0=IEZ3+N
659      IESCAL=IEY0+N
660      IEF1=IESCAL+N
661      IEF2=IEF1+N
662      IEF3=IEF2+N
663      IECON=IEF3+N
664      IEJAC=IECON+4*N
665      IEMAS=IEJAC+N*LDJAC
666      IEE1=IEMAS+NM1*LDMAS
667      IEE2R=IEE1+NM1*LDE1
668      IEE2I=IEE2R+NM1*LDE1
669C ------ TOTAL STORAGE REQUIREMENT -----------
670      ISTORE=IEE2I+NM1*LDE1-1
671      IF(ISTORE.GT.LWORK)THEN
672         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
673         ARRET=.TRUE.
674      END IF
675C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
676      IEIP1=21
677      IEIP2=IEIP1+NM1
678      IEIPH=IEIP2+NM1
679C --------- TOTAL REQUIREMENT ---------------
680      ISTORE=IEIPH+NM1-1
681      IF (ISTORE.GT.LIWORK) THEN
682         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
683         ARRET=.TRUE.
684      END IF
685C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
686      IF (ARRET) THEN
687         IDID=-1
688         RETURN
689      END IF
690C -------- CALL TO CORE INTEGRATOR ------------
691      CALL RADCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
692     &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
693     &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN,
694     &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,
695     &   IMPLCT,JBAND,LDJAC,LDE1,LDMAS2,WORK(IEZ1),WORK(IEZ2),
696     &   WORK(IEZ3),WORK(IEY0),WORK(IESCAL),WORK(IEF1),WORK(IEF2),
697     &   WORK(IEF3),WORK(IEJAC),WORK(IEE1),WORK(IEE2R),WORK(IEE2I),
698     &   WORK(IEMAS),IWORK(IEIP1),IWORK(IEIP2),IWORK(IEIPH),
699     &   WORK(IECON),NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR)
700      IWORK(14)=NFCN
701      IWORK(15)=NJAC
702      IWORK(16)=NSTEP
703      IWORK(17)=NACCPT
704      IWORK(18)=NREJCT
705      IWORK(19)=NDEC
706      IWORK(20)=NSOL
707C ----------- RETURN -----------
708      RETURN
709      END
710C
711C     END OF SUBROUTINE RADAU5
712C
713C ***********************************************************
714C
715      SUBROUTINE RADCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
716     &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
717     &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN,
718     &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,
719     &   IMPLCT,BANDED,LDJAC,LDE1,LDMAS,Z1,Z2,Z3,
720     &   Y0,SCAL,F1,F2,F3,FJAC,E1,E2R,E2I,FMAS,IP1,IP2,IPHES,
721     &   CONT,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR)
722C ----------------------------------------------------------
723C     CORE INTEGRATOR FOR RADAU5
724C     PARAMETERS SAME AS IN RADAU5 WITH WORKSPACE ADDED
725C ----------------------------------------------------------
726C         DECLARATIONS
727C ----------------------------------------------------------
728      IMPLICIT KPP_REAL (A-H,O-Z)
729      DIMENSION Y(N),Z1(N),Z2(N),Z3(N),Y0(N),SCAL(N),F1(N),F2(N),F3(N)
730      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),CONT(4*N)
731      DIMENSION E1(LDE1,NM1),E2R(LDE1,NM1),E2I(LDE1,NM1)
732      DIMENSION AbsTol(*),RelTol(*),RPAR(*),IPAR(*)
733      INTEGER IP1(NM1),IP2(NM1),IPHES(NM1)
734      COMMON /CONRA5/NN,NN2,NN3,NN4,XSOL,HSOL,C2M1,C1M1
735      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
736      LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN,CALHES
737      LOGICAL INDEX1,INDEX2,INDEX3,LAST,PRED
738      EXTERNAL FCN, JAC, MAS, SOLOUT
739C *** *** *** *** *** *** ***
740C  INITIALISATIONS
741C *** *** *** *** *** *** ***
742C --------- DUPLIFY N FOR COMMON BLOCK CONT -----
743      NN=N
744      NN2=2*N
745      NN3=3*N
746      LRC=4*N
747C -------- CHECK THE INDEX OF THE PROBLEM -----
748      INDEX1=NIND1.NE.0
749      INDEX2=NIND2.NE.0
750      INDEX3=NIND3.NE.0
751C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
752      IF (IMPLCT) CALL MAS(NM1,FMAS,LDMAS,RPAR,IPAR)
753C ---------- CONSTANTS ---------
754      SQ6=DSQRT(6.D0)
755      C1=(4.D0-SQ6)/10.D0
756      C2=(4.D0+SQ6)/10.D0
757      C1M1=C1-1.D0
758      C2M1=C2-1.D0
759      C1MC2=C1-C2
760      DD1=-(13.D0+7.D0*SQ6)/3.D0
761      DD2=(-13.D0+7.D0*SQ6)/3.D0
762      DD3=-1.D0/3.D0
763      U1=(6.D0+81.D0**(1.D0/3.D0)-9.D0**(1.D0/3.D0))/30.D0
764      ALPH=(12.D0-81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))/60.D0
765      BETA=(81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))*DSQRT(3.D0)/60.D0
766      CNO=ALPH**2+BETA**2
767      U1=1.0D0/U1
768      ALPH=ALPH/CNO
769      BETA=BETA/CNO
770      T11=9.1232394870892942792D-02
771      T12=-0.14125529502095420843D0
772      T13=-3.0029194105147424492D-02
773      T21=0.24171793270710701896D0
774      T22=0.20412935229379993199D0
775      T23=0.38294211275726193779D0
776      T31=0.96604818261509293619D0
777      TI11=4.3255798900631553510D0
778      TI12=0.33919925181580986954D0
779      TI13=0.54177053993587487119D0
780      TI21=-4.1787185915519047273D0
781      TI22=-0.32768282076106238708D0
782      TI23=0.47662355450055045196D0
783      TI31=-0.50287263494578687595D0
784      TI32=2.5719269498556054292D0
785      TI33=-0.59603920482822492497D0
786      IF (M1.GT.0) IJOB=IJOB+10
787      POSNEG=SIGN(1.D0,XEND-X)
788      HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
789      IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
790      H=MIN(ABS(H),HMAXN)
791      H=SIGN(H,POSNEG)
792      HOLD=H
793      REJECT=.FALSE.
794      FIRST=.TRUE.
795      LAST=.FALSE.
796      IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN
797         H=XEND-X
798         LAST=.TRUE.
799      END IF
800      FACCON=1.D0
801      CFAC=SAFE*(1+2*NIT)
802      NSING=0
803      XOLD=X
804      IF (IOUT.NE.0) THEN
805          IRTRN=1
806          NRSOL=1
807          XOSOL=XOLD
808          XSOL=X
809          DO I=1,N
810             CONT(I)=Y(I)
811          END DO
812          NSOLU=N
813          HSOL=HOLD
814          CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU,
815     &                RPAR,IPAR,IRTRN)
816          IF (IRTRN.LT.0) GOTO 179
817      END IF
818      MLE=MLJAC
819      MUE=MUJAC
820      MBJAC=MLJAC+MUJAC+1
821      MBB=MLMAS+MUMAS+1
822      MDIAG=MLE+MUE+1
823      MDIFF=MLE+MUE-MUMAS
824      MBDIAG=MUMAS+1
825      N2=2*N
826      N3=3*N
827      IF (ITOL.EQ.0) THEN
828          DO I=1,N
829             SCAL(I)=AbsTol(1)+RelTol(1)*ABS(Y(I))
830          END DO
831      ELSE
832          DO I=1,N
833             SCAL(I)=AbsTol(I)+RelTol(I)*ABS(Y(I))
834          END DO
835      END IF
836      HHFAC=H
837      CALL FCN(N,X,Y,Y0,RPAR,IPAR)
838      NFCN=NFCN+1
839C --- BASIC INTEGRATION STEP
840  10  CONTINUE
841C *** *** *** *** *** *** ***
842C  COMPUTATION OF THE JACOBIAN
843C *** *** *** *** *** *** ***
844      NJAC=NJAC+1
845      IF (IJAC.EQ.0) THEN
846C --- COMPUTE JACOBIAN MATRIX NUMERICALLY
847         IF (BANDED) THEN
848C --- JACOBIAN IS BANDED
849            MUJACP=MUJAC+1
850            MD=MIN(MBJAC,M2)
851            DO MM=1,M1/M2+1
852               DO K=1,MD
853                  J=K+(MM-1)*M2
854 12               F1(J)=Y(J)
855                  F2(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
856                  Y(J)=Y(J)+F2(J)
857                  J=J+MD
858                  IF (J.LE.MM*M2) GOTO 12
859                  CALL FCN(N,X,Y,CONT,RPAR,IPAR)
860                  J=K+(MM-1)*M2
861                  J1=K
862                  LBEG=MAX(1,J1-MUJAC)+M1
863 14               LEND=MIN(M2,J1+MLJAC)+M1
864                  Y(J)=F1(J)
865                  MUJACJ=MUJACP-J1-M1
866                  DO L=LBEG,LEND
867                     FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/F2(J)
868                  END DO
869                  J=J+MD
870                  J1=J1+MD
871                  LBEG=LEND+1
872                  IF (J.LE.MM*M2) GOTO 14
873               END DO
874            END DO
875         ELSE
876C --- JACOBIAN IS FULL
877            DO I=1,N
878               YSAFE=Y(I)
879               DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
880               Y(I)=YSAFE+DELT
881               CALL FCN(N,X,Y,CONT,RPAR,IPAR)
882               DO J=M1+1,N
883                 FJAC(J-M1,I)=(CONT(J)-Y0(J))/DELT
884               END DO
885               Y(I)=YSAFE
886            END DO
887         END IF
888      ELSE
889C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
890      CALL JAC_CHEM(N,X,Y,FJAC,LDJAC,RPAR,IPAR)
891      END IF
892      CALJAC=.TRUE.
893      CALHES=.TRUE.
894  20  CONTINUE
895C --- COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS
896      FAC1=U1/H
897      ALPHN=ALPH/H
898      BETAN=BETA/H
899      CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
900     &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
901      IF (IER.NE.0) GOTO 78
902      CALL DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
903     &            M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB)
904      IF (IER.NE.0) GOTO 78
905      NDEC=NDEC+1
906  30  CONTINUE
907      NSTEP=NSTEP+1
908      IF (NSTEP.GT.NMAX) GOTO 178
909      IF (0.1D0*ABS(H).LE.ABS(X)*UROUND) GOTO 177
910          IF (INDEX2) THEN
911             DO I=NIND1+1,NIND1+NIND2
912                SCAL(I)=SCAL(I)/HHFAC
913             END DO
914          END IF
915          IF (INDEX3) THEN
916             DO I=NIND1+NIND2+1,NIND1+NIND2+NIND3
917                SCAL(I)=SCAL(I)/(HHFAC*HHFAC)
918             END DO
919          END IF
920      XPH=X+H
921C *** *** *** *** *** *** ***
922C  STARTING VALUES FOR NEWTON ITERATION
923C *** *** *** *** *** *** ***
924      IF (FIRST.OR.STARTN) THEN
925         DO I=1,N
926            Z1(I)=0.D0
927            Z2(I)=0.D0
928            Z3(I)=0.D0
929            F1(I)=0.D0
930            F2(I)=0.D0
931            F3(I)=0.D0
932         END DO
933      ELSE
934         C3Q=H/HOLD
935         C1Q=C1*C3Q
936         C2Q=C2*C3Q
937         DO I=1,N
938            AK1=CONT(I+N)
939            AK2=CONT(I+N2)
940            AK3=CONT(I+N3)
941            Z1I=C1Q*(AK1+(C1Q-C2M1)*(AK2+(C1Q-C1M1)*AK3))
942            Z2I=C2Q*(AK1+(C2Q-C2M1)*(AK2+(C2Q-C1M1)*AK3))
943            Z3I=C3Q*(AK1+(C3Q-C2M1)*(AK2+(C3Q-C1M1)*AK3))
944            Z1(I)=Z1I
945            Z2(I)=Z2I
946            Z3(I)=Z3I
947            F1(I)=TI11*Z1I+TI12*Z2I+TI13*Z3I
948            F2(I)=TI21*Z1I+TI22*Z2I+TI23*Z3I
949            F3(I)=TI31*Z1I+TI32*Z2I+TI33*Z3I
950         END DO
951      END IF
952C *** *** *** *** *** *** ***
953C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
954C *** *** *** *** *** *** ***
955            NEWT=0
956            FACCON=MAX(FACCON,UROUND)**0.8D0
957            THETA=ABS(THET)
958  40        CONTINUE
959            IF (NEWT.GE.NIT) GOTO 78
960C ---     COMPUTE THE RIGHT-HAND SIDE
961            DO I=1,N
962               CONT(I)=Y(I)+Z1(I)
963            END DO
964            CALL FCN(N,X+C1*H,CONT,Z1,RPAR,IPAR)
965            DO I=1,N
966               CONT(I)=Y(I)+Z2(I)
967            END DO
968            CALL FCN(N,X+C2*H,CONT,Z2,RPAR,IPAR)
969            DO I=1,N
970               CONT(I)=Y(I)+Z3(I)
971            END DO
972            CALL FCN(N,XPH,CONT,Z3,RPAR,IPAR)
973            NFCN=NFCN+3
974C ---     KppSolve THE LINEAR SYSTEMS
975           DO I=1,N
976              A1=Z1(I)
977              A2=Z2(I)
978              A3=Z3(I)
979              Z1(I)=TI11*A1+TI12*A2+TI13*A3
980              Z2(I)=TI21*A1+TI22*A2+TI23*A3
981              Z3(I)=TI31*A1+TI32*A2+TI33*A3
982           END DO
983        CALL SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
984     &          M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3,
985     &          F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB)
986            NSOL=NSOL+1
987            NEWT=NEWT+1
988            DYNO=0.D0
989            DO I=1,N
990               DENOM=SCAL(I)
991               DYNO=DYNO+(Z1(I)/DENOM)**2+(Z2(I)/DENOM)**2
992     &          +(Z3(I)/DENOM)**2
993            END DO
994            DYNO=DSQRT(DYNO/N3)
995C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
996            IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN
997                THQ=DYNO/DYNOLD
998                IF (NEWT.EQ.2) THEN
999                   THETA=THQ
1000                ELSE
1001                   THETA=SQRT(THQ*THQOLD)
1002                END IF
1003                THQOLD=THQ
1004                IF (THETA.LT.0.99D0) THEN
1005                    FACCON=THETA/(1.0D0-THETA)
1006                    DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT
1007                    IF (DYTH.GE.1.0D0) THEN
1008                         QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH))
1009                         HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
1010                         H=HHFAC*H
1011                         REJECT=.TRUE.
1012                         LAST=.FALSE.
1013                         IF (CALJAC) GOTO 20
1014                         GOTO 10
1015                    END IF
1016                ELSE
1017                    GOTO 78
1018                END IF
1019            END IF
1020            DYNOLD=MAX(DYNO,UROUND)
1021            DO I=1,N
1022               F1I=F1(I)+Z1(I)
1023               F2I=F2(I)+Z2(I)
1024               F3I=F3(I)+Z3(I)
1025               F1(I)=F1I
1026               F2(I)=F2I
1027               F3(I)=F3I
1028               Z1(I)=T11*F1I+T12*F2I+T13*F3I
1029               Z2(I)=T21*F1I+T22*F2I+T23*F3I
1030               Z3(I)=T31*F1I+    F2I
1031            END DO
1032            IF (FACCON*DYNO.GT.FNEWT) GOTO 40
1033C --- ERROR ESTIMATION
1034      CALL ESTRAD (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1035     &          H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,
1036     &          E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR,
1037     &          FIRST,REJECT,FAC1,RPAR,IPAR)
1038C --- COMPUTATION OF HNEW
1039C --- WE REQUIRE .2<=HNEW/H<=8.
1040      FAC=MIN(SAFE,CFAC/(NEWT+2*NIT))
1041      QUOT=MAX(FACR,MIN(FACL,ERR**.25D0/FAC))
1042      HNEW=H/QUOT
1043C *** *** *** *** *** *** ***
1044C  IS THE ERROR SMALL ENOUGH ?
1045C *** *** *** *** *** *** ***
1046      IF (ERR.LT.1.D0) THEN
1047C --- STEP IS ACCEPTED
1048         FIRST=.FALSE.
1049         NACCPT=NACCPT+1
1050         IF (PRED) THEN
1051C       --- PREDICTIVE CONTROLLER OF GUSTAFSSON
1052            IF (NACCPT.GT.1) THEN
1053               FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE
1054               FACGUS=MAX(FACR,MIN(FACL,FACGUS))
1055               QUOT=MAX(QUOT,FACGUS)
1056               HNEW=H/QUOT
1057            END IF
1058            HACC=H
1059            ERRACC=MAX(1.0D-2,ERR)
1060         END IF
1061         XOLD=X
1062         HOLD=H
1063         X=XPH
1064         DO I=1,N
1065            Y(I)=Y(I)+Z3(I)
1066            Z2I=Z2(I)
1067            Z1I=Z1(I)
1068            CONT(I+N)=(Z2I-Z3(I))/C2M1
1069            AK=(Z1I-Z2I)/C1MC2
1070            ACONT3=Z1I/C1
1071            ACONT3=(AK-ACONT3)/C2
1072            CONT(I+N2)=(AK-CONT(I+N))/C1M1
1073            CONT(I+N3)=CONT(I+N2)-ACONT3
1074         END DO
1075         IF (ITOL.EQ.0) THEN
1076             DO I=1,N
1077                SCAL(I)=AbsTol(1)+RelTol(1)*ABS(Y(I))
1078             END DO
1079         ELSE
1080             DO I=1,N
1081                SCAL(I)=AbsTol(I)+RelTol(I)*ABS(Y(I))
1082             END DO
1083         END IF
1084         IF (IOUT.NE.0) THEN
1085             NRSOL=NACCPT+1
1086             XSOL=X
1087             XOSOL=XOLD
1088             DO I=1,N
1089                CONT(I)=Y(I)
1090             END DO
1091             NSOLU=N
1092             HSOL=HOLD
1093             CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU,
1094     &                   RPAR,IPAR,IRTRN)
1095             IF (IRTRN.LT.0) GOTO 179
1096         END IF
1097         CALJAC=.FALSE.
1098         IF (LAST) THEN
1099            H=HOPT
1100            IDID=1
1101            RETURN
1102         END IF
1103         CALL FCN(N,X,Y,Y0,RPAR,IPAR)
1104         NFCN=NFCN+1
1105         HNEW=POSNEG*MIN(ABS(HNEW),HMAXN)
1106         HOPT=HNEW
1107         HOPT=MIN(H,HNEW)
1108         IF (REJECT) HNEW=POSNEG*MIN(ABS(HNEW),ABS(H))
1109         REJECT=.FALSE.
1110         IF ((X+HNEW/QUOT1-XEND)*POSNEG.GE.0.D0) THEN
1111            H=XEND-X
1112            LAST=.TRUE.
1113         ELSE
1114            QT=HNEW/H
1115            HHFAC=H
1116            IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30
1117            H=HNEW
1118         END IF
1119         HHFAC=H
1120         IF (THETA.LE.THET) GOTO 20
1121         GOTO 10
1122      ELSE
1123C --- STEP IS REJECTED
1124         REJECT=.TRUE.
1125         LAST=.FALSE.
1126         IF (FIRST) THEN
1127             H=H*0.1D0
1128             HHFAC=0.1D0
1129         ELSE
1130             HHFAC=HNEW/H
1131             H=HNEW
1132         END IF
1133         IF (NACCPT.GE.1) NREJCT=NREJCT+1
1134         IF (CALJAC) GOTO 20
1135         GOTO 10
1136      END IF
1137C --- UNEXPECTED STEP-REJECTION
1138  78  CONTINUE
1139      IF (IER.NE.0) THEN
1140          NSING=NSING+1
1141          IF (NSING.GE.5) GOTO 176
1142      END IF
1143      H=H*0.5D0
1144      HHFAC=0.5D0
1145      REJECT=.TRUE.
1146      LAST=.FALSE.
1147      IF (CALJAC) GOTO 20
1148      GOTO 10
1149C --- FAIL EXIT
1150 176  CONTINUE
1151      WRITE(6,979)X
1152      WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER
1153      IDID=-4
1154      RETURN
1155 177  CONTINUE
1156      WRITE(6,979)X
1157      WRITE(6,*) ' STEP SIZE T0O SMALL, H=',H
1158      IDID=-3
1159      RETURN
1160 178  CONTINUE
1161      WRITE(6,979)X
1162      WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'
1163      IDID=-2
1164      RETURN
1165C --- EXIT CAUSED BY SOLOUT
1166 179  CONTINUE
1167      WRITE(6,979)X
1168 979  FORMAT(' EXIT OF RADAU5 AT X=',E18.4)
1169      IDID=2
1170      RETURN
1171      END
1172C
1173C     END OF SUBROUTINE RADCOR
1174C
1175C ***********************************************************
1176C
1177      KPP_REAL FUNCTION CONTR5(I,X,CONT,LRC)
1178C ----------------------------------------------------------
1179C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
1180C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
1181C     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
1182C     THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU5).
1183C ----------------------------------------------------------
1184      IMPLICIT KPP_REAL (A-H,O-Z)
1185      DIMENSION CONT(LRC)
1186      COMMON /CONRA5/NN,NN2,NN3,NN4,XSOL,HSOL,C2M1,C1M1
1187      S=(X-XSOL)/HSOL
1188      CONTR5=CONT(I)+S*(CONT(I+NN)+(S-C2M1)*(CONT(I+NN2)
1189     &     +(S-C1M1)*CONT(I+NN3)))
1190      RETURN
1191      END
1192C
1193C     END OF FUNCTION CONTR5
1194C
1195C ***********************************************************
1196C ******************************************
1197C     VERSION OF SEPTEMBER 18, 1995
1198C ******************************************
1199C
1200      SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
1201     &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
1202      IMPLICIT KPP_REAL (A-H,O-Z)
1203      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
1204     &          IP1(NM1),IPHES(N)
1205      LOGICAL CALHES
1206      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1207C
1208      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
1209C
1210C -----------------------------------------------------------
1211C
1212   1  CONTINUE
1213C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1214      DO J=1,N
1215         DO  I=1,N
1216            E1(I,J)=-FJAC(I,J)
1217         END DO
1218         E1(J,J)=E1(J,J)+FAC1
1219      END DO
1220      CALL DEC (N,LDE1,E1,IP1,IER)
1221      RETURN
1222C
1223C -----------------------------------------------------------
1224C
1225  11  CONTINUE
1226C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1227      DO J=1,NM1
1228         JM1=J+M1
1229         DO I=1,NM1
1230            E1(I,J)=-FJAC(I,JM1)
1231         END DO
1232         E1(J,J)=E1(J,J)+FAC1
1233      END DO
1234 45   MM=M1/M2
1235      DO J=1,M2
1236         DO I=1,NM1
1237            SUM=0.D0
1238            DO K=0,MM-1
1239               SUM=(SUM+FJAC(I,J+K*M2))/FAC1
1240            END DO
1241            E1(I,J)=E1(I,J)-SUM
1242         END DO
1243      END DO
1244      CALL DEC (NM1,LDE1,E1,IP1,IER)
1245      RETURN
1246C
1247C -----------------------------------------------------------
1248C
1249   2  CONTINUE
1250C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1251      DO J=1,N
1252         DO I=1,MBJAC
1253            E1(I+MLE,J)=-FJAC(I,J)
1254         END DO
1255         E1(MDIAG,J)=E1(MDIAG,J)+FAC1
1256      END DO
1257      CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER)
1258      RETURN
1259C
1260C -----------------------------------------------------------
1261C
1262  12  CONTINUE
1263C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1264      DO J=1,NM1
1265         JM1=J+M1
1266         DO I=1,MBJAC
1267            E1(I+MLE,J)=-FJAC(I,JM1)
1268         END DO
1269         E1(MDIAG,J)=E1(MDIAG,J)+FAC1
1270      END DO
1271  46  MM=M1/M2
1272      DO J=1,M2
1273         DO I=1,MBJAC
1274            SUM=0.D0
1275            DO K=0,MM-1
1276               SUM=(SUM+FJAC(I,J+K*M2))/FAC1
1277            END DO
1278            E1(I+MLE,J)=E1(I+MLE,J)-SUM
1279         END DO
1280      END DO
1281      CALL DECB (NM1,LDE1,E1,MLE,MUE,IP1,IER)
1282      RETURN
1283C
1284C -----------------------------------------------------------
1285C
1286   3  CONTINUE
1287C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1288      DO J=1,N
1289         DO I=1,N
1290            E1(I,J)=-FJAC(I,J)
1291         END DO
1292         DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS)
1293            E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J)
1294         END DO
1295      END DO
1296      CALL DEC (N,LDE1,E1,IP1,IER)
1297      RETURN
1298C
1299C -----------------------------------------------------------
1300C
1301  13  CONTINUE
1302C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1303      DO J=1,NM1
1304         JM1=J+M1
1305         DO I=1,NM1
1306            E1(I,J)=-FJAC(I,JM1)
1307         END DO
1308         DO I=MAX(1,J-MUMAS),MIN(NM1,J+MLMAS)
1309            E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J)
1310         END DO
1311      END DO
1312      GOTO 45
1313C
1314C -----------------------------------------------------------
1315C
1316   4  CONTINUE
1317C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
1318      DO J=1,N
1319         DO I=1,MBJAC
1320            E1(I+MLE,J)=-FJAC(I,J)
1321         END DO
1322         DO I=1,MBB
1323            IB=I+MDIFF
1324            E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J)
1325         END DO
1326      END DO
1327      CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER)
1328      RETURN
1329C
1330C -----------------------------------------------------------
1331C
1332  14  CONTINUE
1333C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
1334      DO J=1,NM1
1335         JM1=J+M1
1336         DO I=1,MBJAC
1337            E1(I+MLE,J)=-FJAC(I,JM1)
1338         END DO
1339         DO I=1,MBB
1340            IB=I+MDIFF
1341            E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J)
1342         END DO
1343      END DO
1344      GOTO 46
1345C
1346C -----------------------------------------------------------
1347C
1348   5  CONTINUE
1349C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
1350      DO J=1,N
1351         DO I=1,N
1352            E1(I,J)=FMAS(I,J)*FAC1-FJAC(I,J)
1353         END DO
1354      END DO
1355      CALL DEC (N,LDE1,E1,IP1,IER)
1356      RETURN
1357C
1358C -----------------------------------------------------------
1359C
1360  15  CONTINUE
1361C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1362      DO J=1,NM1
1363         JM1=J+M1
1364         DO I=1,NM1
1365            E1(I,J)=FMAS(I,J)*FAC1-FJAC(I,JM1)
1366         END DO
1367      END DO
1368      GOTO 45
1369C
1370C -----------------------------------------------------------
1371C
1372   6  CONTINUE
1373C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
1374C ---  THIS OPTION IS NOT PROVIDED
1375      RETURN
1376C
1377C -----------------------------------------------------------
1378C
1379   7  CONTINUE
1380C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
1381      IF (CALHES) CALL Elmhes (LDJAC,N,1,N,FJAC,IPHES)
1382      CALHES=.FALSE.
1383      DO J=1,N-1
1384         J1=J+1
1385         E1(J1,J)=-FJAC(J1,J)
1386      END DO
1387      DO J=1,N
1388         DO I=1,J
1389            E1(I,J)=-FJAC(I,J)
1390         END DO
1391         E1(J,J)=E1(J,J)+FAC1
1392      END DO
1393      CALL DECH(N,LDE1,E1,1,IP1,IER)
1394      RETURN
1395C
1396C -----------------------------------------------------------
1397C
1398  55  CONTINUE
1399      RETURN
1400      END
1401C
1402C     END OF SUBROUTINE DECOMR
1403C
1404C ***********************************************************
1405C
1406      SUBROUTINE DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
1407     &            M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB)
1408      IMPLICIT KPP_REAL (A-H,O-Z)
1409      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),
1410     &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1)
1411      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1412C
1413      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
1414C
1415C -----------------------------------------------------------
1416C
1417   1  CONTINUE
1418C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1419      DO J=1,N
1420         DO I=1,N
1421            E2R(I,J)=-FJAC(I,J)
1422            E2I(I,J)=0.D0
1423         END DO
1424         E2R(J,J)=E2R(J,J)+ALPHN
1425         E2I(J,J)=BETAN
1426      END DO
1427      CALL DECC (N,LDE1,E2R,E2I,IP2,IER)
1428      RETURN
1429C
1430C -----------------------------------------------------------
1431C
1432  11  CONTINUE
1433C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1434      DO J=1,NM1
1435         JM1=J+M1
1436         DO I=1,NM1
1437            E2R(I,J)=-FJAC(I,JM1)
1438            E2I(I,J)=0.D0
1439         END DO
1440         E2R(J,J)=E2R(J,J)+ALPHN
1441         E2I(J,J)=BETAN
1442      END DO
1443  45  MM=M1/M2
1444      ABNO=ALPHN**2+BETAN**2
1445      ALP=ALPHN/ABNO
1446      BET=BETAN/ABNO
1447      DO J=1,M2
1448         DO I=1,NM1
1449            SUMR=0.D0
1450            SUMI=0.D0
1451            DO K=0,MM-1
1452               SUMS=SUMR+FJAC(I,J+K*M2)
1453               SUMR=SUMS*ALP+SUMI*BET
1454               SUMI=SUMI*ALP-SUMS*BET
1455            END DO
1456            E2R(I,J)=E2R(I,J)-SUMR
1457            E2I(I,J)=E2I(I,J)-SUMI
1458         END DO
1459      END DO
1460      CALL DECC (NM1,LDE1,E2R,E2I,IP2,IER)
1461      RETURN
1462C
1463C -----------------------------------------------------------
1464C
1465   2  CONTINUE
1466C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1467      DO J=1,N
1468         DO I=1,MBJAC
1469            IMLE=I+MLE
1470            E2R(IMLE,J)=-FJAC(I,J)
1471            E2I(IMLE,J)=0.D0
1472         END DO
1473         E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN
1474         E2I(MDIAG,J)=BETAN
1475      END DO
1476      CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
1477      RETURN
1478C
1479C -----------------------------------------------------------
1480C
1481  12  CONTINUE
1482C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1483      DO J=1,NM1
1484         JM1=J+M1
1485         DO I=1,MBJAC
1486            E2R(I+MLE,J)=-FJAC(I,JM1)
1487            E2I(I+MLE,J)=0.D0
1488         END DO
1489         E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN
1490         E2I(MDIAG,J)=E2I(MDIAG,J)+BETAN
1491      END DO
1492  46  MM=M1/M2
1493      ABNO=ALPHN**2+BETAN**2
1494      ALP=ALPHN/ABNO
1495      BET=BETAN/ABNO
1496      DO J=1,M2
1497         DO I=1,MBJAC
1498            SUMR=0.D0
1499            SUMI=0.D0
1500            DO K=0,MM-1
1501               SUMS=SUMR+FJAC(I,J+K*M2)
1502               SUMR=SUMS*ALP+SUMI*BET
1503               SUMI=SUMI*ALP-SUMS*BET
1504            END DO
1505            IMLE=I+MLE
1506            E2R(IMLE,J)=E2R(IMLE,J)-SUMR
1507            E2I(IMLE,J)=E2I(IMLE,J)-SUMI
1508         END DO
1509      END DO
1510      CALL DECBC (NM1,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
1511      RETURN
1512C
1513C -----------------------------------------------------------
1514C
1515   3  CONTINUE
1516C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1517      DO  J=1,N
1518         DO  I=1,N
1519            E2R(I,J)=-FJAC(I,J)
1520            E2I(I,J)=0.D0
1521         END DO
1522      END DO
1523      DO J=1,N
1524         DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS)
1525            BB=FMAS(I-J+MBDIAG,J)
1526            E2R(I,J)=E2R(I,J)+ALPHN*BB
1527            E2I(I,J)=BETAN*BB
1528         END DO
1529      END DO
1530      CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
1531      RETURN
1532C
1533C -----------------------------------------------------------
1534C
1535  13  CONTINUE
1536C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1537      DO J=1,NM1
1538         JM1=J+M1
1539         DO I=1,NM1
1540            E2R(I,J)=-FJAC(I,JM1)
1541            E2I(I,J)=0.D0
1542         END DO
1543         DO I=MAX(1,J-MUMAS),MIN(NM1,J+MLMAS)
1544            FFMA=FMAS(I-J+MBDIAG,J)
1545            E2R(I,J)=E2R(I,J)+ALPHN*FFMA
1546            E2I(I,J)=E2I(I,J)+BETAN*FFMA
1547         END DO
1548      END DO
1549      GOTO 45
1550C
1551C -----------------------------------------------------------
1552C
1553   4  CONTINUE
1554C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
1555      DO J=1,N
1556         DO I=1,MBJAC
1557            IMLE=I+MLE
1558            E2R(IMLE,J)=-FJAC(I,J)
1559            E2I(IMLE,J)=0.D0
1560         END DO
1561         DO I=MAX(1,MUMAS+2-J),MIN(MBB,MUMAS+1-J+N)
1562            IB=I+MDIFF
1563            BB=FMAS(I,J)
1564            E2R(IB,J)=E2R(IB,J)+ALPHN*BB
1565            E2I(IB,J)=BETAN*BB
1566         END DO
1567      END DO
1568      CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
1569      RETURN
1570C
1571C -----------------------------------------------------------
1572C
1573  14  CONTINUE
1574C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
1575      DO J=1,NM1
1576         JM1=J+M1
1577         DO I=1,MBJAC
1578            E2R(I+MLE,J)=-FJAC(I,JM1)
1579            E2I(I+MLE,J)=0.D0
1580         END DO
1581         DO I=1,MBB
1582            IB=I+MDIFF
1583            FFMA=FMAS(I,J)
1584            E2R(IB,J)=E2R(IB,J)+ALPHN*FFMA
1585            E2I(IB,J)=E2I(IB,J)+BETAN*FFMA
1586         END DO
1587      END DO
1588      GOTO 46
1589C
1590C -----------------------------------------------------------
1591C
1592   5  CONTINUE
1593C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
1594      DO J=1,N
1595         DO I=1,N
1596            BB=FMAS(I,J)
1597            E2R(I,J)=BB*ALPHN-FJAC(I,J)
1598            E2I(I,J)=BB*BETAN
1599         END DO
1600      END DO
1601      CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
1602      RETURN
1603C
1604C -----------------------------------------------------------
1605C
1606  15  CONTINUE
1607C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1608      DO J=1,NM1
1609         JM1=J+M1
1610         DO I=1,NM1
1611            E2R(I,J)=ALPHN*FMAS(I,J)-FJAC(I,JM1)
1612            E2I(I,J)=BETAN*FMAS(I,J)
1613         END DO
1614      END DO
1615      GOTO 45
1616C
1617C -----------------------------------------------------------
1618C
1619   6  CONTINUE
1620C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
1621C ---  THIS OPTION IS NOT PROVIDED
1622      RETURN
1623C
1624C -----------------------------------------------------------
1625C
1626   7  CONTINUE
1627C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
1628      DO J=1,N-1
1629         J1=J+1
1630         E2R(J1,J)=-FJAC(J1,J)
1631         E2I(J1,J)=0.D0
1632      END DO
1633      DO J=1,N
1634         DO I=1,J
1635            E2I(I,J)=0.D0
1636            E2R(I,J)=-FJAC(I,J)
1637         END DO
1638         E2R(J,J)=E2R(J,J)+ALPHN
1639         E2I(J,J)=BETAN
1640      END DO
1641      CALL DECHC(N,LDE1,E2R,E2I,1,IP2,IER)
1642      RETURN
1643C
1644C -----------------------------------------------------------
1645C
1646  55  CONTINUE
1647      RETURN
1648      END
1649C
1650C     END OF SUBROUTINE DECOMC
1651C
1652C ***********************************************************
1653C
1654      SUBROUTINE SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1655     &          M1,M2,NM1,FAC1,E1,LDE1,Z1,F1,IP1,IPHES,IER,IJOB)
1656      IMPLICIT KPP_REAL (A-H,O-Z)
1657      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
1658     &          IP1(NM1),IPHES(N),Z1(N),F1(N)
1659      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1660C
1661      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB
1662C
1663C -----------------------------------------------------------
1664C
1665   1  CONTINUE
1666C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1667      DO I=1,N
1668         Z1(I)=Z1(I)-F1(I)*FAC1
1669      END DO
1670      CALL SOL (N,LDE1,E1,Z1,IP1)
1671      RETURN
1672C
1673C -----------------------------------------------------------
1674C
1675  11  CONTINUE
1676C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1677      DO I=1,N
1678         Z1(I)=Z1(I)-F1(I)*FAC1
1679      END DO
1680 48   CONTINUE
1681      MM=M1/M2
1682      DO J=1,M2
1683         SUM1=0.D0
1684         DO K=MM-1,0,-1
1685            JKM=J+K*M2
1686            SUM1=(Z1(JKM)+SUM1)/FAC1
1687            DO I=1,NM1
1688               IM1=I+M1
1689               Z1(IM1)=Z1(IM1)+FJAC(I,JKM)*SUM1
1690            END DO
1691         END DO
1692      END DO
1693      CALL SOL (NM1,LDE1,E1,Z1(M1+1),IP1)
1694 49   CONTINUE
1695      DO I=M1,1,-1
1696         Z1(I)=(Z1(I)+Z1(M2+I))/FAC1
1697      END DO
1698      RETURN
1699C
1700C -----------------------------------------------------------
1701C
1702   2  CONTINUE
1703C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1704      DO I=1,N
1705         Z1(I)=Z1(I)-F1(I)*FAC1
1706      END DO
1707      CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
1708      RETURN
1709C
1710C -----------------------------------------------------------
1711C
1712  12  CONTINUE
1713C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1714      DO I=1,N
1715         Z1(I)=Z1(I)-F1(I)*FAC1
1716      END DO
1717  45  CONTINUE
1718      MM=M1/M2
1719      DO J=1,M2
1720         SUM1=0.D0
1721         DO K=MM-1,0,-1
1722            JKM=J+K*M2
1723            SUM1=(Z1(JKM)+SUM1)/FAC1
1724            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
1725               IM1=I+M1
1726               Z1(IM1)=Z1(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM1
1727            END DO
1728         END DO
1729      END DO
1730      CALL SOLB (NM1,LDE1,E1,MLE,MUE,Z1(M1+1),IP1)
1731      GOTO 49
1732C
1733C -----------------------------------------------------------
1734C
1735   3  CONTINUE
1736C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1737      DO I=1,N
1738         S1=0.0D0
1739         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
1740            S1=S1-FMAS(I-J+MBDIAG,J)*F1(J)
1741         END DO
1742         Z1(I)=Z1(I)+S1*FAC1
1743      END DO
1744      CALL SOL (N,LDE1,E1,Z1,IP1)
1745      RETURN
1746C
1747C -----------------------------------------------------------
1748C
1749  13  CONTINUE
1750C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1751      DO I=1,M1
1752         Z1(I)=Z1(I)-F1(I)*FAC1
1753      END DO
1754      DO I=1,NM1
1755         IM1=I+M1
1756         S1=0.0D0
1757         DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
1758            S1=S1-FMAS(I-J+MBDIAG,J)*F1(J+M1)
1759         END DO
1760         Z1(IM1)=Z1(IM1)+S1*FAC1
1761      END DO
1762      IF (IJOB.EQ.14) GOTO 45
1763      GOTO 48
1764C
1765C -----------------------------------------------------------
1766C
1767   4  CONTINUE
1768C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
1769      DO I=1,N
1770         S1=0.0D0
1771         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
1772            S1=S1-FMAS(I-J+MBDIAG,J)*F1(J)
1773         END DO
1774         Z1(I)=Z1(I)+S1*FAC1
1775      END DO
1776      CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
1777      RETURN
1778C
1779C -----------------------------------------------------------
1780C
1781   5  CONTINUE
1782C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
1783      DO I=1,N
1784         S1=0.0D0
1785         DO J=1,N
1786            S1=S1-FMAS(I,J)*F1(J)
1787         END DO
1788         Z1(I)=Z1(I)+S1*FAC1
1789      END DO
1790      CALL SOL (N,LDE1,E1,Z1,IP1)
1791      RETURN
1792C
1793C -----------------------------------------------------------
1794C
1795  15  CONTINUE
1796C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1797      DO I=1,M1
1798         Z1(I)=Z1(I)-F1(I)*FAC1
1799      END DO
1800      DO I=1,NM1
1801         IM1=I+M1
1802         S1=0.0D0
1803         DO J=1,NM1
1804            S1=S1-FMAS(I,J)*F1(J+M1)
1805         END DO
1806         Z1(IM1)=Z1(IM1)+S1*FAC1
1807      END DO
1808      GOTO 48
1809C
1810C -----------------------------------------------------------
1811C
1812   6  CONTINUE
1813C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
1814C ---  THIS OPTION IS NOT PROVIDED
1815      RETURN
1816C
1817C -----------------------------------------------------------
1818C
1819   7  CONTINUE
1820C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
1821      DO I=1,N
1822         Z1(I)=Z1(I)-F1(I)*FAC1
1823      END DO
1824      DO MM=N-2,1,-1
1825          MP=N-MM
1826          MP1=MP-1
1827          I=IPHES(MP)
1828          IF (I.EQ.MP) GOTO 746
1829          ZSAFE=Z1(MP)
1830          Z1(MP)=Z1(I)
1831          Z1(I)=ZSAFE
1832 746      CONTINUE
1833          DO I=MP+1,N
1834             Z1(I)=Z1(I)-FJAC(I,MP1)*Z1(MP)
1835          END DO
1836       END DO
1837       CALL SOLH(N,LDE1,E1,1,Z1,IP1)
1838       DO MM=1,N-2
1839          MP=N-MM
1840          MP1=MP-1
1841          DO I=MP+1,N
1842             Z1(I)=Z1(I)+FJAC(I,MP1)*Z1(MP)
1843          END DO
1844          I=IPHES(MP)
1845          IF (I.EQ.MP) GOTO 750
1846          ZSAFE=Z1(MP)
1847          Z1(MP)=Z1(I)
1848          Z1(I)=ZSAFE
1849 750      CONTINUE
1850      END DO
1851      RETURN
1852C
1853C -----------------------------------------------------------
1854C
1855  55  CONTINUE
1856      RETURN
1857      END
1858C
1859C     END OF SUBROUTINE SLVRAR
1860C
1861C ***********************************************************
1862C
1863      SUBROUTINE SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1864     &          M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,Z2,Z3,
1865     &          F2,F3,CONT,IP2,IPHES,IER,IJOB)
1866      IMPLICIT KPP_REAL (A-H,O-Z)
1867      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),
1868     &          IP2(NM1),IPHES(N),Z2(N),Z3(N),F2(N),F3(N)
1869      DIMENSION E2R(LDE1,NM1),E2I(LDE1,NM1)
1870      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1871C
1872      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB
1873C
1874C -----------------------------------------------------------
1875C
1876   1  CONTINUE
1877C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1878      DO I=1,N
1879         S2=-F2(I)
1880         S3=-F3(I)
1881         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1882         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1883      END DO
1884      CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2)
1885      RETURN
1886C
1887C -----------------------------------------------------------
1888C
1889  11  CONTINUE
1890C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1891      DO I=1,N
1892         S2=-F2(I)
1893         S3=-F3(I)
1894         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1895         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1896      END DO
1897 48   ABNO=ALPHN**2+BETAN**2
1898      MM=M1/M2
1899      DO J=1,M2
1900         SUM2=0.D0
1901         SUM3=0.D0
1902         DO K=MM-1,0,-1
1903            JKM=J+K*M2
1904            SUMH=(Z2(JKM)+SUM2)/ABNO
1905            SUM3=(Z3(JKM)+SUM3)/ABNO
1906            SUM2=SUMH*ALPHN+SUM3*BETAN
1907            SUM3=SUM3*ALPHN-SUMH*BETAN
1908            DO I=1,NM1
1909               IM1=I+M1
1910               Z2(IM1)=Z2(IM1)+FJAC(I,JKM)*SUM2
1911               Z3(IM1)=Z3(IM1)+FJAC(I,JKM)*SUM3
1912            END DO
1913         END DO
1914      END DO
1915      CALL SOLC (NM1,LDE1,E2R,E2I,Z2(M1+1),Z3(M1+1),IP2)
1916 49   CONTINUE
1917      DO I=M1,1,-1
1918         MPI=M2+I
1919         Z2I=Z2(I)+Z2(MPI)
1920         Z3I=Z3(I)+Z3(MPI)
1921         Z3(I)=(Z3I*ALPHN-Z2I*BETAN)/ABNO
1922         Z2(I)=(Z2I*ALPHN+Z3I*BETAN)/ABNO
1923      END DO
1924      RETURN
1925C
1926C -----------------------------------------------------------
1927C
1928   2  CONTINUE
1929C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1930      DO I=1,N
1931         S2=-F2(I)
1932         S3=-F3(I)
1933         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1934         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1935      END DO
1936      CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
1937      RETURN
1938C
1939C -----------------------------------------------------------
1940C
1941  12  CONTINUE
1942C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1943      DO I=1,N
1944         S2=-F2(I)
1945         S3=-F3(I)
1946         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1947         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1948      END DO
1949  45  ABNO=ALPHN**2+BETAN**2
1950      MM=M1/M2
1951      DO J=1,M2
1952         SUM2=0.D0
1953         SUM3=0.D0
1954         DO K=MM-1,0,-1
1955            JKM=J+K*M2
1956            SUMH=(Z2(JKM)+SUM2)/ABNO
1957            SUM3=(Z3(JKM)+SUM3)/ABNO
1958            SUM2=SUMH*ALPHN+SUM3*BETAN
1959            SUM3=SUM3*ALPHN-SUMH*BETAN
1960            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
1961               IM1=I+M1
1962               IIMU=I+MUJAC+1-J
1963               Z2(IM1)=Z2(IM1)+FJAC(IIMU,JKM)*SUM2
1964               Z3(IM1)=Z3(IM1)+FJAC(IIMU,JKM)*SUM3
1965            END DO
1966         END DO
1967      END DO
1968      CALL SOLBC (NM1,LDE1,E2R,E2I,MLE,MUE,Z2(M1+1),Z3(M1+1),IP2)
1969      GOTO 49
1970C
1971C -----------------------------------------------------------
1972C
1973   3  CONTINUE
1974C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1975      DO I=1,N
1976         S2=0.0D0
1977         S3=0.0D0
1978         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
1979            BB=FMAS(I-J+MBDIAG,J)
1980            S2=S2-BB*F2(J)
1981            S3=S3-BB*F3(J)
1982         END DO
1983         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1984         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1985      END DO
1986      CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
1987      RETURN
1988C
1989C -----------------------------------------------------------
1990C
1991  13  CONTINUE
1992C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1993      DO I=1,M1
1994         S2=-F2(I)
1995         S3=-F3(I)
1996         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1997         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1998      END DO
1999      DO I=1,NM1
2000         IM1=I+M1
2001         S2=0.0D0
2002         S3=0.0D0
2003         DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2004            JM1=J+M1
2005            BB=FMAS(I-J+MBDIAG,J)
2006            S2=S2-BB*F2(JM1)
2007            S3=S3-BB*F3(JM1)
2008         END DO
2009         Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2010         Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2011      END DO
2012      IF (IJOB.EQ.14) GOTO 45
2013      GOTO 48
2014C
2015C -----------------------------------------------------------
2016C
2017   4  CONTINUE
2018C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2019      DO I=1,N
2020         S2=0.0D0
2021         S3=0.0D0
2022         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2023            BB=FMAS(I-J+MBDIAG,J)
2024            S2=S2-BB*F2(J)
2025            S3=S3-BB*F3(J)
2026         END DO
2027         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2028         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2029      END DO
2030      CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
2031      RETURN
2032C
2033C -----------------------------------------------------------
2034C
2035   5  CONTINUE
2036C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2037      DO I=1,N
2038         S2=0.0D0
2039         S3=0.0D0
2040         DO J=1,N
2041            BB=FMAS(I,J)
2042            S2=S2-BB*F2(J)
2043            S3=S3-BB*F3(J)
2044         END DO
2045         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2046         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2047      END DO
2048      CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
2049      RETURN
2050C
2051C -----------------------------------------------------------
2052C
2053  15  CONTINUE
2054C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2055      DO I=1,M1
2056         S2=-F2(I)
2057         S3=-F3(I)
2058         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2059         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2060      END DO
2061      DO I=1,NM1
2062         IM1=I+M1
2063         S2=0.0D0
2064         S3=0.0D0
2065         DO J=1,NM1
2066            JM1=J+M1
2067            BB=FMAS(I,J)
2068            S2=S2-BB*F2(JM1)
2069            S3=S3-BB*F3(JM1)
2070         END DO
2071         Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2072         Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2073      END DO
2074      GOTO 48
2075C
2076C -----------------------------------------------------------
2077C
2078   6  CONTINUE
2079C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
2080C ---  THIS OPTION IS NOT PROVIDED
2081      RETURN
2082C
2083C -----------------------------------------------------------
2084C
2085   7  CONTINUE
2086C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
2087      DO I=1,N
2088         S2=-F2(I)
2089         S3=-F3(I)
2090         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2091         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2092      END DO
2093      DO MM=N-2,1,-1
2094          MP=N-MM
2095          MP1=MP-1
2096          I=IPHES(MP)
2097          IF (I.EQ.MP) GOTO 746
2098          ZSAFE=Z2(MP)
2099          Z2(MP)=Z2(I)
2100          Z2(I)=ZSAFE
2101          ZSAFE=Z3(MP)
2102          Z3(MP)=Z3(I)
2103          Z3(I)=ZSAFE
2104 746      CONTINUE
2105          DO I=MP+1,N
2106             E1IMP=FJAC(I,MP1)
2107             Z2(I)=Z2(I)-E1IMP*Z2(MP)
2108             Z3(I)=Z3(I)-E1IMP*Z3(MP)
2109          END DO
2110       END DO
2111       CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2)
2112       DO MM=1,N-2
2113          MP=N-MM
2114          MP1=MP-1
2115          DO I=MP+1,N
2116             E1IMP=FJAC(I,MP1)
2117             Z2(I)=Z2(I)+E1IMP*Z2(MP)
2118             Z3(I)=Z3(I)+E1IMP*Z3(MP)
2119          END DO
2120          I=IPHES(MP)
2121          IF (I.EQ.MP) GOTO 750
2122          ZSAFE=Z2(MP)
2123          Z2(MP)=Z2(I)
2124          Z2(I)=ZSAFE
2125          ZSAFE=Z3(MP)
2126          Z3(MP)=Z3(I)
2127          Z3(I)=ZSAFE
2128 750      CONTINUE
2129      END DO
2130      RETURN
2131C
2132C -----------------------------------------------------------
2133C
2134  55  CONTINUE
2135      RETURN
2136      END
2137C
2138C     END OF SUBROUTINE SLVRAI
2139C
2140C ***********************************************************
2141C
2142      SUBROUTINE SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
2143     &          M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3,
2144     &          F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB)
2145      IMPLICIT KPP_REAL (A-H,O-Z)
2146      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
2147     &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP1(NM1),IP2(NM1),
2148     &          IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N)
2149      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
2150C
2151      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB
2152C
2153C -----------------------------------------------------------
2154C
2155   1  CONTINUE
2156C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
2157      DO I=1,N
2158         S2=-F2(I)
2159         S3=-F3(I)
2160         Z1(I)=Z1(I)-F1(I)*FAC1
2161         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2162         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2163      END DO
2164      CALL SOL (N,LDE1,E1,Z1,IP1)
2165      CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2)
2166      RETURN
2167C
2168C -----------------------------------------------------------
2169C
2170  11  CONTINUE
2171C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
2172      DO I=1,N
2173         S2=-F2(I)
2174         S3=-F3(I)
2175         Z1(I)=Z1(I)-F1(I)*FAC1
2176         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2177         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2178      END DO
2179 48   ABNO=ALPHN**2+BETAN**2
2180      MM=M1/M2
2181      DO J=1,M2
2182         SUM1=0.D0
2183         SUM2=0.D0
2184         SUM3=0.D0
2185         DO K=MM-1,0,-1
2186            JKM=J+K*M2
2187            SUM1=(Z1(JKM)+SUM1)/FAC1
2188            SUMH=(Z2(JKM)+SUM2)/ABNO
2189            SUM3=(Z3(JKM)+SUM3)/ABNO
2190            SUM2=SUMH*ALPHN+SUM3*BETAN
2191            SUM3=SUM3*ALPHN-SUMH*BETAN
2192            DO I=1,NM1
2193               IM1=I+M1
2194               Z1(IM1)=Z1(IM1)+FJAC(I,JKM)*SUM1
2195               Z2(IM1)=Z2(IM1)+FJAC(I,JKM)*SUM2
2196               Z3(IM1)=Z3(IM1)+FJAC(I,JKM)*SUM3
2197            END DO
2198         END DO
2199      END DO
2200      CALL SOL (NM1,LDE1,E1,Z1(M1+1),IP1)
2201      CALL SOLC (NM1,LDE1,E2R,E2I,Z2(M1+1),Z3(M1+1),IP2)
2202 49   CONTINUE
2203      DO I=M1,1,-1
2204         MPI=M2+I
2205         Z1(I)=(Z1(I)+Z1(MPI))/FAC1
2206         Z2I=Z2(I)+Z2(MPI)
2207         Z3I=Z3(I)+Z3(MPI)
2208         Z3(I)=(Z3I*ALPHN-Z2I*BETAN)/ABNO
2209         Z2(I)=(Z2I*ALPHN+Z3I*BETAN)/ABNO
2210      END DO
2211      RETURN
2212C
2213C -----------------------------------------------------------
2214C
2215   2  CONTINUE
2216C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
2217      DO I=1,N
2218         S2=-F2(I)
2219         S3=-F3(I)
2220         Z1(I)=Z1(I)-F1(I)*FAC1
2221         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2222         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2223      END DO
2224      CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
2225      CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
2226      RETURN
2227C
2228C -----------------------------------------------------------
2229C
2230  12  CONTINUE
2231C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
2232      DO I=1,N
2233         S2=-F2(I)
2234         S3=-F3(I)
2235         Z1(I)=Z1(I)-F1(I)*FAC1
2236         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2237         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2238      END DO
2239  45  ABNO=ALPHN**2+BETAN**2
2240      MM=M1/M2
2241      DO J=1,M2
2242         SUM1=0.D0
2243         SUM2=0.D0
2244         SUM3=0.D0
2245         DO K=MM-1,0,-1
2246            JKM=J+K*M2
2247            SUM1=(Z1(JKM)+SUM1)/FAC1
2248            SUMH=(Z2(JKM)+SUM2)/ABNO
2249            SUM3=(Z3(JKM)+SUM3)/ABNO
2250            SUM2=SUMH*ALPHN+SUM3*BETAN
2251            SUM3=SUM3*ALPHN-SUMH*BETAN
2252            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2253               IM1=I+M1
2254               FFJA=FJAC(I+MUJAC+1-J,JKM)
2255               Z1(IM1)=Z1(IM1)+FFJA*SUM1
2256               Z2(IM1)=Z2(IM1)+FFJA*SUM2
2257               Z3(IM1)=Z3(IM1)+FFJA*SUM3
2258            END DO
2259         END DO
2260      END DO
2261      CALL SOLB (NM1,LDE1,E1,MLE,MUE,Z1(M1+1),IP1)
2262      CALL SOLBC (NM1,LDE1,E2R,E2I,MLE,MUE,Z2(M1+1),Z3(M1+1),IP2)
2263      GOTO 49
2264C
2265C -----------------------------------------------------------
2266C
2267   3  CONTINUE
2268C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
2269      DO I=1,N
2270         S1=0.0D0
2271         S2=0.0D0
2272         S3=0.0D0
2273         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2274            BB=FMAS(I-J+MBDIAG,J)
2275            S1=S1-BB*F1(J)
2276            S2=S2-BB*F2(J)
2277            S3=S3-BB*F3(J)
2278         END DO
2279         Z1(I)=Z1(I)+S1*FAC1
2280         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2281         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2282      END DO
2283      CALL SOL (N,LDE1,E1,Z1,IP1)
2284      CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
2285      RETURN
2286C
2287C -----------------------------------------------------------
2288C
2289  13  CONTINUE
2290C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2291      DO I=1,M1
2292         S2=-F2(I)
2293         S3=-F3(I)
2294         Z1(I)=Z1(I)-F1(I)*FAC1
2295         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2296         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2297      END DO
2298      DO I=1,NM1
2299         IM1=I+M1
2300         S1=0.0D0
2301         S2=0.0D0
2302         S3=0.0D0
2303         J1B=MAX(1,I-MLMAS)
2304         J2B=MIN(NM1,I+MUMAS)
2305         DO J=J1B,J2B
2306            JM1=J+M1
2307            BB=FMAS(I-J+MBDIAG,J)
2308            S1=S1-BB*F1(JM1)
2309            S2=S2-BB*F2(JM1)
2310            S3=S3-BB*F3(JM1)
2311         END DO
2312         Z1(IM1)=Z1(IM1)+S1*FAC1
2313         Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2314         Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2315      END DO
2316      IF (IJOB.EQ.14) GOTO 45
2317      GOTO 48
2318C
2319C -----------------------------------------------------------
2320C
2321   4  CONTINUE
2322C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2323      DO I=1,N
2324         S1=0.0D0
2325         S2=0.0D0
2326         S3=0.0D0
2327         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2328            BB=FMAS(I-J+MBDIAG,J)
2329            S1=S1-BB*F1(J)
2330            S2=S2-BB*F2(J)
2331            S3=S3-BB*F3(J)
2332         END DO
2333         Z1(I)=Z1(I)+S1*FAC1
2334         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2335         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2336      END DO
2337      CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
2338      CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
2339      RETURN
2340C
2341C -----------------------------------------------------------
2342C
2343   5  CONTINUE
2344C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2345      DO I=1,N
2346         S1=0.0D0
2347         S2=0.0D0
2348         S3=0.0D0
2349         DO J=1,N
2350            BB=FMAS(I,J)
2351            S1=S1-BB*F1(J)
2352            S2=S2-BB*F2(J)
2353            S3=S3-BB*F3(J)
2354         END DO
2355         Z1(I)=Z1(I)+S1*FAC1
2356         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2357         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2358      END DO
2359      CALL SOL (N,LDE1,E1,Z1,IP1)
2360      CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
2361      RETURN
2362C
2363C -----------------------------------------------------------
2364C
2365  15  CONTINUE
2366C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2367      DO I=1,M1
2368         S2=-F2(I)
2369         S3=-F3(I)
2370         Z1(I)=Z1(I)-F1(I)*FAC1
2371         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2372         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2373      END DO
2374      DO I=1,NM1
2375         IM1=I+M1
2376         S1=0.0D0
2377         S2=0.0D0
2378         S3=0.0D0
2379         DO J=1,NM1
2380            JM1=J+M1
2381            BB=FMAS(I,J)
2382            S1=S1-BB*F1(JM1)
2383            S2=S2-BB*F2(JM1)
2384            S3=S3-BB*F3(JM1)
2385         END DO
2386         Z1(IM1)=Z1(IM1)+S1*FAC1
2387         Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2388         Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2389      END DO
2390      GOTO 48
2391C
2392C -----------------------------------------------------------
2393C
2394   6  CONTINUE
2395C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
2396C ---  THIS OPTION IS NOT PROVIDED
2397      RETURN
2398C
2399C -----------------------------------------------------------
2400C
2401   7  CONTINUE
2402C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
2403      DO I=1,N
2404         S2=-F2(I)
2405         S3=-F3(I)
2406         Z1(I)=Z1(I)-F1(I)*FAC1
2407         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2408         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2409      END DO
2410      DO MM=N-2,1,-1
2411          MP=N-MM
2412          MP1=MP-1
2413          I=IPHES(MP)
2414          IF (I.EQ.MP) GOTO 746
2415          ZSAFE=Z1(MP)
2416          Z1(MP)=Z1(I)
2417          Z1(I)=ZSAFE
2418          ZSAFE=Z2(MP)
2419          Z2(MP)=Z2(I)
2420          Z2(I)=ZSAFE
2421          ZSAFE=Z3(MP)
2422          Z3(MP)=Z3(I)
2423          Z3(I)=ZSAFE
2424 746      CONTINUE
2425          DO I=MP+1,N
2426             E1IMP=FJAC(I,MP1)
2427             Z1(I)=Z1(I)-E1IMP*Z1(MP)
2428             Z2(I)=Z2(I)-E1IMP*Z2(MP)
2429             Z3(I)=Z3(I)-E1IMP*Z3(MP)
2430          END DO
2431       END DO
2432       CALL SOLH(N,LDE1,E1,1,Z1,IP1)
2433       CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2)
2434       DO MM=1,N-2
2435          MP=N-MM
2436          MP1=MP-1
2437          DO I=MP+1,N
2438             E1IMP=FJAC(I,MP1)
2439             Z1(I)=Z1(I)+E1IMP*Z1(MP)
2440             Z2(I)=Z2(I)+E1IMP*Z2(MP)
2441             Z3(I)=Z3(I)+E1IMP*Z3(MP)
2442          END DO
2443          I=IPHES(MP)
2444          IF (I.EQ.MP) GOTO 750
2445          ZSAFE=Z1(MP)
2446          Z1(MP)=Z1(I)
2447          Z1(I)=ZSAFE
2448          ZSAFE=Z2(MP)
2449          Z2(MP)=Z2(I)
2450          Z2(I)=ZSAFE
2451          ZSAFE=Z3(MP)
2452          Z3(MP)=Z3(I)
2453          Z3(I)=ZSAFE
2454 750      CONTINUE
2455      END DO
2456      RETURN
2457C
2458C -----------------------------------------------------------
2459C
2460  55  CONTINUE
2461      RETURN
2462      END
2463C
2464C     END OF SUBROUTINE SLVRAD
2465C
2466C ***********************************************************
2467C
2468      SUBROUTINE ESTRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
2469     &          H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,
2470     &          E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR,
2471     &          FIRST,REJECT,FAC1,RPAR,IPAR)
2472      IMPLICIT KPP_REAL (A-H,O-Z)
2473      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1),
2474     &     SCAL(N),IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),Y0(N),Y(N)
2475      DIMENSION CONT(N),RPAR(1),IPAR(1)
2476      LOGICAL FIRST,REJECT
2477      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
2478      EXTERNAL FCN
2479      HEE1=DD1/H
2480      HEE2=DD2/H
2481      HEE3=DD3/H
2482      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
2483C
2484   1  CONTINUE
2485C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
2486      DO  I=1,N
2487         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2488         CONT(I)=F2(I)+Y0(I)
2489      END DO
2490      CALL SOL (N,LDE1,E1,CONT,IP1)
2491      GOTO 77
2492C
2493  11  CONTINUE
2494C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
2495      DO I=1,N
2496         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2497         CONT(I)=F2(I)+Y0(I)
2498      END DO
2499  48  MM=M1/M2
2500      DO J=1,M2
2501         SUM1=0.D0
2502         DO K=MM-1,0,-1
2503            SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2504            DO I=1,NM1
2505               IM1=I+M1
2506               CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
2507            END DO
2508         END DO
2509      END DO
2510      CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1)
2511      DO I=M1,1,-1
2512         CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2513      END DO
2514      GOTO 77
2515C
2516   2  CONTINUE
2517C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX
2518      DO I=1,N
2519         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2520         CONT(I)=F2(I)+Y0(I)
2521      END DO
2522      CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2523      GOTO 77
2524C
2525  12  CONTINUE
2526C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
2527      DO I=1,N
2528         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2529         CONT(I)=F2(I)+Y0(I)
2530      END DO
2531  45  MM=M1/M2
2532      DO J=1,M2
2533         SUM1=0.D0
2534         DO K=MM-1,0,-1
2535            SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2536            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2537               IM1=I+M1
2538               CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
2539            END DO
2540         END DO
2541      END DO
2542      CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
2543      DO I=M1,1,-1
2544         CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2545      END DO
2546      GOTO 77
2547C
2548   3  CONTINUE
2549C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
2550      DO I=1,N
2551         F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2552      END DO
2553      DO I=1,N
2554         SUM=0.D0
2555         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2556            SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J)
2557         END DO
2558         F2(I)=SUM
2559         CONT(I)=SUM+Y0(I)
2560      END DO
2561      CALL SOL (N,LDE1,E1,CONT,IP1)
2562      GOTO 77
2563C
2564  13  CONTINUE
2565C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2566      DO I=1,M1
2567         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2568         CONT(I)=F2(I)+Y0(I)
2569      END DO
2570      DO I=M1+1,N
2571         F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2572      END DO
2573      DO I=1,NM1
2574         SUM=0.D0
2575         DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2576            SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J+M1)
2577         END DO
2578         IM1=I+M1
2579         F2(IM1)=SUM
2580         CONT(IM1)=SUM+Y0(IM1)
2581      END DO
2582      GOTO 48
2583C
2584   4  CONTINUE
2585C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2586      DO I=1,N
2587         F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2588      END DO
2589      DO I=1,N
2590         SUM=0.D0
2591         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2592            SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J)
2593         END DO
2594         F2(I)=SUM
2595         CONT(I)=SUM+Y0(I)
2596      END DO
2597      CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2598      GOTO 77
2599C
2600  14  CONTINUE
2601C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
2602      DO I=1,M1
2603         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2604         CONT(I)=F2(I)+Y0(I)
2605      END DO
2606      DO I=M1+1,N
2607         F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2608      END DO
2609      DO I=1,NM1
2610         SUM=0.D0
2611         DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2612            SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J+M1)
2613         END DO
2614         IM1=I+M1
2615         F2(IM1)=SUM
2616         CONT(IM1)=SUM+Y0(IM1)
2617      END DO
2618      GOTO 45
2619C
2620   5  CONTINUE
2621C ------  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2622      DO I=1,N
2623         F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2624      END DO
2625      DO I=1,N
2626         SUM=0.D0
2627         DO J=1,N
2628            SUM=SUM+FMAS(I,J)*F1(J)
2629         END DO
2630         F2(I)=SUM
2631         CONT(I)=SUM+Y0(I)
2632      END DO
2633      CALL SOL (N,LDE1,E1,CONT,IP1)
2634      GOTO 77
2635C
2636  15  CONTINUE
2637C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2638      DO I=1,M1
2639         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2640         CONT(I)=F2(I)+Y0(I)
2641      END DO
2642      DO I=M1+1,N
2643         F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2644      END DO
2645      DO I=1,NM1
2646         SUM=0.D0
2647         DO J=1,NM1
2648            SUM=SUM+FMAS(I,J)*F1(J+M1)
2649         END DO
2650         IM1=I+M1
2651         F2(IM1)=SUM
2652         CONT(IM1)=SUM+Y0(IM1)
2653      END DO
2654      GOTO 48
2655C
2656   6  CONTINUE
2657C ------  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
2658C ------  THIS OPTION IS NOT PROVIDED
2659      RETURN
2660C
2661   7  CONTINUE
2662C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
2663      DO I=1,N
2664         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2665         CONT(I)=F2(I)+Y0(I)
2666      END DO
2667      DO MM=N-2,1,-1
2668         MP=N-MM
2669         I=IPHES(MP)
2670         IF (I.EQ.MP) GOTO 310
2671         ZSAFE=CONT(MP)
2672         CONT(MP)=CONT(I)
2673         CONT(I)=ZSAFE
2674 310     CONTINUE
2675         DO I=MP+1,N
2676            CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
2677         END DO
2678      END DO
2679      CALL SOLH(N,LDE1,E1,1,CONT,IP1)
2680      DO MM=1,N-2
2681         MP=N-MM
2682         DO I=MP+1,N
2683            CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
2684         END DO
2685         I=IPHES(MP)
2686         IF (I.EQ.MP) GOTO 440
2687         ZSAFE=CONT(MP)
2688         CONT(MP)=CONT(I)
2689         CONT(I)=ZSAFE
2690 440     CONTINUE
2691      END DO
2692C
2693C --------------------------------------
2694C
2695  77  CONTINUE
2696      ERR=0.D0
2697      DO  I=1,N
2698         ERR=ERR+(CONT(I)/SCAL(I))**2
2699      END DO
2700      ERR=MAX(SQRT(ERR/N),1.D-10)
2701C
2702      IF (ERR.LT.1.D0) RETURN
2703      IF (FIRST.OR.REJECT) THEN
2704          DO I=1,N
2705             CONT(I)=Y(I)+CONT(I)
2706          END DO
2707          CALL FCN(N,X,CONT,F1,RPAR,IPAR)
2708          NFCN=NFCN+1
2709          DO I=1,N
2710             CONT(I)=F1(I)+F2(I)
2711          END DO
2712          GOTO (31,32,31,32,31,32,33,55,55,55,41,42,41,42,41), IJOB
2713C ------ FULL MATRIX OPTION
2714  31      CONTINUE
2715          CALL SOL(N,LDE1,E1,CONT,IP1)
2716          GOTO 88
2717C ------ FULL MATRIX OPTION, SECOND ORDER
2718 41      CONTINUE
2719         DO J=1,M2
2720            SUM1=0.D0
2721            DO K=MM-1,0,-1
2722               SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2723               DO I=1,NM1
2724                  IM1=I+M1
2725                  CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
2726               END DO
2727            END DO
2728         END DO
2729         CALL SOL(NM1,LDE1,E1,CONT(M1+1),IP1)
2730         DO I=M1,1,-1
2731            CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2732         END DO
2733         GOTO 88
2734C ------ BANDED MATRIX OPTION
2735 32      CONTINUE
2736         CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2737         GOTO 88
2738C ------ BANDED MATRIX OPTION, SECOND ORDER
2739 42      CONTINUE
2740         DO J=1,M2
2741            SUM1=0.D0
2742            DO K=MM-1,0,-1
2743               SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2744               DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2745                  IM1=I+M1
2746                  CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
2747               END DO
2748            END DO
2749         END DO
2750         CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
2751         DO I=M1,1,-1
2752            CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2753         END DO
2754          GOTO 88
2755C ------ HESSENBERG MATRIX OPTION
2756  33      CONTINUE
2757          DO MM=N-2,1,-1
2758             MP=N-MM
2759             I=IPHES(MP)
2760             IF (I.EQ.MP) GOTO 510
2761             ZSAFE=CONT(MP)
2762             CONT(MP)=CONT(I)
2763             CONT(I)=ZSAFE
2764 510         CONTINUE
2765             DO I=MP+1,N
2766                CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
2767             END DO
2768          END DO
2769          CALL SOLH(N,LDE1,E1,1,CONT,IP1)
2770          DO MM=1,N-2
2771             MP=N-MM
2772             DO I=MP+1,N
2773                CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
2774             END DO
2775             I=IPHES(MP)
2776             IF (I.EQ.MP) GOTO 640
2777             ZSAFE=CONT(MP)
2778             CONT(MP)=CONT(I)
2779             CONT(I)=ZSAFE
2780 640         CONTINUE
2781          END DO
2782C -----------------------------------
2783   88     CONTINUE
2784          ERR=0.D0
2785          DO I=1,N
2786             ERR=ERR+(CONT(I)/SCAL(I))**2
2787          END DO
2788          ERR=MAX(SQRT(ERR/N),1.D-10)
2789       END IF
2790       RETURN
2791C -----------------------------------------------------------
2792  55   CONTINUE
2793       RETURN
2794       END
2795C
2796C     END OF SUBROUTINE ESTRAD
2797C
2798C ***********************************************************
2799C
2800      SUBROUTINE ESTRAV(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
2801     &          H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS,
2802     &          E1,LDE1,ZZ,CONT,FF,IP1,IPHES,SCAL,ERR,
2803     &          FIRST,REJECT,FAC1,RPAR,IPAR)
2804      IMPLICIT KPP_REAL (A-H,O-Z)
2805      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1),
2806     &     SCAL(N),IPHES(N),ZZ(NNS),FF(NNS),Y0(N),Y(N)
2807      DIMENSION DD(NS),CONT(N),RPAR(1),IPAR(1)
2808      LOGICAL FIRST,REJECT
2809      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
2810      EXTERNAL FCN
2811
2812      GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
2813C
2814   1  CONTINUE
2815C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
2816      DO  I=1,N
2817         SUM=0.D0
2818         DO K=1,NS
2819            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2820         END DO
2821         FF(I+N)=SUM/H
2822         CONT(I)=FF(I+N)+Y0(I)
2823      END DO
2824      CALL SOL (N,LDE1,E1,CONT,IP1)
2825      GOTO 77
2826C
2827  11  CONTINUE
2828C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
2829      DO  I=1,N
2830         SUM=0.D0
2831         DO K=1,NS
2832            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2833         END DO
2834         FF(I+N)=SUM/H
2835         CONT(I)=FF(I+N)+Y0(I)
2836      END DO
2837  48  MM=M1/M2
2838      DO J=1,M2
2839         SUM1=0.D0
2840         DO K=MM-1,0,-1
2841            SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2842            DO I=1,NM1
2843               IM1=I+M1
2844               CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
2845            END DO
2846         END DO
2847      END DO
2848      CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1)
2849      DO I=M1,1,-1
2850         CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2851      END DO
2852      GOTO 77
2853C
2854   2  CONTINUE
2855C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX
2856      DO  I=1,N
2857         SUM=0.D0
2858         DO K=1,NS
2859            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2860         END DO
2861         FF(I+N)=SUM/H
2862         CONT(I)=FF(I+N)+Y0(I)
2863      END DO
2864      CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2865      GOTO 77
2866C
2867  12  CONTINUE
2868C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
2869      DO  I=1,N
2870         SUM=0.D0
2871         DO K=1,NS
2872            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2873         END DO
2874         FF(I+N)=SUM/H
2875         CONT(I)=FF(I+N)+Y0(I)
2876      END DO
2877  45  MM=M1/M2
2878      DO J=1,M2
2879         SUM1=0.D0
2880         DO K=MM-1,0,-1
2881            SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2882            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2883               IM1=I+M1
2884               CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
2885            END DO
2886         END DO
2887      END DO
2888      CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
2889      DO I=M1,1,-1
2890         CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2891      END DO
2892      GOTO 77
2893C
2894   3  CONTINUE
2895C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
2896      DO  I=1,N
2897         SUM=0.D0
2898         DO K=1,NS
2899            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2900         END DO
2901         FF(I)=SUM/H
2902      END DO
2903      DO I=1,N
2904         SUM=0.D0
2905         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2906            SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J)
2907         END DO
2908         FF(I+N)=SUM
2909         CONT(I)=SUM+Y0(I)
2910      END DO
2911      CALL SOL (N,LDE1,E1,CONT,IP1)
2912      GOTO 77
2913C
2914  13  CONTINUE
2915C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2916      DO  I=1,M1
2917         SUM=0.D0
2918         DO K=1,NS
2919            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2920         END DO
2921         FF(I+N)=SUM/H
2922         CONT(I)=FF(I+N)+Y0(I)
2923      END DO
2924      DO I=M1+1,N
2925         SUM=0.D0
2926         DO K=1,NS
2927            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2928         END DO
2929         FF(I)=SUM/H
2930      END DO
2931      DO I=1,NM1
2932         SUM=0.D0
2933         DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2934            SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J+M1)
2935         END DO
2936         IM1=I+M1
2937         FF(IM1+N)=SUM
2938         CONT(IM1)=SUM+Y0(IM1)
2939      END DO
2940      GOTO 48
2941C
2942   4  CONTINUE
2943C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2944      DO  I=1,N
2945         SUM=0.D0
2946         DO K=1,NS
2947            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2948         END DO
2949         FF(I)=SUM/H
2950      END DO
2951      DO I=1,N
2952         SUM=0.D0
2953         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2954            SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J)
2955         END DO
2956         FF(I+N)=SUM
2957         CONT(I)=SUM+Y0(I)
2958      END DO
2959      CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2960      GOTO 77
2961C
2962  14  CONTINUE
2963C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
2964      DO  I=1,M1
2965         SUM=0.D0
2966         DO K=1,NS
2967            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2968         END DO
2969         FF(I+N)=SUM/H
2970         CONT(I)=FF(I+N)+Y0(I)
2971      END DO
2972      DO I=M1+1,N
2973         SUM=0.D0
2974         DO K=1,NS
2975            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2976         END DO
2977         FF(I)=SUM/H
2978      END DO
2979      DO I=1,NM1
2980         SUM=0.D0
2981         DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2982            SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J+M1)
2983         END DO
2984         IM1=I+M1
2985         FF(IM1+N)=SUM
2986         CONT(IM1)=SUM+Y0(IM1)
2987      END DO
2988      GOTO 45
2989C
2990   5  CONTINUE
2991C ------  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2992      DO I=1,N
2993         SUM=0.D0
2994         DO K=1,NS
2995            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2996         END DO
2997         FF(I)=SUM/H
2998      END DO
2999      DO I=1,N
3000         SUM=0.D0
3001         DO J=1,N
3002            SUM=SUM+FMAS(I,J)*FF(J)
3003         END DO
3004         FF(I+N)=SUM
3005         CONT(I)=SUM+Y0(I)
3006      END DO
3007      CALL SOL (N,LDE1,E1,CONT,IP1)
3008      GOTO 77
3009C
3010  15  CONTINUE
3011C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
3012      DO  I=1,M1
3013         SUM=0.D0
3014         DO K=1,NS
3015            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
3016         END DO
3017         FF(I+N)=SUM/H
3018         CONT(I)=FF(I+N)+Y0(I)
3019      END DO
3020      DO I=M1+1,N
3021         SUM=0.D0
3022         DO K=1,NS
3023            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
3024         END DO
3025         FF(I)=SUM/H
3026      END DO
3027      DO I=1,NM1
3028         SUM=0.D0
3029         DO J=1,NM1
3030            SUM=SUM+FMAS(I,J)*FF(J+M1)
3031         END DO
3032         IM1=I+M1
3033         FF(IM1+N)=SUM
3034         CONT(IM1)=SUM+Y0(IM1)
3035      END DO
3036      GOTO 48
3037C
3038   6  CONTINUE
3039C ------  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
3040C ------  THIS OPTION IS NOT PROVIDED
3041      RETURN
3042C
3043   7  CONTINUE
3044C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
3045      DO  I=1,N
3046         SUM=0.D0
3047         DO K=1,NS
3048            SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
3049         END DO
3050         FF(I+N)=SUM/H
3051         CONT(I)=FF(I+N)+Y0(I)
3052      END DO
3053      DO MM=N-2,1,-1
3054         MP=N-MM
3055         I=IPHES(MP)
3056         IF (I.EQ.MP) GOTO 310
3057         ZSAFE=CONT(MP)
3058         CONT(MP)=CONT(I)
3059         CONT(I)=ZSAFE
3060 310     CONTINUE
3061         DO I=MP+1,N
3062            CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
3063         END DO
3064      END DO
3065      CALL SOLH(N,LDE1,E1,1,CONT,IP1)
3066      DO MM=1,N-2
3067         MP=N-MM
3068         DO I=MP+1,N
3069            CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
3070         END DO
3071         I=IPHES(MP)
3072         IF (I.EQ.MP) GOTO 440
3073         ZSAFE=CONT(MP)
3074         CONT(MP)=CONT(I)
3075         CONT(I)=ZSAFE
3076 440     CONTINUE
3077      END DO
3078C
3079C --------------------------------------
3080C
3081  77  CONTINUE
3082      ERR=0.D0
3083      DO  I=1,N
3084         ERR=ERR+(CONT(I)/SCAL(I))**2
3085      END DO
3086      ERR=MAX(SQRT(ERR/N),1.D-10)
3087C
3088      IF (ERR.LT.1.D0) RETURN
3089      IF (FIRST.OR.REJECT) THEN
3090          DO I=1,N
3091             CONT(I)=Y(I)+CONT(I)
3092          END DO
3093          CALL FCN(N,X,CONT,FF,RPAR,IPAR)
3094          NFCN=NFCN+1
3095          DO I=1,N
3096             CONT(I)=FF(I)+FF(I+N)
3097          END DO
3098          GOTO (31,32,31,32,31,32,33,55,55,55,41,42,41,42,41), IJOB
3099C ------ FULL MATRIX OPTION
3100 31      CONTINUE
3101         CALL SOL (N,LDE1,E1,CONT,IP1)
3102          GOTO 88
3103C ------ FULL MATRIX OPTION, SECOND ORDER
3104 41      CONTINUE
3105         DO J=1,M2
3106            SUM1=0.D0
3107            DO K=MM-1,0,-1
3108               SUM1=(CONT(J+K*M2)+SUM1)/FAC1
3109               DO I=1,NM1
3110                  IM1=I+M1
3111                  CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
3112               END DO
3113            END DO
3114         END DO
3115         CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1)
3116         DO I=M1,1,-1
3117            CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
3118         END DO
3119          GOTO 88
3120C ------ BANDED MATRIX OPTION
3121 32      CONTINUE
3122         CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
3123          GOTO 88
3124C ------ BANDED MATRIX OPTION, SECOND ORDER
3125 42      CONTINUE
3126         DO J=1,M2
3127            SUM1=0.D0
3128            DO K=MM-1,0,-1
3129               SUM1=(CONT(J+K*M2)+SUM1)/FAC1
3130               DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
3131                  IM1=I+M1
3132                  CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
3133               END DO
3134            END DO
3135         END DO
3136         CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
3137         DO I=M1,1,-1
3138            CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
3139         END DO
3140          GOTO 88
3141C ------ HESSENBERG MATRIX OPTION
3142  33      CONTINUE
3143          DO MM=N-2,1,-1
3144             MP=N-MM
3145             I=IPHES(MP)
3146             IF (I.EQ.MP) GOTO 510
3147             ZSAFE=CONT(MP)
3148             CONT(MP)=CONT(I)
3149             CONT(I)=ZSAFE
3150 510         CONTINUE
3151             DO I=MP+1,N
3152                CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
3153             END DO
3154          END DO
3155          CALL SOLH(N,LDE1,E1,1,CONT,IP1)
3156          DO MM=1,N-2
3157             MP=N-MM
3158             DO I=MP+1,N
3159                CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
3160             END DO
3161             I=IPHES(MP)
3162             IF (I.EQ.MP) GOTO 640
3163             ZSAFE=CONT(MP)
3164             CONT(MP)=CONT(I)
3165             CONT(I)=ZSAFE
3166 640         CONTINUE
3167          END DO
3168C -----------------------------------
3169  88      CONTINUE
3170          ERR=0.D0
3171          DO I=1,N
3172             ERR=ERR+(CONT(I)/SCAL(I))**2
3173          END DO
3174          ERR=MAX(SQRT(ERR/N),1.D-10)
3175       END IF
3176       RETURN
3177C
3178C -----------------------------------------------------------
3179C
3180  55  CONTINUE
3181      RETURN
3182       END
3183C
3184C     END OF SUBROUTINE ESTRAV
3185C
3186C ***********************************************************
3187C
3188      SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
3189     &          M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1)
3190      IMPLICIT KPP_REAL (A-H,O-Z)
3191      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),
3192     &          IP(NM1),DY(N),AK(N),FX(N),YNEW(N)
3193      LOGICAL STAGE1
3194      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
3195C
3196      IF (HD.EQ.0.D0) THEN
3197         DO  I=1,N
3198           AK(I)=DY(I)
3199         END DO
3200      ELSE
3201         DO I=1,N
3202            AK(I)=DY(I)+HD*FX(I)
3203         END DO
3204      END IF
3205C
3206      GOTO (1,2,3,4,5,6,55,55,55,55,11,12,13,13,15), IJOB
3207C
3208C -----------------------------------------------------------
3209C
3210   1  CONTINUE
3211C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
3212      IF (STAGE1) THEN
3213         DO I=1,N
3214            AK(I)=AK(I)+YNEW(I)
3215         END DO
3216      END IF
3217      CALL SOL (N,LDE,E,AK,IP)
3218      RETURN
3219C
3220C -----------------------------------------------------------
3221C
3222  11  CONTINUE
3223C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
3224      IF (STAGE1) THEN
3225         DO I=1,N
3226            AK(I)=AK(I)+YNEW(I)
3227         END DO
3228      END IF
3229 48   MM=M1/M2
3230      DO J=1,M2
3231         SUM=0.D0
3232         DO K=MM-1,0,-1
3233            JKM=J+K*M2
3234            SUM=(AK(JKM)+SUM)/FAC1
3235            DO I=1,NM1
3236               IM1=I+M1
3237               AK(IM1)=AK(IM1)+FJAC(I,JKM)*SUM
3238            END DO
3239         END DO
3240      END DO
3241      CALL SOL (NM1,LDE,E,AK(M1+1),IP)
3242      DO I=M1,1,-1
3243         AK(I)=(AK(I)+AK(M2+I))/FAC1
3244      END DO
3245      RETURN
3246C
3247C -----------------------------------------------------------
3248C
3249   2  CONTINUE
3250C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
3251      IF (STAGE1) THEN
3252         DO I=1,N
3253            AK(I)=AK(I)+YNEW(I)
3254         END DO
3255      END IF
3256      CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
3257      RETURN
3258C
3259C -----------------------------------------------------------
3260C
3261  12  CONTINUE
3262C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
3263      IF (STAGE1) THEN
3264         DO I=1,N
3265            AK(I)=AK(I)+YNEW(I)
3266         END DO
3267      END IF
3268  45  MM=M1/M2
3269      DO J=1,M2
3270         SUM=0.D0
3271         DO K=MM-1,0,-1
3272            JKM=J+K*M2
3273            SUM=(AK(JKM)+SUM)/FAC1
3274            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
3275               IM1=I+M1
3276               AK(IM1)=AK(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM
3277            END DO
3278         END DO
3279      END DO
3280      CALL SOLB (NM1,LDE,E,MLE,MUE,AK(M1+1),IP)
3281      DO I=M1,1,-1
3282         AK(I)=(AK(I)+AK(M2+I))/FAC1
3283      END DO
3284      RETURN
3285C
3286C -----------------------------------------------------------
3287C
3288   3  CONTINUE
3289C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
3290      IF (STAGE1) THEN
3291      DO  I=1,N
3292         SUM=0.D0
3293         DO  J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
3294            SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
3295         END DO
3296         AK(I)=AK(I)+SUM
3297      END DO
3298      END IF
3299      CALL SOL (N,LDE,E,AK,IP)
3300      RETURN
3301C
3302C -----------------------------------------------------------
3303C
3304  13  CONTINUE
3305C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
3306      IF (STAGE1) THEN
3307         DO I=1,M1
3308            AK(I)=AK(I)+YNEW(I)
3309         END DO
3310         DO I=1,NM1
3311            SUM=0.D0
3312            DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
3313                SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J+M1)
3314            END DO
3315            IM1=I+M1
3316            AK(IM1)=AK(IM1)+SUM
3317         END DO
3318      END IF
3319      IF (IJOB.EQ.14) GOTO 45
3320      GOTO 48
3321C
3322C -----------------------------------------------------------
3323C
3324   4  CONTINUE
3325C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
3326      IF (STAGE1) THEN
3327      DO I=1,N
3328         SUM=0.D0
3329         DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
3330            SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
3331         END DO
3332         AK(I)=AK(I)+SUM
3333      END DO
3334      END IF
3335      CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
3336      RETURN
3337C
3338C -----------------------------------------------------------
3339C
3340   5  CONTINUE
3341C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
3342      IF (STAGE1) THEN
3343      DO I=1,N
3344         SUM=0.D0
3345         DO J=1,N
3346            SUM=SUM+FMAS(I,J)*YNEW(J)
3347         END DO
3348         AK(I)=AK(I)+SUM
3349      END DO
3350      END IF
3351      CALL SOL (N,LDE,E,AK,IP)
3352      RETURN
3353C
3354C -----------------------------------------------------------
3355C
3356  15  CONTINUE
3357C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
3358      IF (STAGE1) THEN
3359         DO I=1,M1
3360            AK(I)=AK(I)+YNEW(I)
3361         END DO
3362         DO I=1,NM1
3363            SUM=0.D0
3364            DO J=1,NM1
3365               SUM=SUM+FMAS(I,J)*YNEW(J+M1)
3366            END DO
3367            IM1=I+M1
3368            AK(IM1)=AK(IM1)+SUM
3369         END DO
3370      END IF
3371      GOTO 48
3372C
3373C -----------------------------------------------------------
3374C
3375   6  CONTINUE
3376C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
3377C ---  THIS OPTION IS NOT PROVIDED
3378      IF (STAGE1) THEN
3379      DO 624 I=1,N
3380         SUM=0.D0
3381         DO 623 J=1,N
3382  623       SUM=SUM+FMAS(I,J)*YNEW(J)
3383  624    AK(I)=AK(I)+SUM
3384      CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
3385      END IF
3386      RETURN
3387C
3388C -----------------------------------------------------------
3389C
3390  55  CONTINUE
3391      RETURN
3392      END
3393C
3394C     END OF SUBROUTINE SLVROD
3395C
3396C
3397C ***********************************************************
3398C
3399      SUBROUTINE SLVSEU(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
3400     &          M1,M2,NM1,FAC1,E,LDE,IP,IPHES,DEL,IJOB)
3401      IMPLICIT KPP_REAL (A-H,O-Z)
3402      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),DEL(N)
3403      DIMENSION IP(NM1),IPHES(N)
3404      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
3405C
3406      GOTO (1,2,1,2,1,55,7,55,55,55,11,12,11,12,11), IJOB
3407C
3408C -----------------------------------------------------------
3409C
3410   1  CONTINUE
3411C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
3412      CALL SOL (N,LDE,E,DEL,IP)
3413      RETURN
3414C
3415C -----------------------------------------------------------
3416C
3417  11  CONTINUE
3418C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
3419      MM=M1/M2
3420      DO J=1,M2
3421         SUM=0.D0
3422         DO K=MM-1,0,-1
3423            JKM=J+K*M2
3424            SUM=(DEL(JKM)+SUM)/FAC1
3425            DO I=1,NM1
3426               IM1=I+M1
3427               DEL(IM1)=DEL(IM1)+FJAC(I,JKM)*SUM
3428            END DO
3429         END DO
3430      END DO
3431      CALL SOL (NM1,LDE,E,DEL(M1+1),IP)
3432      DO I=M1,1,-1
3433         DEL(I)=(DEL(I)+DEL(M2+I))/FAC1
3434      END DO
3435      RETURN
3436C
3437C -----------------------------------------------------------
3438C
3439   2  CONTINUE
3440C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
3441      CALL SOLB (N,LDE,E,MLE,MUE,DEL,IP)
3442      RETURN
3443C
3444C -----------------------------------------------------------
3445C
3446  12  CONTINUE
3447C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
3448      MM=M1/M2
3449      DO J=1,M2
3450         SUM=0.D0
3451         DO K=MM-1,0,-1
3452            JKM=J+K*M2
3453            SUM=(DEL(JKM)+SUM)/FAC1
3454            DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
3455               IM1=I+M1
3456               DEL(IM1)=DEL(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM
3457            END DO
3458         END DO
3459      END DO
3460      CALL SOLB (NM1,LDE,E,MLE,MUE,DEL(M1+1),IP)
3461      DO I=M1,1,-1
3462         DEL(I)=(DEL(I)+DEL(M2+I))/FAC1
3463      END DO
3464      RETURN
3465C
3466C -----------------------------------------------------------
3467C
3468   7  CONTINUE
3469C ---  HESSENBERG OPTION
3470      DO MMM=N-2,1,-1
3471         MP=N-MMM
3472         MP1=MP-1
3473         I=IPHES(MP)
3474         IF (I.EQ.MP) GOTO 110
3475         ZSAFE=DEL(MP)
3476         DEL(MP)=DEL(I)
3477         DEL(I)=ZSAFE
3478 110     CONTINUE
3479         DO I=MP+1,N
3480            DEL(I)=DEL(I)-FJAC(I,MP1)*DEL(MP)
3481         END DO
3482      END DO
3483      CALL SOLH(N,LDE,E,1,DEL,IP)
3484      DO MMM=1,N-2
3485         MP=N-MMM
3486         MP1=MP-1
3487         DO I=MP+1,N
3488            DEL(I)=DEL(I)+FJAC(I,MP1)*DEL(MP)
3489         END DO
3490         I=IPHES(MP)
3491         IF (I.EQ.MP) GOTO 240
3492         ZSAFE=DEL(MP)
3493         DEL(MP)=DEL(I)
3494         DEL(I)=ZSAFE
3495 240     CONTINUE
3496      END DO
3497      RETURN
3498C
3499C -----------------------------------------------------------
3500C
3501  55  CONTINUE
3502      RETURN
3503      END
3504C
3505C     END OF SUBROUTINE SLVSEU
3506C
3507      SUBROUTINE DEC (N, NDIM, A, IP, IER)
3508C VERSION REAL KPP_REAL
3509      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
3510      KPP_REAL A,T
3511      DIMENSION A(NDIM,N), IP(N)
3512C-----------------------------------------------------------------------
3513C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
3514C  INPUT..
3515C     N = ORDER OF MATRIX.
3516C     NDIM = DECLARED DIMENSION OF ARRAY  A .
3517C     A = MATRIX TO BE TRIANGULARIZED.
3518C  OUTPUT..
3519C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
3520C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3521C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3522C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3523C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3524C           SINGULAR AT STAGE K.
3525C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3526C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
3527C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3528C
3529C  REFERENCE..
3530C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3531C     C.A.C.M. 15 (1972), P. 274.
3532C-----------------------------------------------------------------------
3533      IER = 0
3534      IP(N) = 1
3535      IF (N .EQ. 1) GO TO 70
3536      NM1 = N - 1
3537      DO 60 K = 1,NM1
3538        KP1 = K + 1
3539        M = K
3540        DO 10 I = KP1,N
3541          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
3542 10     CONTINUE
3543        IP(K) = M
3544        T = A(M,K)
3545        IF (M .EQ. K) GO TO 20
3546        IP(N) = -IP(N)
3547        A(M,K) = A(K,K)
3548        A(K,K) = T
3549 20     CONTINUE
3550        IF (T .EQ. 0.D0) GO TO 80
3551        T = 1.D0/T
3552        DO 30 I = KP1,N
3553 30       A(I,K) = -A(I,K)*T
3554        DO 50 J = KP1,N
3555          T = A(M,J)
3556          A(M,J) = A(K,J)
3557          A(K,J) = T
3558          IF (T .EQ. 0.D0) GO TO 45
3559          DO 40 I = KP1,N
3560 40         A(I,J) = A(I,J) + A(I,K)*T
3561 45       CONTINUE
3562 50       CONTINUE
3563 60     CONTINUE
3564 70   K = N
3565      IF (A(N,N) .EQ. 0.D0) GO TO 80
3566      RETURN
3567 80   IER = K
3568      IP(N) = 0
3569      RETURN
3570C----------------------- END OF SUBROUTINE DEC -------------------------
3571      END
3572C
3573C
3574      SUBROUTINE SOL (N, NDIM, A, B, IP)
3575C VERSION REAL KPP_REAL
3576      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
3577      KPP_REAL A,B,T
3578      DIMENSION A(NDIM,N), B(N), IP(N)
3579C-----------------------------------------------------------------------
3580C  SOLUTION OF LINEAR SYSTEM, A*X = B .
3581C  INPUT..
3582C    N = ORDER OF MATRIX.
3583C    NDIM = DECLARED DIMENSION OF ARRAY  A .
3584C    A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
3585C    B = RIGHT HAND SIDE VECTOR.
3586C    IP = PIVOT VECTOR OBTAINED FROM DEC.
3587C  DO NOT USE IF DEC HAS SET IER .NE. 0.
3588C  OUTPUT..
3589C    B = SOLUTION VECTOR, X .
3590C-----------------------------------------------------------------------
3591      IF (N .EQ. 1) GO TO 50
3592      NM1 = N - 1
3593      DO 20 K = 1,NM1
3594        KP1 = K + 1
3595        M = IP(K)
3596        T = B(M)
3597        B(M) = B(K)
3598        B(K) = T
3599        DO 10 I = KP1,N
3600 10       B(I) = B(I) + A(I,K)*T
3601 20     CONTINUE
3602      DO 40 KB = 1,NM1
3603        KM1 = N - KB
3604        K = KM1 + 1
3605        B(K) = B(K)/A(K,K)
3606        T = -B(K)
3607        DO 30 I = 1,KM1
3608 30       B(I) = B(I) + A(I,K)*T
3609 40     CONTINUE
3610 50   B(1) = B(1)/A(1,1)
3611      RETURN
3612C----------------------- END OF SUBROUTINE SOL -------------------------
3613      END
3614c
3615c
3616      SUBROUTINE DECH (N, NDIM, A, LB, IP, IER)
3617C VERSION REAL KPP_REAL
3618      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J,LB,NA
3619      KPP_REAL A,T
3620      DIMENSION A(NDIM,N), IP(N)
3621C-----------------------------------------------------------------------
3622C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A HESSENBERG
3623C  MATRIX WITH LOWER BANDWIDTH LB
3624C  INPUT..
3625C     N = ORDER OF MATRIX A.
3626C     NDIM = DECLARED DIMENSION OF ARRAY  A .
3627C     A = MATRIX TO BE TRIANGULARIZED.
3628C     LB = LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED, LB.GE.1).
3629C  OUTPUT..
3630C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
3631C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3632C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3633C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3634C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3635C           SINGULAR AT STAGE K.
3636C  USE  SOLH  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3637C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
3638C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3639C
3640C  REFERENCE..
3641C     THIS IS A SLIGHT MODIFICATION OF
3642C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3643C     C.A.C.M. 15 (1972), P. 274.
3644C-----------------------------------------------------------------------
3645      IER = 0
3646      IP(N) = 1
3647      IF (N .EQ. 1) GO TO 70
3648      NM1 = N - 1
3649      DO 60 K = 1,NM1
3650        KP1 = K + 1
3651        M = K
3652        NA = MIN0(N,LB+K)
3653        DO 10 I = KP1,NA
3654          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
3655 10     CONTINUE
3656        IP(K) = M
3657        T = A(M,K)
3658        IF (M .EQ. K) GO TO 20
3659        IP(N) = -IP(N)
3660        A(M,K) = A(K,K)
3661        A(K,K) = T
3662 20     CONTINUE
3663        IF (T .EQ. 0.D0) GO TO 80
3664        T = 1.D0/T
3665        DO 30 I = KP1,NA
3666 30       A(I,K) = -A(I,K)*T
3667        DO 50 J = KP1,N
3668          T = A(M,J)
3669          A(M,J) = A(K,J)
3670          A(K,J) = T
3671          IF (T .EQ. 0.D0) GO TO 45
3672          DO 40 I = KP1,NA
3673 40         A(I,J) = A(I,J) + A(I,K)*T
3674 45       CONTINUE
3675 50       CONTINUE
3676 60     CONTINUE
3677 70   K = N
3678      IF (A(N,N) .EQ. 0.D0) GO TO 80
3679      RETURN
3680 80   IER = K
3681      IP(N) = 0
3682      RETURN
3683C----------------------- END OF SUBROUTINE DECH ------------------------
3684      END
3685C
3686C
3687      SUBROUTINE SOLH (N, NDIM, A, LB, B, IP)
3688C VERSION REAL KPP_REAL
3689      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1,LB,NA
3690      KPP_REAL A,B,T
3691      DIMENSION A(NDIM,N), B(N), IP(N)
3692C-----------------------------------------------------------------------
3693C  SOLUTION OF LINEAR SYSTEM, A*X = B .
3694C  INPUT..
3695C    N = ORDER OF MATRIX A.
3696C    NDIM = DECLARED DIMENSION OF ARRAY  A .
3697C    A = TRIANGULARIZED MATRIX OBTAINED FROM DECH.
3698C    LB = LOWER BANDWIDTH OF A.
3699C    B = RIGHT HAND SIDE VECTOR.
3700C    IP = PIVOT VECTOR OBTAINED FROM DEC.
3701C  DO NOT USE IF DECH HAS SET IER .NE. 0.
3702C  OUTPUT..
3703C    B = SOLUTION VECTOR, X .
3704C-----------------------------------------------------------------------
3705      IF (N .EQ. 1) GO TO 50
3706      NM1 = N - 1
3707      DO 20 K = 1,NM1
3708        KP1 = K + 1
3709        M = IP(K)
3710        T = B(M)
3711        B(M) = B(K)
3712        B(K) = T
3713        NA = MIN0(N,LB+K)
3714        DO 10 I = KP1,NA
3715 10       B(I) = B(I) + A(I,K)*T
3716 20     CONTINUE
3717      DO 40 KB = 1,NM1
3718        KM1 = N - KB
3719        K = KM1 + 1
3720        B(K) = B(K)/A(K,K)
3721        T = -B(K)
3722        DO 30 I = 1,KM1
3723 30       B(I) = B(I) + A(I,K)*T
3724 40     CONTINUE
3725 50   B(1) = B(1)/A(1,1)
3726      RETURN
3727C----------------------- END OF SUBROUTINE SOLH ------------------------
3728      END
3729C
3730      SUBROUTINE DECC (N, NDIM, AR, AI, IP, IER)
3731C VERSION COMPLEX KPP_REAL
3732      IMPLICIT KPP_REAL (A-H,O-Z)
3733      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
3734      DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N)
3735C-----------------------------------------------------------------------
3736C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION
3737C  ------ MODIFICATION FOR COMPLEX MATRICES --------
3738C  INPUT..
3739C     N = ORDER OF MATRIX.
3740C     NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI .
3741C     (AR, AI) = MATRIX TO BE TRIANGULARIZED.
3742C  OUTPUT..
3743C     AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART.
3744C     AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART.
3745C     AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3746C                                                    REAL PART.
3747C     AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3748C                                                    IMAGINARY PART.
3749C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3750C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3751C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3752C           SINGULAR AT STAGE K.
3753C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3754C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3755C
3756C  REFERENCE..
3757C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3758C     C.A.C.M. 15 (1972), P. 274.
3759C-----------------------------------------------------------------------
3760      IER = 0
3761      IP(N) = 1
3762      IF (N .EQ. 1) GO TO 70
3763      NM1 = N - 1
3764      DO 60 K = 1,NM1
3765        KP1 = K + 1
3766        M = K
3767        DO 10 I = KP1,N
3768          IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT.
3769     &          DABS(AR(M,K))+DABS(AI(M,K))) M = I
3770 10     CONTINUE
3771        IP(K) = M
3772        TR = AR(M,K)
3773        TI = AI(M,K)
3774        IF (M .EQ. K) GO TO 20
3775        IP(N) = -IP(N)
3776        AR(M,K) = AR(K,K)
3777        AI(M,K) = AI(K,K)
3778        AR(K,K) = TR
3779        AI(K,K) = TI
3780 20     CONTINUE
3781        IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80
3782        DEN=TR*TR+TI*TI
3783        TR=TR/DEN
3784        TI=-TI/DEN
3785        DO 30 I = KP1,N
3786          PRODR=AR(I,K)*TR-AI(I,K)*TI
3787          PRODI=AI(I,K)*TR+AR(I,K)*TI
3788          AR(I,K)=-PRODR
3789          AI(I,K)=-PRODI
3790 30       CONTINUE
3791        DO 50 J = KP1,N
3792          TR = AR(M,J)
3793          TI = AI(M,J)
3794          AR(M,J) = AR(K,J)
3795          AI(M,J) = AI(K,J)
3796          AR(K,J) = TR
3797          AI(K,J) = TI
3798          IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48
3799          IF (TI .EQ. 0.D0) THEN
3800            DO 40 I = KP1,N
3801            PRODR=AR(I,K)*TR
3802            PRODI=AI(I,K)*TR
3803            AR(I,J) = AR(I,J) + PRODR
3804            AI(I,J) = AI(I,J) + PRODI
3805 40         CONTINUE
3806            GO TO 48
3807          END IF
3808          IF (TR .EQ. 0.D0) THEN
3809            DO 45 I = KP1,N
3810            PRODR=-AI(I,K)*TI
3811            PRODI=AR(I,K)*TI
3812            AR(I,J) = AR(I,J) + PRODR
3813            AI(I,J) = AI(I,J) + PRODI
3814 45         CONTINUE
3815            GO TO 48
3816          END IF
3817          DO 47 I = KP1,N
3818            PRODR=AR(I,K)*TR-AI(I,K)*TI
3819            PRODI=AI(I,K)*TR+AR(I,K)*TI
3820            AR(I,J) = AR(I,J) + PRODR
3821            AI(I,J) = AI(I,J) + PRODI
3822 47         CONTINUE
3823 48       CONTINUE
3824 50       CONTINUE
3825 60     CONTINUE
3826 70   K = N
3827      IF (DABS(AR(N,N))+DABS(AI(N,N)) .EQ. 0.D0) GO TO 80
3828      RETURN
3829 80   IER = K
3830      IP(N) = 0
3831      RETURN
3832C----------------------- END OF SUBROUTINE DECC ------------------------
3833      END
3834C
3835C
3836      SUBROUTINE SOLC (N, NDIM, AR, AI, BR, BI, IP)
3837C VERSION COMPLEX KPP_REAL
3838      IMPLICIT KPP_REAL (A-H,O-Z)
3839      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
3840      DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N)
3841C-----------------------------------------------------------------------
3842C  SOLUTION OF LINEAR SYSTEM, A*X = B .
3843C  INPUT..
3844C    N = ORDER OF MATRIX.
3845C    NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI.
3846C    (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
3847C    (BR,BI) = RIGHT HAND SIDE VECTOR.
3848C    IP = PIVOT VECTOR OBTAINED FROM DEC.
3849C  DO NOT USE IF DEC HAS SET IER .NE. 0.
3850C  OUTPUT..
3851C    (BR,BI) = SOLUTION VECTOR, X .
3852C-----------------------------------------------------------------------
3853      IF (N .EQ. 1) GO TO 50
3854      NM1 = N - 1
3855      DO 20 K = 1,NM1
3856        KP1 = K + 1
3857        M = IP(K)
3858        TR = BR(M)
3859        TI = BI(M)
3860        BR(M) = BR(K)
3861        BI(M) = BI(K)
3862        BR(K) = TR
3863        BI(K) = TI
3864        DO 10 I = KP1,N
3865          PRODR=AR(I,K)*TR-AI(I,K)*TI
3866          PRODI=AI(I,K)*TR+AR(I,K)*TI
3867          BR(I) = BR(I) + PRODR
3868          BI(I) = BI(I) + PRODI
3869 10       CONTINUE
3870 20     CONTINUE
3871      DO 40 KB = 1,NM1
3872        KM1 = N - KB
3873        K = KM1 + 1
3874        DEN=AR(K,K)*AR(K,K)+AI(K,K)*AI(K,K)
3875        PRODR=BR(K)*AR(K,K)+BI(K)*AI(K,K)
3876        PRODI=BI(K)*AR(K,K)-BR(K)*AI(K,K)
3877        BR(K)=PRODR/DEN
3878        BI(K)=PRODI/DEN
3879        TR = -BR(K)
3880        TI = -BI(K)
3881        DO 30 I = 1,KM1
3882          PRODR=AR(I,K)*TR-AI(I,K)*TI
3883          PRODI=AI(I,K)*TR+AR(I,K)*TI
3884          BR(I) = BR(I) + PRODR
3885          BI(I) = BI(I) + PRODI
3886 30       CONTINUE
3887 40     CONTINUE
3888 50     CONTINUE
3889        DEN=AR(1,1)*AR(1,1)+AI(1,1)*AI(1,1)
3890        PRODR=BR(1)*AR(1,1)+BI(1)*AI(1,1)
3891        PRODI=BI(1)*AR(1,1)-BR(1)*AI(1,1)
3892        BR(1)=PRODR/DEN
3893        BI(1)=PRODI/DEN
3894      RETURN
3895C----------------------- END OF SUBROUTINE SOLC ------------------------
3896      END
3897C
3898C
3899      SUBROUTINE DECHC (N, NDIM, AR, AI, LB, IP, IER)
3900C VERSION COMPLEX KPP_REAL
3901      IMPLICIT KPP_REAL (A-H,O-Z)
3902      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
3903      DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N)
3904C-----------------------------------------------------------------------
3905C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION
3906C  ------ MODIFICATION FOR COMPLEX MATRICES --------
3907C  INPUT..
3908C     N = ORDER OF MATRIX.
3909C     NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI .
3910C     (AR, AI) = MATRIX TO BE TRIANGULARIZED.
3911C  OUTPUT..
3912C     AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART.
3913C     AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART.
3914C     AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3915C                                                    REAL PART.
3916C     AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3917C                                                    IMAGINARY PART.
3918C     LB = LOWER BANDWIDTH OF A (DIAGONAL NOT COUNTED), LB.GE.1.
3919C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3920C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3921C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3922C           SINGULAR AT STAGE K.
3923C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3924C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3925C
3926C  REFERENCE..
3927C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3928C     C.A.C.M. 15 (1972), P. 274.
3929C-----------------------------------------------------------------------
3930      IER = 0
3931      IP(N) = 1
3932      IF (LB .EQ. 0) GO TO 70
3933      IF (N .EQ. 1) GO TO 70
3934      NM1 = N - 1
3935      DO 60 K = 1,NM1
3936        KP1 = K + 1
3937        M = K
3938        NA = MIN0(N,LB+K)
3939        DO 10 I = KP1,NA
3940          IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT.
3941     &          DABS(AR(M,K))+DABS(AI(M,K))) M = I
3942 10     CONTINUE
3943        IP(K) = M
3944        TR = AR(M,K)
3945        TI = AI(M,K)
3946        IF (M .EQ. K) GO TO 20
3947        IP(N) = -IP(N)
3948        AR(M,K) = AR(K,K)
3949        AI(M,K) = AI(K,K)
3950        AR(K,K) = TR
3951        AI(K,K) = TI
3952 20     CONTINUE
3953        IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80
3954        DEN=TR*TR+TI*TI
3955        TR=TR/DEN
3956        TI=-TI/DEN
3957        DO 30 I = KP1,NA
3958          PRODR=AR(I,K)*TR-AI(I,K)*TI
3959          PRODI=AI(I,K)*TR+AR(I,K)*TI
3960          AR(I,K)=-PRODR
3961          AI(I,K)=-PRODI
3962 30       CONTINUE
3963        DO 50 J = KP1,N
3964          TR = AR(M,J)
3965          TI = AI(M,J)
3966          AR(M,J) = AR(K,J)
3967          AI(M,J) = AI(K,J)
3968          AR(K,J) = TR
3969          AI(K,J) = TI
3970          IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48
3971          IF (TI .EQ. 0.D0) THEN
3972            DO 40 I = KP1,NA
3973            PRODR=AR(I,K)*TR
3974            PRODI=AI(I,K)*TR
3975            AR(I,J) = AR(I,J) + PRODR
3976            AI(I,J) = AI(I,J) + PRODI
3977 40         CONTINUE
3978            GO TO 48
3979          END IF
3980          IF (TR .EQ. 0.D0) THEN
3981            DO 45 I = KP1,NA
3982            PRODR=-AI(I,K)*TI
3983            PRODI=AR(I,K)*TI
3984            AR(I,J) = AR(I,J) + PRODR
3985            AI(I,J) = AI(I,J) + PRODI
3986 45         CONTINUE
3987            GO TO 48
3988          END IF
3989          DO 47 I = KP1,NA
3990            PRODR=AR(I,K)*TR-AI(I,K)*TI
3991            PRODI=AI(I,K)*TR+AR(I,K)*TI
3992            AR(I,J) = AR(I,J) + PRODR
3993            AI(I,J) = AI(I,J) + PRODI
3994 47         CONTINUE
3995 48       CONTINUE
3996 50       CONTINUE
3997 60     CONTINUE
3998 70   K = N
3999      IF (DABS(AR(N,N))+DABS(AI(N,N)) .EQ. 0.D0) GO TO 80
4000      RETURN
4001 80   IER = K
4002      IP(N) = 0
4003      RETURN
4004C----------------------- END OF SUBROUTINE DECHC -----------------------
4005      END
4006C
4007C
4008      SUBROUTINE SOLHC (N, NDIM, AR, AI, LB, BR, BI, IP)
4009C VERSION COMPLEX KPP_REAL
4010      IMPLICIT KPP_REAL (A-H,O-Z)
4011      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
4012      DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N)
4013C-----------------------------------------------------------------------
4014C  SOLUTION OF LINEAR SYSTEM, A*X = B .
4015C  INPUT..
4016C    N = ORDER OF MATRIX.
4017C    NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI.
4018C    (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
4019C    (BR,BI) = RIGHT HAND SIDE VECTOR.
4020C    LB = LOWER BANDWIDTH OF A.
4021C    IP = PIVOT VECTOR OBTAINED FROM DEC.
4022C  DO NOT USE IF DEC HAS SET IER .NE. 0.
4023C  OUTPUT..
4024C    (BR,BI) = SOLUTION VECTOR, X .
4025C-----------------------------------------------------------------------
4026      IF (N .EQ. 1) GO TO 50
4027      NM1 = N - 1
4028      IF (LB .EQ. 0) GO TO 25
4029      DO 20 K = 1,NM1
4030        KP1 = K + 1
4031        M = IP(K)
4032        TR = BR(M)
4033        TI = BI(M)
4034        BR(M) = BR(K)
4035        BI(M) = BI(K)
4036        BR(K) = TR
4037        BI(K) = TI
4038        DO 10 I = KP1,MIN0(N,LB+K)
4039          PRODR=AR(I,K)*TR-AI(I,K)*TI
4040          PRODI=AI(I,K)*TR+AR(I,K)*TI
4041          BR(I) = BR(I) + PRODR
4042          BI(I) = BI(I) + PRODI
4043 10       CONTINUE
4044 20     CONTINUE
4045 25     CONTINUE
4046      DO 40 KB = 1,NM1
4047        KM1 = N - KB
4048        K = KM1 + 1
4049        DEN=AR(K,K)*AR(K,K)+AI(K,K)*AI(K,K)
4050        PRODR=BR(K)*AR(K,K)+BI(K)*AI(K,K)
4051        PRODI=BI(K)*AR(K,K)-BR(K)*AI(K,K)
4052        BR(K)=PRODR/DEN
4053        BI(K)=PRODI/DEN
4054        TR = -BR(K)
4055        TI = -BI(K)
4056        DO 30 I = 1,KM1
4057          PRODR=AR(I,K)*TR-AI(I,K)*TI
4058          PRODI=AI(I,K)*TR+AR(I,K)*TI
4059          BR(I) = BR(I) + PRODR
4060          BI(I) = BI(I) + PRODI
4061 30       CONTINUE
4062 40     CONTINUE
4063 50     CONTINUE
4064        DEN=AR(1,1)*AR(1,1)+AI(1,1)*AI(1,1)
4065        PRODR=BR(1)*AR(1,1)+BI(1)*AI(1,1)
4066        PRODI=BI(1)*AR(1,1)-BR(1)*AI(1,1)
4067        BR(1)=PRODR/DEN
4068        BI(1)=PRODI/DEN
4069      RETURN
4070C----------------------- END OF SUBROUTINE SOLHC -----------------------
4071      END
4072C
4073      SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER)
4074      KPP_REAL A,T
4075      DIMENSION A(NDIM,N), IP(N)
4076C-----------------------------------------------------------------------
4077C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
4078C  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
4079C  INPUT..
4080C     N       ORDER OF THE ORIGINAL MATRIX A.
4081C     NDIM    DECLARED DIMENSION OF ARRAY  A.
4082C     A       CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
4083C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  A  AND
4084C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
4085C                ML+1 THROUGH 2*ML+MU+1 OF  A.
4086C     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4087C     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4088C  OUTPUT..
4089C     A       AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
4090C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
4091C     IP      INDEX VECTOR OF PIVOT INDICES.
4092C     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
4093C     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
4094C                SINGULAR AT STAGE K.
4095C  USE  SOLB  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
4096C  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
4097C  IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
4098C
4099C  REFERENCE..
4100C     THIS IS A MODIFICATION OF
4101C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
4102C     C.A.C.M. 15 (1972), P. 274.
4103C-----------------------------------------------------------------------
4104      IER = 0
4105      IP(N) = 1
4106      MD = ML + MU + 1
4107      MD1 = MD + 1
4108      JU = 0
4109      IF (ML .EQ. 0) GO TO 70
4110      IF (N .EQ. 1) GO TO 70
4111      IF (N .LT. MU+2) GO TO 7
4112      DO 5 J = MU+2,N
4113      DO 5 I = 1,ML
4114  5   A(I,J) = 0.D0
4115  7   NM1 = N - 1
4116      DO 60 K = 1,NM1
4117        KP1 = K + 1
4118        M = MD
4119        MDL = MIN(ML,N-K) + MD
4120        DO 10 I = MD1,MDL
4121          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
4122 10     CONTINUE
4123        IP(K) = M + K - MD
4124        T = A(M,K)
4125        IF (M .EQ. MD) GO TO 20
4126        IP(N) = -IP(N)
4127        A(M,K) = A(MD,K)
4128        A(MD,K) = T
4129 20     CONTINUE
4130        IF (T .EQ. 0.D0) GO TO 80
4131        T = 1.D0/T
4132        DO 30 I = MD1,MDL
4133 30       A(I,K) = -A(I,K)*T
4134        JU = MIN0(MAX0(JU,MU+IP(K)),N)
4135        MM = MD
4136        IF (JU .LT. KP1) GO TO 55
4137        DO 50 J = KP1,JU
4138          M = M - 1
4139          MM = MM - 1
4140          T = A(M,J)
4141          IF (M .EQ. MM) GO TO 35
4142          A(M,J) = A(MM,J)
4143          A(MM,J) = T
4144 35       CONTINUE
4145          IF (T .EQ. 0.D0) GO TO 45
4146          JK = J - K
4147          DO 40 I = MD1,MDL
4148            IJK = I - JK
4149 40         A(IJK,J) = A(IJK,J) + A(I,K)*T
4150 45       CONTINUE
4151 50       CONTINUE
4152 55     CONTINUE
4153 60     CONTINUE
4154 70   K = N
4155      IF (A(MD,N) .EQ. 0.D0) GO TO 80
4156      RETURN
4157 80   IER = K
4158      IP(N) = 0
4159      RETURN
4160C----------------------- END OF SUBROUTINE DECB ------------------------
4161      END
4162C
4163C
4164      SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP)
4165      KPP_REAL A,B,T
4166      DIMENSION A(NDIM,N), B(N), IP(N)
4167C-----------------------------------------------------------------------
4168C  SOLUTION OF LINEAR SYSTEM, A*X = B .
4169C  INPUT..
4170C    N      ORDER OF MATRIX A.
4171C    NDIM   DECLARED DIMENSION OF ARRAY  A .
4172C    A      TRIANGULARIZED MATRIX OBTAINED FROM DECB.
4173C    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4174C    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4175C    B      RIGHT HAND SIDE VECTOR.
4176C    IP     PIVOT VECTOR OBTAINED FROM DECB.
4177C  DO NOT USE IF DECB HAS SET IER .NE. 0.
4178C  OUTPUT..
4179C    B      SOLUTION VECTOR, X .
4180C-----------------------------------------------------------------------
4181      MD = ML + MU + 1
4182      MD1 = MD + 1
4183      MDM = MD - 1
4184      NM1 = N - 1
4185      IF (ML .EQ. 0) GO TO 25
4186      IF (N .EQ. 1) GO TO 50
4187      DO 20 K = 1,NM1
4188        M = IP(K)
4189        T = B(M)
4190        B(M) = B(K)
4191        B(K) = T
4192        MDL = MIN(ML,N-K) + MD
4193        DO 10 I = MD1,MDL
4194          IMD = I + K - MD
4195 10       B(IMD) = B(IMD) + A(I,K)*T
4196 20     CONTINUE
4197 25   CONTINUE
4198      DO 40 KB = 1,NM1
4199        K = N + 1 - KB
4200        B(K) = B(K)/A(MD,K)
4201        T = -B(K)
4202        KMD = MD - K
4203        LM = MAX0(1,KMD+1)
4204        DO 30 I = LM,MDM
4205          IMD = I - KMD
4206 30       B(IMD) = B(IMD) + A(I,K)*T
4207 40     CONTINUE
4208 50   B(1) = B(1)/A(MD,1)
4209      RETURN
4210C----------------------- END OF SUBROUTINE SOLB ------------------------
4211      END
4212C
4213      SUBROUTINE DECBC (N, NDIM, AR, AI, ML, MU, IP, IER)
4214      IMPLICIT KPP_REAL (A-H,O-Z)
4215      DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N)
4216C-----------------------------------------------------------------------
4217C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED COMPLEX
4218C  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
4219C  INPUT..
4220C     N       ORDER OF THE ORIGINAL MATRIX A.
4221C     NDIM    DECLARED DIMENSION OF ARRAY  A.
4222C     AR, AI     CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
4223C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  AR (REAL
4224C                PART) AND AI (IMAGINARY PART)  AND
4225C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
4226C                ML+1 THROUGH 2*ML+MU+1 OF  AR AND AI.
4227C     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4228C     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4229C  OUTPUT..
4230C     AR, AI  AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
4231C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
4232C     IP      INDEX VECTOR OF PIVOT INDICES.
4233C     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
4234C     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
4235C                SINGULAR AT STAGE K.
4236C  USE  SOLBC  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
4237C  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
4238C  IF IP(N)=O, A IS SINGULAR, SOLBC WILL DIVIDE BY ZERO.
4239C
4240C  REFERENCE..
4241C     THIS IS A MODIFICATION OF
4242C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
4243C     C.A.C.M. 15 (1972), P. 274.
4244C-----------------------------------------------------------------------
4245      IER = 0
4246      IP(N) = 1
4247      MD = ML + MU + 1
4248      MD1 = MD + 1
4249      JU = 0
4250      IF (ML .EQ. 0) GO TO 70
4251      IF (N .EQ. 1) GO TO 70
4252      IF (N .LT. MU+2) GO TO 7
4253      DO 5 J = MU+2,N
4254      DO 5 I = 1,ML
4255      AR(I,J) = 0.D0
4256      AI(I,J) = 0.D0
4257  5   CONTINUE
4258  7   NM1 = N - 1
4259      DO 60 K = 1,NM1
4260        KP1 = K + 1
4261        M = MD
4262        MDL = MIN(ML,N-K) + MD
4263        DO 10 I = MD1,MDL
4264          IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT.
4265     &          DABS(AR(M,K))+DABS(AI(M,K))) M = I
4266 10     CONTINUE
4267        IP(K) = M + K - MD
4268        TR = AR(M,K)
4269        TI = AI(M,K)
4270        IF (M .EQ. MD) GO TO 20
4271        IP(N) = -IP(N)
4272        AR(M,K) = AR(MD,K)
4273        AI(M,K) = AI(MD,K)
4274        AR(MD,K) = TR
4275        AI(MD,K) = TI
4276 20     IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80
4277        DEN=TR*TR+TI*TI
4278        TR=TR/DEN
4279        TI=-TI/DEN
4280        DO 30 I = MD1,MDL
4281          PRODR=AR(I,K)*TR-AI(I,K)*TI
4282          PRODI=AI(I,K)*TR+AR(I,K)*TI
4283          AR(I,K)=-PRODR
4284          AI(I,K)=-PRODI
4285 30       CONTINUE
4286        JU = MIN0(MAX0(JU,MU+IP(K)),N)
4287        MM = MD
4288        IF (JU .LT. KP1) GO TO 55
4289        DO 50 J = KP1,JU
4290          M = M - 1
4291          MM = MM - 1
4292          TR = AR(M,J)
4293          TI = AI(M,J)
4294          IF (M .EQ. MM) GO TO 35
4295          AR(M,J) = AR(MM,J)
4296          AI(M,J) = AI(MM,J)
4297          AR(MM,J) = TR
4298          AI(MM,J) = TI
4299 35       CONTINUE
4300          IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48
4301          JK = J - K
4302          IF (TI .EQ. 0.D0) THEN
4303            DO 40 I = MD1,MDL
4304            IJK = I - JK
4305            PRODR=AR(I,K)*TR
4306            PRODI=AI(I,K)*TR
4307            AR(IJK,J) = AR(IJK,J) + PRODR
4308            AI(IJK,J) = AI(IJK,J) + PRODI
4309 40         CONTINUE
4310            GO TO 48
4311          END IF
4312          IF (TR .EQ. 0.D0) THEN
4313            DO 45 I = MD1,MDL
4314            IJK = I - JK
4315            PRODR=-AI(I,K)*TI
4316            PRODI=AR(I,K)*TI
4317            AR(IJK,J) = AR(IJK,J) + PRODR
4318            AI(IJK,J) = AI(IJK,J) + PRODI
4319 45         CONTINUE
4320            GO TO 48
4321          END IF
4322          DO 47 I = MD1,MDL
4323            IJK = I - JK
4324            PRODR=AR(I,K)*TR-AI(I,K)*TI
4325            PRODI=AI(I,K)*TR+AR(I,K)*TI
4326            AR(IJK,J) = AR(IJK,J) + PRODR
4327            AI(IJK,J) = AI(IJK,J) + PRODI
4328 47         CONTINUE
4329 48       CONTINUE
4330 50       CONTINUE
4331 55     CONTINUE
4332 60     CONTINUE
4333 70   K = N
4334      IF (DABS(AR(MD,N))+DABS(AI(MD,N)) .EQ. 0.D0) GO TO 80
4335      RETURN
4336 80   IER = K
4337      IP(N) = 0
4338      RETURN
4339C----------------------- END OF SUBROUTINE DECBC ------------------------
4340      END
4341C
4342C
4343      SUBROUTINE SOLBC (N, NDIM, AR, AI, ML, MU, BR, BI, IP)
4344      IMPLICIT KPP_REAL (A-H,O-Z)
4345      DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N)
4346C-----------------------------------------------------------------------
4347C  SOLUTION OF LINEAR SYSTEM, A*X = B ,
4348C                  VERSION BANDED AND COMPLEX-KPP_REAL.
4349C  INPUT..
4350C    N      ORDER OF MATRIX A.
4351C    NDIM   DECLARED DIMENSION OF ARRAY  A .
4352C    AR, AI TRIANGULARIZED MATRIX OBTAINED FROM DECB (REAL AND IMAG. PART).
4353C    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4354C    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4355C    BR, BI RIGHT HAND SIDE VECTOR (REAL AND IMAG. PART).
4356C    IP     PIVOT VECTOR OBTAINED FROM DECBC.
4357C  DO NOT USE IF DECB HAS SET IER .NE. 0.
4358C  OUTPUT..
4359C    BR, BI SOLUTION VECTOR, X (REAL AND IMAG. PART).
4360C-----------------------------------------------------------------------
4361      MD = ML + MU + 1
4362      MD1 = MD + 1
4363      MDM = MD - 1
4364      NM1 = N - 1
4365      IF (ML .EQ. 0) GO TO 25
4366      IF (N .EQ. 1) GO TO 50
4367      DO 20 K = 1,NM1
4368        M = IP(K)
4369        TR = BR(M)
4370        TI = BI(M)
4371        BR(M) = BR(K)
4372        BI(M) = BI(K)
4373        BR(K) = TR
4374        BI(K) = TI
4375        MDL = MIN(ML,N-K) + MD
4376        DO 10 I = MD1,MDL
4377          IMD = I + K - MD
4378          PRODR=AR(I,K)*TR-AI(I,K)*TI
4379          PRODI=AI(I,K)*TR+AR(I,K)*TI
4380          BR(IMD) = BR(IMD) + PRODR
4381          BI(IMD) = BI(IMD) + PRODI
4382 10     CONTINUE
4383 20     CONTINUE
4384 25     CONTINUE
4385      DO 40 KB = 1,NM1
4386        K = N + 1 - KB
4387        DEN=AR(MD,K)*AR(MD,K)+AI(MD,K)*AI(MD,K)
4388        PRODR=BR(K)*AR(MD,K)+BI(K)*AI(MD,K)
4389        PRODI=BI(K)*AR(MD,K)-BR(K)*AI(MD,K)
4390        BR(K)=PRODR/DEN
4391        BI(K)=PRODI/DEN
4392        TR = -BR(K)
4393        TI = -BI(K)
4394        KMD = MD - K
4395        LM = MAX0(1,KMD+1)
4396        DO 30 I = LM,MDM
4397          IMD = I - KMD
4398          PRODR=AR(I,K)*TR-AI(I,K)*TI
4399          PRODI=AI(I,K)*TR+AR(I,K)*TI
4400          BR(IMD) = BR(IMD) + PRODR
4401          BI(IMD) = BI(IMD) + PRODI
4402 30       CONTINUE
4403 40     CONTINUE
4404        DEN=AR(MD,1)*AR(MD,1)+AI(MD,1)*AI(MD,1)
4405        PRODR=BR(1)*AR(MD,1)+BI(1)*AI(MD,1)
4406        PRODI=BI(1)*AR(MD,1)-BR(1)*AI(MD,1)
4407        BR(1)=PRODR/DEN
4408        BI(1)=PRODI/DEN
4409 50   CONTINUE
4410      RETURN
4411C----------------------- END OF SUBROUTINE SOLBC ------------------------
4412      END
4413c
4414C
4415      subroutine Elmhes(nm,n,low,igh,a,int)
4416C
4417      integer i,j,m,n,la,nm,igh,kp1,low,mm1,mp1
4418      real*8 a(nm,n)
4419      real*8 x,y
4420      real*8 dabs
4421      integer int(igh)
4422C
4423C     this subroutine is a translation of the algol procedure elmhes,
4424C     num. math. 12, 349-368(1968) by martin and wilkinson.
4425C     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
4426C
4427C     given a real general matrix, this subroutine
4428C     reduces a submatrix situated in rows and columns
4429C     low through igh to upper hessenberg form by
4430C     stabilized elementary similarity transformations.
4431C
4432C     on input:
4433C
4434C      nm must be set to the row dimension of two-dimensional
4435C        array parameters as declared in the calling program
4436C        dimension statement;
4437C
4438C      n is the order of the matrix;
4439C
4440C      low and igh are integers determined by the balancing
4441C        subroutine  balanc.      if  balanc  has not been used,
4442C        set low=1, igh=n;
4443C
4444C      a contains the input matrix.
4445C
4446C     on output:
4447C
4448C      a contains the hessenberg matrix.  the multipliers
4449C        which were used in the reduction are stored in the
4450C        remaining triangle under the hessenberg matrix;
4451C
4452C      int contains information on the rows and columns
4453C        interchanged in the reduction.
4454C        only elements low through igh are used.
4455C
4456C     questions and comments should be directed to b. s. garbow,
4457C     applied mathematics division, argonne national laboratory
4458C
4459C     ------------------------------------------------------------------
4460C
4461      la = igh - 1
4462      kp1 = low + 1
4463      if (la .lt. kp1) go to 200
4464C
4465      do 180 m = kp1, la
4466       mm1 = m - 1
4467       x = 0.0d0
4468       i = m
4469C
4470       do 100 j = m, igh
4471          if (dabs(a(j,mm1)) .le. dabs(x)) go to 100
4472          x = a(j,mm1)
4473          i = j
4474  100   continue
4475C
4476       int(m) = i
4477       if (i .eq. m) go to 130
4478C    :::::::::: interchange rows and columns of a ::::::::::
4479       do 110 j = mm1, n
4480          y = a(i,j)
4481          a(i,j) = a(m,j)
4482          a(m,j) = y
4483  110   continue
4484C
4485       do 120 j = 1, igh
4486          y = a(j,i)
4487          a(j,i) = a(j,m)
4488          a(j,m) = y
4489  120   continue
4490C    :::::::::: end interchange ::::::::::
4491  130   if (x .eq. 0.0d0) go to 180
4492       mp1 = m + 1
4493C
4494       do 160 i = mp1, igh
4495          y = a(i,mm1)
4496          if (y .eq. 0.0d0) go to 160
4497          y = y / x
4498          a(i,mm1) = y
4499C
4500          do 140 j = m, n
4501  140      a(i,j) = a(i,j) - y * a(m,j)
4502C
4503          do 150 j = 1, igh
4504  150      a(j,m) = a(j,m) + y * a(j,i)
4505C
4506  160   continue
4507C
4508  180 continue
4509C
4510  200 return
4511C    :::::::::: last card of elmhes ::::::::::
4512      end
4513
4514        SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN)
4515C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5"
4516        IMPLICIT KPP_REAL (A-H,O-Z)
4517        DIMENSION Y(N),CONT(LRC)
4518        COMMON /INTERN/XOUT
4519        IF (NR.EQ.1) THEN
4520           WRITE (6,99) X,Y(1),Y(2),NR-1
4521           XOUT=0.2D0
4522        ELSE
4523 10        CONTINUE
4524           IF (X.GE.XOUT) THEN
4525C --- CONTINUOUS OUTPUT FOR RADAU5
4526              WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC),
4527     &                     CONTR5(2,XOUT,CONT,LRC),NR-1
4528              XOUT=XOUT+0.2D0
4529              GOTO 10
4530           END IF
4531        END IF
4532 99     FORMAT(1X,'X =',F5.2,'    Y =',2E18.10,'    NSTEP =',I4)
4533        RETURN
4534        END
4535
4536        SUBROUTINE FUNC_CHEM (N, T, V, FCT, RPAR, IPAR)
4537        IMPLICIT NONE
4538        INCLUDE 'KPP_ROOT_Parameters.h'
4539        INCLUDE 'KPP_ROOT_Global.h'
4540        KPP_REAL V(NVAR), FCT(NVAR)
4541        KPP_REAL T, TOLD, RPAR(1)
4542        INTEGER N, IPAR(1)
4543        TOLD = TIME
4544        TIME = T
4545        CALL Update_SUN()
4546        CALL Update_RCONST()
4547        TIME = TOLD
4548        CALL Fun(V,  FIX, RCONST, FCT)
4549        RETURN
4550        END
4551
4552        SUBROUTINE JAC_CHEM (N, T, V, JV, NROWPD, RPAR, IPAR)
4553        IMPLICIT NONE
4554        INCLUDE 'KPP_ROOT_Parameters.h'
4555        INCLUDE 'KPP_ROOT_Global.h'
4556        KPP_REAL V(NVAR), JV(NVAR,NVAR)
4557        KPP_REAL T, TOLD, RPAR(1)
4558        INTEGER N, IPAR(1), NROWPD, i, j
4559        TOLD = TIME
4560        TIME = T
4561        CALL Update_SUN()
4562        CALL Update_RCONST()
4563        TIME = TOLD
4564        DO i=1,NVAR
4565           DO j=1,NVAR
4566              JV(i,j) = 0.D0
4567           END DO
4568        END DO
4569        CALL Jac(V,  FIX, RCONST, JV)
4570        RETURN
4571        END
4572
Note: See TracBrowser for help on using the repository browser.