source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/kpp_ros4.f @ 2696

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

Merge of branch palm4u into trunk

File size: 36.0 KB
Line 
1      SUBROUTINE INTEGRATE( TIN, TOUT )
2
3      IMPLICIT KPP_REAL (A-H,O-Z)       
4      INCLUDE 'KPP_ROOT_params.h'
5      INCLUDE 'KPP_ROOT_global.h'
6
7C TIN - Start Time
8      KPP_REAL TIN
9C TOUT - End Time
10      KPP_REAL TOUT
11      INTEGER i
12
13      PARAMETER (LWORK=2*NVAR*NVAR+14*NVAR+20,LIWORK=NVAR+20)
14      KPP_REAL WORK(LWORK)
15      INTEGER IWORK(LIWORK)
16      EXTERNAL FUNC_CHEM,JAC_CHEM,SOLOUT
17
18      ITOL=1     ! --- VECTOR TOLERANCES
19      IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
20      MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
21      MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
22      IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
23      IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
24      IDFX=0     ! --- INTERNAL TIME DERIVATIVE
25
26      DO i=1,20
27        IWORK(i) = 0
28        WORK(i) = 0.D0
29      ENDDO
30      IWORK(2) = 6
31
32      CALL KPP_ROS4(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT,
33     &                  STEPMIN,RTOL,ATOL,ITOL,
34     &                  JAC_CHEM,IJAC,MLJAC,MUJAC,FUNC_CHEM,IDFX,
35     &                  FUNC_CHEM,IMAS,MLJAC,MUJAC,
36     &                  SOLOUT,IOUT,
37     &                  WORK,LWORK,IWORK,LIWORK,IDID)
38
39      IF (IDID.LT.0) THEN
40        print *,'KPP_ROS4: Unsucessful step at T=',
41     &          TIN,' (IDID=',IDID,')'
42      ENDIF
43
44      RETURN
45      END
46
47
48      SUBROUTINE KPP_ROS4(N,FCN,IFCN,X,Y,XEND,H,
49     &                  RelTol,AbsTol,ITOL,
50     &                  JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX,
51     &                  MAS ,IMAS,MLMAS,MUMAS,
52     &                  SOLOUT,IOUT,
53     &                  WORK,LWORK,IWORK,LIWORK,IDID)
54C ----------------------------------------------------------
55C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
56C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
57C     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4 
58C     (WITH STEP SIZE CONTROL).
59C     C.F. SECTION IV.7
60C
61C     AUTHORS: E. HAIRER AND G. WANNER
62C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
63C              CH-1211 GENEVE 24, SWITZERLAND
64C              E-MAIL:  HAIRER@CGEUGE51.BITNET,  WANNER@CGEUGE51.BITNET
65C     
66C     THIS CODE IS PART OF THE BOOK:
67C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
68C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
69C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
70C         SPRINGER-VERLAG (1990)               
71C     
72C     VERSION OF NOVEMBER 17, 1992
73C
74C     INPUT PARAMETERS 
75C     ---------------- 
76C     N           DIMENSION OF THE SYSTEM
77C
78C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
79C                 VALUE OF F(X,Y):
80C                    SUBROUTINE FCN(N,X,Y,F)
81C                    KPP_REAL X,Y(N),F(N)
82C                    F(1)=...   ETC.
83C
84C     IFCN        GIVES INFORMATION ON FCN:
85C                    IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS)
86C                    IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS)
87C
88C     X           INITIAL X-VALUE
89C
90C     Y(N)        INITIAL VALUES FOR Y
91C
92C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
93C
94C     H           INITIAL STEP SIZE GUESS;
95C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
96C                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
97C                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
98C                 ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
99C                 STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
100C                 (IF H=0.D0, THE CODE PUTS H=1.D-6).
101C
102C     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
103C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
104C
105C     ITOL        SWITCH FOR RelTol AND AbsTol:
106C                   ITOL=0: BOTH RelTol AND AbsTol ARE SCALARS.
107C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
108C                     Y(I) BELOW RelTol*ABS(Y(I))+AbsTol
109C                   ITOL=1: BOTH RelTol AND AbsTol ARE VECTORS.
110C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
111C                     RelTol(I)*ABS(Y(I))+AbsTol(I).
112C
113C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
114C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
115C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
116C                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
117C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
118C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
119C                    KPP_REAL X,Y(N),DFY(LDFY,N)
120C                    DFY(1,1)= ...
121C                 LDFY, THE COLOMN-LENGTH OF THE ARRAY, IS
122C                 FURNISHED BY THE CALLING PROGRAM.
123C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
124C                    BE FULL AND THE PARTIAL DERIVATIVES ARE
125C                    STORED IN DFY AS
126C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
127C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
128C                    THE PARTIAL DERIVATIVES ARE STORED
129C                    DIAGONAL-WISE AS
130C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
131C
132C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
133C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
134C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
135C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
136C
137C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
138C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
139C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
140C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN 
141C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
142C                       THE MAIN DIAGONAL).
143C
144C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
145C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
146C                 NEED NOT BE DEFINED IF MLJAC=N.
147C
148C     DFX         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
149C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X
150C                 (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1;
151C                 SUPPLY A DUMMY SUBROUTINE IN THE CASE IDFX=0 OR IFCN=0).
152C                 OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM
153C                    SUBROUTINE DFX(N,X,Y,FX)
154C                    KPP_REAL X,Y(N),FX(N)
155C                    FX(1)= ...
156C               
157C     IDFX        SWITCH FOR THE COMPUTATION OF THE DF/DX:
158C                    IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE
159C                       DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED.
160C                    IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX.
161C
162C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
163C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
164C
165C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
166C                 MATRIX M.
167C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
168C                 MATRIX AND NEEDS NOT TO BE DEFINED;
169C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
170C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
171C                    SUBROUTINE MAS(N,AM,LMAS)
172C                    KPP_REAL AM(LMAS,N)
173C                    AM(1,1)= ....
174C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
175C                    AS FULL MATRIX LIKE
176C                         AM(I,J) = M(I,J)
177C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
178C                    DIAGONAL-WISE AS
179C                         AM(I-J+MUMAS+1,J) = M(I,J).
180C
181C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
182C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
183C                       MATRIX, MAS IS NEVER CALLED.
184C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
185C
186C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
187C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
188C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
189C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
190C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
191C                       THE MAIN DIAGONAL).
192C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
193C
194C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
195C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
196C                 NEED NOT BE DEFINED IF MLMAS=N.
197C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
198C
199C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
200C                 NUMERICAL SOLUTION DURING INTEGRATION. 
201C                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
202C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0. 
203C                 IT MUST HAVE THE FORM
204C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
205C                    KPP_REAL X,Y(N)
206C                    .... 
207C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
208C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
209C                    THE FIRST GRID-POINT).
210C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
211C                    IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM.
212C           
213C     IOUT        GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
214C                    IOUT=0: SUBROUTINE IS NEVER CALLED
215C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
216C
217C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
218C                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
219C                 "LWORK" MUST BE AT LEAST
220C                             N*(LJAC+LMAS+LE1+8)+5
221C                 WHERE
222C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
223C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
224C                 AND                 
225C                    LMAS=0              IF IMAS=0
226C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
227C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
228C                 AND
229C                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
230C                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.). 
231c
232C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
233C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
234C                 STORAGE REQUIREMENT IS 
235C                             LWORK = 2*N*N+8*N+5.
236C
237C     LWORK       DECLARED LENGHT OF ARRAY "WORK".
238C
239C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
240C                 "LIWORK" MUST BE AT LEAST N+2.
241C
242C     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
243C
244C ----------------------------------------------------------------------
245C 
246C     SOPHISTICATED SETTING OF PARAMETERS
247C     -----------------------------------
248C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK 
249C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5)
250C              AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO.
251C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
252C
253C    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
254C              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
255C
256C    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
257C              IF IWORK(2).EQ.1  METHOD OF SHAMPINE
258C              IF IWORK(2).EQ.2  METHOD GRK4T OF KAPS-RENTROP
259C              IF IWORK(2).EQ.3  METHOD GRK4A OF KAPS-RENTROP
260C              IF IWORK(2).EQ.4  METHOD OF VAN VELDHUIZEN (GAMMA=1/2)
261C              IF IWORK(2).EQ.5  METHOD OF VAN VELDHUIZEN ("D-STABLE")
262C              IF IWORK(2).EQ.6  AN L-STABLE METHOD
263C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2.
264C
265C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
266C
267C    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
268C
269C    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
270C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
271C                 WORK(3) <= HNEW/HOLD <= WORK(4)
272C              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0
273C
274C    WORK(5)   AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS
275C              THE STEP SIZE IS MULTIPLIED BY WORK(5)
276C              DEFAULT VALUES: WORK(5)=0.1D0
277C
278C-----------------------------------------------------------------------
279C
280C     OUTPUT PARAMETERS 
281C     ----------------- 
282C     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
283C                 (AFTER SUCCESSFUL RETURN X=XEND)
284C
285C     Y(N)        SOLUTION AT X
286C 
287C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
288C
289C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
290C                   IDID=1  COMPUTATION SUCCESSFUL,
291C                   IDID=-1 COMPUTATION UNSUCCESSFUL.
292C
293C --------------------------------------------------------- 
294C *** *** *** *** *** *** *** *** *** *** *** *** ***
295C          DECLARATIONS 
296C *** *** *** *** *** *** *** *** *** *** *** *** ***
297      IMPLICIT KPP_REAL (A-H,O-Z)
298      DIMENSION Y(N),AbsTol(1),RelTol(1),WORK(LWORK),IWORK(LIWORK)
299      LOGICAL AUTNMS,IMPLCT,JBAND,ARRET
300      EXTERNAL FCN,JAC,DFX,MAS,SOLOUT
301      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
302C -----------------------------------------------------
303C --- COMMON STAT CAN BE USED FOR STATISTICS
304C ---    NFCN      NUMBER OF FUNC_CHEMCTION EVALUATIONS (THOSE FOR NUMERICAL
305C                  EVALUATION OF THE JACOBIAN ARE NOT COUNTED) 
306C ---    NJAC      NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
307C                  OR NUMERICALLY)
308C ---    NSTEP     NUMBER OF COMPUTED STEPS
309C ---    NACCPT    NUMBER OF ACCEPTED STEPS
310C ---    NREJCT    NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
311C                  HAS BEEN ACCEPTED)
312C ---    NDEC      NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
313C ---    NSOL      NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
314C -----------------------------------------------------
315C *** *** *** *** *** *** ***
316C        SETTING THE PARAMETERS 
317C *** *** *** *** *** *** ***
318      NFCN=0
319      NJAC=0
320      NSTEP=0
321      NACCPT=0
322      NREJCT=0
323      NDEC=0
324      NSOL=0
325      ARRET=.FALSE.
326C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
327      IF(IWORK(1).EQ.0)THEN
328         NMAX=100000
329      ELSE
330         NMAX=IWORK(1)
331         IF(NMAX.LE.0)THEN
332            WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1)
333            ARRET=.TRUE.
334         END IF
335      END IF
336C -------- METH   COEFFICIENTS OF THE METHOD
337      IF(IWORK(2).EQ.0)THEN
338         METH=2
339      ELSE
340         METH=IWORK(2)
341         IF(METH.LE.0.OR.METH.GE.7)THEN
342            WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2)
343            ARRET=.TRUE.
344         END IF
345      END IF
346C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 
347      IF(WORK(1).EQ.0.D0)THEN
348         UROUND=1.D-16
349      ELSE
350         UROUND=WORK(1)
351         IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN
352            WRITE(6,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1)
353            ARRET=.TRUE.
354         END IF
355      END IF
356C -------- MAXIMAL STEP SIZE
357      IF(WORK(2).EQ.0.D0)THEN
358         HMAX=XEND-X
359      ELSE
360         HMAX=WORK(2)
361      END IF
362C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
363      IF(WORK(3).EQ.0.D0)THEN
364         FAC1=5.D0
365      ELSE
366         FAC1=1.D0/WORK(3)
367      END IF
368      IF(WORK(4).EQ.0.D0)THEN
369         FAC2=1.D0/6.0D0
370      ELSE
371         FAC2=1.D0/WORK(4)
372      END IF
373C -------  FACREJ    FOR THE HUMP
374      IF(WORK(5).EQ.0.D0)THEN
375         FACREJ=0.1D0
376      ELSE
377         FACREJ=WORK(5)
378      END IF
379C --------- CHECK IF TOLERANCES ARE O.K.
380      IF (ITOL.EQ.0) THEN
381          IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
382              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
383              ARRET=.TRUE.
384          END IF
385      ELSE
386          DO 15 I=1,N
387          IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
388              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
389              ARRET=.TRUE.
390          END IF
391  15      CONTINUE
392      END IF
393C *** *** *** *** *** *** *** *** *** *** *** *** ***
394C         COMPUTATION OF ARRAY ENTRIES
395C *** *** *** *** *** *** *** *** *** *** *** *** ***
396C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
397      AUTNMS=IFCN.EQ.0
398      IMPLCT=IMAS.NE.0
399      JBAND=MLJAC.NE.N
400      ARRET=.FALSE.
401C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
402C -- JACOBIAN 
403      IF(JBAND)THEN
404         LDJAC=MLJAC+MUJAC+1
405      ELSE
406         LDJAC=N
407      END IF
408C -- MATRIX E FOR LINEAR ALGEBRA
409      IF(JBAND)THEN
410         LDE=2*MLJAC+MUJAC+1
411      ELSE
412         LDE=N
413      END IF
414C -- MASS MATRIX
415      IF (IMPLCT) THEN
416          IF (MLMAS.NE.N) THEN
417              LDMAS=MLMAS+MUMAS+1
418          ELSE
419              LDMAS=N
420          END IF
421C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC"
422          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
423              WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF
424     & "JAC"'
425            ARRET=.TRUE.
426          END IF
427      ELSE
428          LDMAS=0
429      END IF
430      LDMAS2=MAX(1,LDMAS)
431C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
432      IEYNEW=6
433      IEDY1=IEYNEW+N
434      IEDY=IEDY1+N
435      IEAK1=IEDY+N
436      IEAK2=IEAK1+N
437      IEAK3=IEAK2+N
438      IEAK4=IEAK3+N
439      IEFX =IEAK4+N
440      IEJAC=IEFX +N
441      IEMAS=IEJAC+N*LDJAC
442      IEE  =IEMAS+N*LDMAS
443C ------ TOTAL STORAGE REQUIREMENT -----------
444      ISTORE=IEE+N*LDE-1
445      IF(ISTORE.GT.LWORK)THEN
446         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
447         ARRET=.TRUE.
448      END IF
449C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
450      IEIP=3
451C --------- TOTAL REQUIREMENT ---------------
452      ISTORE=IEIP+N-1
453      IF(ISTORE.GT.LIWORK)THEN
454         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
455         ARRET=.TRUE.
456      END IF
457C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
458      IF (ARRET) THEN
459         IDID=-1
460         RETURN
461      END IF
462C -------- CALL TO CORE INTEGRATOR ------------
463      CALL RO4COR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,IJAC,
464     &   MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
465     &   NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS,IMPLCT,JBAND,
466     &   LDJAC,LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),
467     &   WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),
468     &   WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP))
469C ----------- RETURN -----------
470      RETURN
471      END
472C
473C
474C
475C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
476C
477      SUBROUTINE RO4COR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,
478     &  IJAC,MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
479     &  NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS,IMPLCT,BANDED,
480     &  LFJAC,LE,LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,FMAS,IP)
481C ----------------------------------------------------------
482C     CORE INTEGRATOR FOR ROS4
483C     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED 
484C ---------------------------------------------------------- 
485C         DECLARATIONS 
486C ---------------------------------------------------------- 
487      IMPLICIT KPP_REAL (A-H,O-Z)       
488      INCLUDE 'KPP_ROOT_params.h'
489      INCLUDE 'KPP_ROOT_sparse.h'
490      KPP_REAL Y(N),YNEW(N),DY1(N),DY(N),AK1(N),
491     *  AK2(N),AK3(N),AK4(N),FX(N),
492     *  FJAC(LU_NONZERO),E(LU_NONZERO),
493     *  FMAS(LDMAS,N),AbsTol(1),RelTol(1)
494      INTEGER IP(N)
495      LOGICAL REJECT,RJECT2,AUTNMS,IMPLCT,BANDED
496      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
497      EXTERNAL FCN,JAC,MAS,SOLOUT,DFX
498
499
500C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
501      IF (IMPLCT) CALL MAS(N,FMAS,LDMAS)
502C ---- PREPARE BANDWIDTHS -----
503      IF (BANDED) THEN
504          MLE=MLJAC
505          MUE=MUJAC
506          MBJAC=MLJAC+MUJAC+1
507          MBB=MLMAS+MUMAS+1
508          MDIAG=MLE+MUE+1
509          MBDIAG=MUMAS+1
510          MDIFF=MLE+MUE-MUMAS
511      END IF
512C *** *** *** *** *** *** ***
513C  INITIALISATIONS
514C *** *** *** *** *** *** ***
515      IF (METH.EQ.1) CALL SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43,
516     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
517      IF (METH.EQ.2) CALL GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43,
518     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
519      IF (METH.EQ.3) CALL GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43,
520     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
521      IF (METH.EQ.4) CALL VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43,
522     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
523      IF (METH.EQ.5) CALL VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43,
524     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
525      IF (METH.EQ.6) CALL LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43,
526     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
527C --- INITIAL PREPARATIONS
528      POSNEG=SIGN(1.D0,XEND-X)
529      HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
530      IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
531      H=MIN(ABS(H),HMAXN) 
532      H=SIGN(H,POSNEG) 
533      REJECT=.FALSE.
534      NSING=0
535      IRTRN=1 
536      XXOLD=X
537      IF (IRTRN.LT.0) GOTO 79
538C --- BASIC INTEGRATION STEP 
539   1  IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
540      IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
541          H=HOPT
542          IDID=1
543          RETURN
544      END IF
545      HOPT=H
546      IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X
547      CALL FCN(N,X,Y,DY1)
548      NFCN=NFCN+1
549
550C *** *** *** *** *** *** ***
551C  COMPUTATION OF THE JACOBIAN
552C *** *** *** *** *** *** ***
553      NJAC=NJAC+1
554      CALL JAC(N,X,Y,FJAC)
555
556      IF (.NOT.AUTNMS) THEN
557          IF (IDFX.EQ.0) THEN
558C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X
559              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(X)))
560              XDELT=X+DELT
561              CALL FCN(N,XDELT,Y,AK1)
562              DO 19 J=1,N
563  19          FX(J)=(AK1(J)-DY1(J))/DELT
564          ELSE
565C --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X
566              CALL DFX(N,X,Y,FX)
567          END IF
568      END IF
569   2  CONTINUE
570
571C *** *** *** *** *** *** ***
572C  COMPUTE THE STAGES
573C *** *** *** *** *** *** ***
574      NDEC=NDEC+1
575      HC21=C21/H
576      HC31=C31/H
577      HC32=C32/H
578      HC41=C41/H
579      HC42=C42/H
580      HC43=C43/H
581      FAC=1.D0/(H*GAMMA)
582C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
583              DO 800 I=1,LU_NONZERO
584  800         E(I)=-FJAC(I)
585              DO 801 J=1,N
586  801         E(LU_DIAG(J))=E(LU_DIAG(J))+FAC
587              CALL KppDecomp (E,IER)
588              IF (IER.NE.0) GOTO 80
589              IF (AUTNMS) THEN
590C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
591C ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
592C ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
593C ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
594                  DO 803 I=1,N
595  803             AK1(I)=DY1(I)
596                  CALL KppSolve(E,AK1)
597                  DO 810 I=1,N
598  810             YNEW(I)=Y(I)+A21*AK1(I) 
599                  CALL FCN(N,X,YNEW,DY)
600                  DO 811 I=1,N
601  811             AK2(I)=DY(I)+HC21*AK1(I)
602                  CALL KppSolve(E,AK2)
603                  DO 820 I=1,N
604  820             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) 
605                  CALL FCN(N,X,YNEW,DY)
606                  DO 821 I=1,N
607  821             AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
608                  CALL KppSolve(E,AK3)
609                  DO 831 I=1,N
610  831             AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
611                  CALL KppSolve(E,AK4)
612              ELSE
613C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
614C ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
615C ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
616C ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
617                  HD1=H*D1
618                  HD2=H*D2
619                  HD3=H*D3
620                  HD4=H*D4
621                  DO 903 I=1,N
622  903             AK1(I)=DY1(I)+HD1*FX(I)
623                  CALL KppSolve(E,AK1)
624                  DO 910 I=1,N
625  910             YNEW(I)=Y(I)+A21*AK1(I) 
626                  CALL FCN(N,X+C2*H,YNEW,DY)
627                  DO 911 I=1,N
628  911             AK2(I)=DY(I)+HD2*FX(I)+HC21*AK1(I)
629                  CALL KppSolve(E,AK2)
630                  DO 920 I=1,N
631  920             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) 
632                  CALL FCN(N,X+C3*H,YNEW,DY)
633                  DO 921 I=1,N
634  921             AK3(I)=DY(I)+HD3*FX(I)+HC31*AK1(I)+HC32*AK2(I)
635                  CALL KppSolve(E,AK3)
636                  DO 931 I=1,N
637  931             AK4(I)=DY(I)+HD4*FX(I)+HC41*AK1(I)+HC42*AK2(I)
638     &                  +HC43*AK3(I) 
639                  CALL KppSolve(E,AK4)
640              END IF
641      NSOL=NSOL+4
642      NFCN=NFCN+2
643C *** *** *** *** *** *** ***
644C  ERROR ESTIMATION 
645C *** *** *** *** *** *** ***
646      NSTEP=NSTEP+1
647C ------------ NEW SOLUTION ---------------
648      DO 240 I=1,N
649  240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)+B4*AK4(I) 
650C ------------ COMPUTE ERROR ESTIMATION ----------------
651      ERR=0.D0
652      DO 300 I=1,N
653      S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)+E4*AK4(I) 
654      IF (ITOL.EQ.0) THEN
655         SK=AbsTol(1)+RelTol(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
656      ELSE
657         SK=AbsTol(I)+RelTol(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
658      END IF
659  300 ERR=ERR+(S/SK)**2
660      ERR=DSQRT(ERR/N)
661C --- COMPUTATION OF HNEW
662C --- WE REQUIRE .2<=HNEW/H<=6.
663      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0))
664      HNEW=H/FAC 
665C *** *** *** *** *** *** ***
666C  IS THE ERROR SMALL ENOUGH ?
667C *** *** *** *** *** *** ***
668      IF (ERR.LE.1.D0) THEN
669C --- STEP IS ACCEPTED 
670         NACCPT=NACCPT+1
671         DO 44 I=1,N 
672  44     Y(I)=YNEW(I) 
673         XXOLD=X 
674         X=X+H
675         IF (IRTRN.LT.0) GOTO 79
676         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
677         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) 
678         REJECT=.FALSE.
679         RJECT2=.FALSE.
680         H=HNEW
681         GOTO 1
682      ELSE
683C --- STEP IS REJECTED 
684         IF (RJECT2) HNEW=H*FACREJ
685         IF (REJECT) RJECT2=.TRUE.
686         REJECT=.TRUE.
687         H=HNEW
688         IF (NACCPT.GE.1) NREJCT=NREJCT+1
689         GOTO 2
690      END IF
691C --- EXIT
692  80  WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
693      NSING=NSING+1
694      IF (NSING.GE.5) GOTO 79
695      H=H*0.5D0
696      GOTO 2
697  79  WRITE(6,979)X,H
698 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7)
699      IDID=-1
700      RETURN
701      END
702C
703      SUBROUTINE SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43,
704     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
705      IMPLICIT KPP_REAL (A-H,O-Z)
706         A21=2.D0
707         A31=48.D0/25.D0
708         A32=6.D0/25.D0
709         C21=-8.D0
710         C31=372.D0/25.D0
711         C32=12.D0/5.D0
712         C41=-112.D0/125.D0
713         C42=-54.D0/125.D0
714         C43=-2.D0/5.D0
715         B1=19.D0/9.D0
716         B2=1.D0/2.D0
717         B3=25.D0/108.D0
718         B4=125.D0/108.D0
719         E1=17.D0/54.D0
720         E2=7.D0/36.D0
721         E3=0.D0
722         E4=125.D0/108.D0
723         GAMMA=.5D0
724         C2= 0.1000000000000000D+01
725         C3= 0.6000000000000000D+00
726         D1= 0.5000000000000000D+00
727         D2=-0.1500000000000000D+01
728         D3= 0.2420000000000000D+01
729         D4= 0.1160000000000000D+00
730      RETURN
731      END
732C
733      SUBROUTINE GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43,
734     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
735      IMPLICIT KPP_REAL (A-H,O-Z)
736       A21= 0.1108860759493671D+01
737       A31= 0.2377085261983360D+01
738       A32= 0.1850114988899692D+00
739       C21=-0.4920188402397641D+01
740       C31= 0.1055588686048583D+01
741       C32= 0.3351817267668938D+01
742       C41= 0.3846869007049313D+01
743       C42= 0.3427109241268180D+01
744       C43=-0.2162408848753263D+01
745       B1= 0.1845683240405840D+01
746       B2= 0.1369796894360503D+00
747       B3= 0.7129097783291559D+00
748       B4= 0.6329113924050632D+00
749       E1= 0.4831870177201765D-01
750       E2=-0.6471108651049505D+00
751       E3= 0.2186876660500240D+00
752       E4=-0.6329113924050632D+00
753       GAMMA= 0.3950000000000000D+00
754       C2= 0.4380000000000000D+00
755       C3= 0.8700000000000000D+00
756       D1= 0.3950000000000000D+00
757       D2=-0.3726723954840920D+00
758       D3= 0.6629196544571492D-01
759       D4= 0.4340946962568634D+00
760      RETURN
761      END
762C
763      SUBROUTINE GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43,
764     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
765      IMPLICIT KPP_REAL (A-H,O-Z)
766       A21= 0.2000000000000000D+01
767       A31= 0.4524708207373116D+01
768       A32= 0.4163528788597648D+01
769       C21=-0.5071675338776316D+01
770       C31= 0.6020152728650786D+01
771       C32= 0.1597506846727117D+00
772       C41=-0.1856343618686113D+01
773       C42=-0.8505380858179826D+01
774       C43=-0.2084075136023187D+01
775       B1= 0.3957503746640777D+01
776       B2= 0.4624892388363313D+01
777       B3= 0.6174772638750108D+00
778       B4= 0.1282612945269037D+01
779       E1= 0.2302155402932996D+01
780       E2= 0.3073634485392623D+01
781       E3=-0.8732808018045032D+00
782       E4=-0.1282612945269037D+01
783       GAMMA= 0.2310000000000000D+00
784       C2= 0.4620000000000000D+00
785       C3= 0.8802083333333334D+00
786       D1= 0.2310000000000000D+00
787       D2=-0.3962966775244303D-01
788       D3= 0.5507789395789127D+00
789       D4=-0.5535098457052764D-01
790      RETURN
791      END
792C
793      SUBROUTINE VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43,
794     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
795C --- METHOD GIVEN BY VAN VELDHUIZEN
796      IMPLICIT KPP_REAL (A-H,O-Z)
797       A21= 0.2000000000000000D+01
798       A31= 0.1750000000000000D+01
799       A32= 0.2500000000000000D+00
800       C21=-0.8000000000000000D+01
801       C31=-0.8000000000000000D+01
802       C32=-0.1000000000000000D+01
803       C41= 0.5000000000000000D+00
804       C42=-0.5000000000000000D+00
805       C43= 0.2000000000000000D+01
806       B1= 0.1333333333333333D+01
807       B2= 0.6666666666666667D+00
808       B3=-0.1333333333333333D+01
809       B4= 0.1333333333333333D+01
810       E1=-0.3333333333333333D+00
811       E2=-0.3333333333333333D+00
812       E3=-0.0000000000000000D+00
813       E4=-0.1333333333333333D+01
814       GAMMA= 0.5000000000000000D+00
815       C2= 0.1000000000000000D+01
816       C3= 0.5000000000000000D+00
817       D1= 0.5000000000000000D+00
818       D2=-0.1500000000000000D+01
819       D3=-0.7500000000000000D+00
820       D4= 0.2500000000000000D+00
821      RETURN
822      END
823C
824      SUBROUTINE VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43,
825     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
826C --- METHOD GIVEN BY VAN VELDHUIZEN
827      IMPLICIT KPP_REAL (A-H,O-Z)
828       A21= 0.2000000000000000D+01
829       A31= 0.4812234362695436D+01
830       A32= 0.4578146956747842D+01
831       C21=-0.5333333333333331D+01
832       C31= 0.6100529678848254D+01
833       C32= 0.1804736797378427D+01
834       C41=-0.2540515456634749D+01
835       C42=-0.9443746328915205D+01
836       C43=-0.1988471753215993D+01
837       B1= 0.4289339254654537D+01
838       B2= 0.5036098482851414D+01
839       B3= 0.6085736420673917D+00
840       B4= 0.1355958941201148D+01
841       E1= 0.2175672787531755D+01
842       E2= 0.2950911222575741D+01
843       E3=-0.7859744544887430D+00
844       E4=-0.1355958941201148D+01
845       GAMMA= 0.2257081148225682D+00
846       C2= 0.4514162296451364D+00
847       C3= 0.8755928946018455D+00
848       D1= 0.2257081148225682D+00
849       D2=-0.4599403502680582D-01
850       D3= 0.5177590504944076D+00
851       D4=-0.3805623938054428D-01
852      RETURN
853      END
854C
855      SUBROUTINE LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43,
856     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
857C --- AN L-STABLE METHOD
858      IMPLICIT KPP_REAL (A-H,O-Z)
859       A21= 0.2000000000000000D+01
860       A31= 0.1867943637803922D+01
861       A32= 0.2344449711399156D+00
862       C21=-0.7137615036412310D+01
863       C31= 0.2580708087951457D+01
864       C32= 0.6515950076447975D+00
865       C41=-0.2137148994382534D+01
866       C42=-0.3214669691237626D+00
867       C43=-0.6949742501781779D+00
868C  B_i = M_i
869       B1= 0.2255570073418735D+01
870       B2= 0.2870493262186792D+00
871       B3= 0.4353179431840180D+00
872       B4= 0.1093502252409163D+01
873C E_i = error estimator       
874       E1=-0.2815431932141155D+00
875       E2=-0.7276199124938920D-01
876       E3=-0.1082196201495311D+00
877       E4=-0.1093502252409163D+01
878C gamma = gamma
879       GAMMA= 0.5728200000000000D+00
880C  C_i = alpha_i
881       C2= 0.1145640000000000D+01
882       C3= 0.6552168638155900D+00
883C D_i = gamma_i       
884       D1= 0.5728200000000000D+00
885       D2=-0.1769193891319233D+01
886       D3= 0.7592633437920482D+00
887       D4=-0.1049021087100450D+00
888      RETURN
889      END
890
891        SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
892C --- PRINTS SOLUTION
893        IMPLICIT KPP_REAL (A-H,O-Z)
894        DIMENSION Y(N)
895        COMMON /INTERN/XOUT
896        IF (NR.EQ.1) THEN
897           WRITE (6,99) X,Y(1),Y(2),NR-1
898           XOUT=0.1D0
899        ELSE
900           IF (X.GE.XOUT) THEN
901              WRITE (6,99) X,Y(1),Y(2),NR-1
902              XOUT=DMAX1(XOUT+0.1D0,X)
903           END IF
904        END IF
905 99     FORMAT(1X,'X =',F5.2,'    Y =',2E18.10,'    NSTEP =',I4)
906        RETURN
907        END
908
909      SUBROUTINE DEC (N, NDIM, A, IP, IER)
910C VERSION REAL KPP_REAL
911      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
912      KPP_REAL A,T
913      DIMENSION A(NDIM,N), IP(N)
914C-----------------------------------------------------------------------
915C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
916C  INPUT..
917C     N = ORDER OF MATRIX.
918C     NDIM = DECLARED DIMENSION OF ARRAY  A .
919C     A = MATRIX TO BE TRIANGULARIZED.
920C  OUTPUT..
921C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
922C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
923C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
924C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
925C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
926C           SINGULAR AT STAGE K.
927C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
928C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
929C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
930C
931C  REFERENCE..
932C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
933C     C.A.C.M. 15 (1972), P. 274.
934C-----------------------------------------------------------------------
935      IER = 0
936      IP(N) = 1
937      IF (N .EQ. 1) GO TO 70
938      NM1 = N - 1
939      DO 60 K = 1,NM1
940        KP1 = K + 1
941        M = K
942        DO 10 I = KP1,N
943          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I 
944 10     CONTINUE
945        IP(K) = M
946        T = A(M,K)
947        IF (M .EQ. K) GO TO 20
948        IP(N) = -IP(N)
949        A(M,K) = A(K,K)
950        A(K,K) = T
951 20     CONTINUE
952        IF (T .EQ. 0.D0) GO TO 80
953        T = 1.D0/T
954        DO 30 I = KP1,N
955 30       A(I,K) = -A(I,K)*T
956        DO 50 J = KP1,N
957          T = A(M,J)
958          A(M,J) = A(K,J)
959          A(K,J) = T
960          IF (T .EQ. 0.D0) GO TO 45
961          DO 40 I = KP1,N
962 40         A(I,J) = A(I,J) + A(I,K)*T
963 45       CONTINUE
964 50       CONTINUE
965 60     CONTINUE
966 70   K = N
967      IF (A(N,N) .EQ. 0.D0) GO TO 80
968      RETURN
969 80   IER = K
970      IP(N) = 0
971      RETURN
972C----------------------- END OF SUBROUTINE DEC -------------------------
973      END
974C
975C
976      SUBROUTINE SOL (N, NDIM, A, B, IP)
977C VERSION REAL KPP_REAL
978      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
979      KPP_REAL A,B,T
980      DIMENSION A(NDIM,N), B(N), IP(N)
981C-----------------------------------------------------------------------
982C  SOLUTION OF LINEAR SYSTEM, A*X = B .
983C  INPUT..
984C    N = ORDER OF MATRIX.
985C    NDIM = DECLARED DIMENSION OF ARRAY  A .
986C    A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
987C    B = RIGHT HAND SIDE VECTOR.
988C    IP = PIVOT VECTOR OBTAINED FROM DEC.
989C  DO NOT USE IF DEC HAS SET IER .NE. 0.
990C  OUTPUT..
991C    B = SOLUTION VECTOR, X .
992C-----------------------------------------------------------------------
993      IF (N .EQ. 1) GO TO 50
994      NM1 = N - 1
995      DO 20 K = 1,NM1
996        KP1 = K + 1
997        M = IP(K)
998        T = B(M)
999        B(M) = B(K)
1000        B(K) = T
1001        DO 10 I = KP1,N
1002 10       B(I) = B(I) + A(I,K)*T
1003 20     CONTINUE
1004      DO 40 KB = 1,NM1
1005        KM1 = N - KB
1006        K = KM1 + 1
1007        B(K) = B(K)/A(K,K)
1008        T = -B(K)
1009        DO 30 I = 1,KM1
1010 30       B(I) = B(I) + A(I,K)*T
1011 40     CONTINUE
1012 50   B(1) = B(1)/A(1,1)
1013      RETURN
1014C----------------------- END OF SUBROUTINE SOL -------------------------
1015      END
1016
1017
1018 
1019 
1020      SUBROUTINE FUNC_CHEM(N, T, Y, P)
1021      INCLUDE 'KPP_ROOT_params.h'
1022      INCLUDE 'KPP_ROOT_global.h'
1023      INTEGER N
1024      KPP_REAL   T, Told
1025      KPP_REAL   Y(NVAR), P(NVAR)
1026      Told = TIME
1027      TIME = T
1028      CALL Update_SUN()
1029      CALL Update_RCONST()
1030      CALL Fun( Y,  FIX, RCONST, P )
1031      TIME = Told
1032      RETURN
1033      END
1034
1035 
1036      SUBROUTINE JAC_CHEM(N, T, Y, J)
1037      INCLUDE 'KPP_ROOT_params.h'
1038      INCLUDE 'KPP_ROOT_global.h'
1039      INTEGER N
1040      KPP_REAL   Told, T
1041      KPP_REAL   Y(NVAR), J(LU_NONZERO)
1042      Told = TIME
1043      TIME = T
1044      CALL Update_SUN()
1045      CALL Update_RCONST()
1046      CALL Jac_SP( Y,  FIX, RCONST, J )
1047      TIME = Told
1048      RETURN
1049      END                                                                                                                 
1050
1051
1052
Note: See TracBrowser for help on using the repository browser.