source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/kpp_seulex.f @ 4072

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

Merge of branch palm4u into trunk

File size: 40.5 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
11      INTEGER i
12
13      PARAMETER (KM=12,KM2=2+KM*(KM+3)/2,NRDENS=NVAR)
14      PARAMETER (LWORK=2*NVAR*NVAR+(KM+8)*NVAR+4*KM+20+KM2*NRDENS)
15      PARAMETER (LIWORK=2*NVAR+KM+20+NRDENS)
16
17      KPP_REAL WORK(LWORK)
18      INTEGER IWORK(LIWORK)
19      EXTERNAL FUNC_CHEM,JAC_CHEM
20
21      ITOL=1     ! --- VECTOR TOLERANCES
22      IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
23      MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
24      MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
25      IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
26      IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
27      IDFX=0     ! --- INTERNAL TIME DERIVATIVE
28
29      DO i=1,20
30        IWORK(i) = 0
31        WORK(i) = 0.D0
32      ENDDO
33
34      CALL ATMSEULEX(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT,
35     &                  STEPMIN,RTOL,ATOL,ITOL,
36     &                  JAC_CHEM,IJAC,MLJAC,MUJAC,
37     &                  FUNC_CHEM,IMAS,MLJAC,MUJAC,
38     &                  WORK,LWORK,IWORK,LIWORK,IDID)
39
40      IF (IDID.LT.0) THEN
41        print *,'ATMSEULEX: Unsucessfull exit at T=',
42     &          TIN,' (IDID=',IDID,')'
43      ENDIF
44
45      RETURN
46      END
47
48
49      SUBROUTINE ATMSEULEX(N,FCN,IFCN,X,Y,XEND,H,
50     &                  RelTol,AbsTol,ITOL,
51     &                  JAC ,IJAC,MLJAC,MUJAC,
52     &                  MAS,IMAS,MLMAS,MUMAS,
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 EXTRAPOLATION-ALGORITHM, BASED ON THE
58C     LINEARLY IMPLICIT EULER METHOD (WITH STEP SIZE CONTROL
59C     AND ORDER SELECTION).
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@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
65C              INCLUSION OF DENSE OUTPUT BY E. HAIRER AND A. OSTERMANN
66C
67C     THIS CODE IS PART OF THE BOOK:
68C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
69C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
70C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
71C         SPRINGER-VERLAG (1991)
72C
73C     VERSION OF SEPTEMBER 30, 1995
74C
75C     INPUT PARAMETERS
76C     ----------------
77C     N           DIMENSION OF THE SYSTEM
78C
79C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
80C                 VALUE OF F(X,Y):
81C                    SUBROUTINE FCN(N,X,Y,F)
82C                    KPP_REAL X,Y(N),F(N)
83C                    F(1)=...   ETC.
84C                 RPAR, IPAR (SEE BELOW)
85C
86C     IFCN        GIVES INFORMATION ON FCN:
87C                    IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS)
88C                    IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS)
89C
90C     X           INITIAL X-VALUE
91C
92C     Y(N)        INITIAL VALUES FOR Y
93C
94C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
95C
96C     H           INITIAL STEP SIZE GUESS;
97C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
98C                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
99C                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
100C                 ADAPTS ITS STEP SIZE (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     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
149C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
150C
151C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
152C                 MATRIX M.
153C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
154C                 MATRIX AND NEEDS NOT TO BE DEFINED;
155C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
156C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
157C                    SUBROUTINE MAS(N,AM,LMAS)
158C                    KPP_REAL AM(LMAS,N)
159C                    AM(1,1)= ....
160C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
161C                    AS FULL MATRIX LIKE
162C                         AM(I,J) = M(I,J)
163C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
164C                    DIAGONAL-WISE AS
165C                         AM(I-J+MUMAS+1,J) = M(I,J).
166C
167C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
168C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
169C                       MATRIX, MAS IS NEVER CALLED.
170C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
171C
172C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
173C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
174C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
175C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
176C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
177C                       THE MAIN DIAGONAL).
178C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
179C
180C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
181C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
182C                 NEED NOT BE DEFINED IF MLMAS=N.
183C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
184C
185C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
186C                 NUMERICAL SOLUTION DURING INTEGRATION.
187C                 IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
188C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
189C                 IT MUST HAVE THE FORM
190C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,RC,LRC,IC,LIC,N,
191C                                       RPAR,IPAR,IRTRN)
192C                    KPP_REAL X,Y(N),RC(LRC),IC(LIC)
193C                    ....
194C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
195C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
196C                    THE FIRST GRID-POINT).
197C                 "XOLD" IS THE PRECEEDING GRID-POINT.
198C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
199C                    IS SET <0, SEULEX RETURNS TO THE CALLING PROGRAM.
200C                 DO NOT CHANGE THE ENTRIES OF RC(LRC),IC(LIC)!
201C
202C          -----  CONTINUOUS OUTPUT (IF IOUT=2): -----
203C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
204C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
205C                 THE KPP_REAL FUNCTION
206C                        >>>   CONTEX(I,S,RC,LRC,IC,LIC)   <<<
207C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
208C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
209C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
210C
211C     IOUT        GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
212C                    IOUT=0: SUBROUTINE IS NEVER CALLED
213C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
214C                    IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT
215C
216C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
217C                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
218C                 "LWORK" MUST BE AT LEAST
219C                        N*(LJAC+LMAS+LE1+KM+8)+4*KM+20+KM2*NRDENS
220C                 WHERE
221C                    KM2=2+KM*(KM+3)/2  AND  NRDENS=IWORK(6) (SEE BELOW)
222C                 AND
223C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
224C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
225C                 AND
226C                    LMAS=0              IF IMAS=0
227C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
228C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
229C                 AND
230C                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
231C                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.).
232C                 AND
233C                    KM=12               IF IWORK(3)=0
234C                    KM=IWORK(3)         IF IWORK(3).GT.0
235C
236C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
237C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
238C                 STORAGE REQUIREMENT IS
239C                         LWORK = 2*N*N+(KM+8)*N+4*KM+13+KM2*NRDENS.
240C                 IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST
241C                    N*(LJAC+KM+8)+(N-M1)*(LMAS+LE1)+4*KM+20+KM2*NRDENS
242C                 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE1 THE
243C                 NUMBER N CAN BE REPLACED BY N-M1.
244C
245C     LWORK       DECLARED LENGTH OF ARRAY "WORK".
246C
247C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
248C                 "LIWORK" MUST BE AT LEAST  2*N+KM+20+NRDENS.
249C
250C     LIWORK      DECLARED LENGTH OF ARRAY "IWORK".
251C
252C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH
253C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
254C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES.
255C
256C ----------------------------------------------------------------------
257C
258C     SOPHISTICATED SETTING OF PARAMETERS
259C     -----------------------------------
260C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
261C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(13)
262C              AS WELL AS IWORK(1),..,IWORK(4) DIFFERENT FROM ZERO.
263C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
264C
265C    IWORK(1)  IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN
266C              MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY
267C              ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN.
268C              IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N)
269C              AND NOT FOR IMPLICIT SYSTEMS (IMAS=1).
270C
271C    IWORK(2)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
272C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000.
273C
274C    IWORK(3)  THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION
275C              TABLE. THE DEFAULT VALUE (FOR IWORK(3)=0) IS 12.
276C              IF IWORK(3).NE.0 THEN IWORK(3) SHOULD BE .GE.3.
277C
278C    IWORK(4)  SWITCH FOR THE STEP SIZE SEQUENCE
279C              IF IWORK(4).EQ.1 THEN 1,2,3,4,6,8,12,16,24,32,48,...
280C              IF IWORK(4).EQ.2 THEN 2,3,4,6,8,12,16,24,32,48,64,...
281C              IF IWORK(4).EQ.3 THEN 1,2,3,4,5,6,7,8,9,10,...
282C              IF IWORK(4).EQ.4 THEN 2,3,4,5,6,7,8,9,10,11,...
283C              THE DEFAULT VALUE (FOR IWORK(4)=0) IS IWORK(4)=2.
284C
285C    IWORK(5)  PARAMETER "LAMBDA" OF DENSE OUTPUT; POSSIBLE VALUES
286C              ARE 0 AND 1; DEFAULT IWORK(5)=0.
287C
288C    IWORK(6)  = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT
289C              IS REQUIRED
290C
291C    IWORK(21),...,IWORK(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH
292C              DENSE OUTPUT IS REQUIRED
293C
294C       IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT
295C            Y(I)' = Y(I+M2)   FOR  I=1,...,M1,
296C       WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME
297C       CAN BE ACHIEVED BY SETTING THE FOLLOWING TWO PARAMETERS. E.G.,
298C       FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE
299C       VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2.
300C       FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS:
301C       - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE
302C              JACOBIAN HAVE TO BE STORED
303C              IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL
304C                 DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J)
305C                FOR I=1,N-M1 AND J=1,N.
306C              ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM )
307C                 DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2)
308C                FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM.
309C       - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL
310C                0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM)
311C                     PARTIAL F(I+M1) / PARTIAL Y(J+K*M2),  I,J=1,M2
312C                    ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH
313C                    OF THESE MM+1 SUBMATRICES
314C       - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES
315C                NEED NOT BE DEFINED IF MLJAC=N-M1
316C       - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND
317C              NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
318C              IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF
319C              DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX.
320C              IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL
321C                 AM(I,J) = M(I+M1,J+M1)     FOR I=1,N-M1 AND J=1,N-M1.
322C              ELSE, THE MASS MATRIX IS BANDED
323C                 AM(I-J+MUMAS+1,J) = M(I+M1,J+M1)
324C       - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL
325C                0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX
326C       - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX
327C                NEED NOT BE DEFINED IF MLMAS=N-M1
328C
329C    IWORK(9)  THE VALUE OF M1.  DEFAULT M1=0.
330C
331C    IWORK(10) THE VALUE OF M2.  DEFAULT M2=M1.
332C
333C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
334C
335C    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
336C
337C    WORK(3)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
338C              INCREASE WORK(3), TO 0.01 SAY, WHEN JACOBIAN EVALUATIONS
339C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER.
340C              DEFAULT MIN(1.0D-4,RelTol(1))
341C
342C    WORK(4), WORK(5)   PARAMETERS FOR STEP SIZE SELECTION
343C              THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS
344C              CHOSEN SUBJECT TO THE RESTRICTION
345C                 FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN
346C              WHERE FACMIN=WORK(4)**(1/(J-1))
347C              DEFAULT VALUES: WORK(4)=0.1D0, WORK(5)=4.D0
348C
349C    WORK(6), WORK(7)   PARAMETERS FOR THE ORDER SELECTION
350C              ORDER IS DECREASED IF    W(K-1) <= W(K)*WORK(6)
351C              ORDER IS INCREASED IF    W(K) <= W(K-1)*WORK(7)
352C              DEFAULT VALUES: WORK(6)=0.7D0, WORK(7)=0.9D0
353C
354C    WORK(8), WORK(9)   SAFETY FACTORS FOR STEP CONTROL ALGORITHM
355C             HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1))
356C              DEFAULT VALUES: WORK(8)=0.8D0, WORK(9)=0.93D0
357C
358C    WORK(10), WORK(11), WORK(12), WORK(13)   ESTIMATED WORKS FOR
359C             A CALL TO  FCN, JAC, DEC, SOL, RESPECTIVELY.
360C             DEFAULT VALUES ARE: WORK(10)=1.D0, WORK(11)=5.D0,
361C             WORK(12)=1.D0, WORK(13)=1.D0.
362C
363C-----------------------------------------------------------------------
364C
365C     OUTPUT PARAMETERS
366C     -----------------
367C     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
368C                 (AFTER SUCCESSFUL RETURN X=XEND)
369C
370C     Y(N)        SOLUTION AT X
371C
372C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
373C
374C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
375C                   IDID=1  COMPUTATION SUCCESSFUL,
376C                   IDID=-1 COMPUTATION UNSUCCESSFUL.
377C
378C   IWORK(14)  NFCN    NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
379C                      EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
380C   IWORK(15)  NJAC    NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
381C                      OR NUMERICALLY)
382C   IWORK(16)  NSTEP   NUMBER OF COMPUTED STEPS
383C   IWORK(17)  NACCPT  NUMBER OF ACCEPTED STEPS
384C   IWORK(18)  NREJCT  NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
385C                      HAS BEEN ACCEPTED)
386C   IWORK(19)  NDEC    NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
387C   IWORK(20)  NSOL    NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
388C-----------------------------------------------------------------------
389C *** *** *** *** *** *** *** *** *** *** *** *** ***
390C          DECLARATIONS
391C *** *** *** *** *** *** *** *** *** *** *** *** ***
392      IMPLICIT KPP_REAL (A-H,O-Z)
393      DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK)
394      LOGICAL AUTNMS,IMPLCT,ARRET,JBAND
395      EXTERNAL FCN,JAC,MAS
396C *** *** *** *** *** *** ***
397C        SETTING THE PARAMETERS
398C *** *** *** *** *** *** ***
399      IOUT = 0
400      NFCN=0
401      NJAC=0
402      NSTEP=0
403      NACCPT=0
404      NREJCT=0
405      NDEC=0
406      NSOL=0
407      ARRET=.FALSE.
408C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
409      IF(IWORK(2).EQ.0)THEN
410         NMAX=100000
411      ELSE
412         NMAX=IWORK(2)
413         IF(NMAX.LE.0)THEN
414            WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
415            ARRET=.TRUE.
416         END IF
417      END IF
418C -------- KM     MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION
419      IF(IWORK(3).EQ.0)THEN
420         KM=12
421      ELSE
422         KM=IWORK(3)
423         IF(KM.LE.2)THEN
424            WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
425            ARRET=.TRUE.
426         END IF
427      END IF
428C -------- NSEQU     CHOICE OF STEP SIZE SEQUENCE
429      NSEQU=IWORK(4)
430      IF(IWORK(4).EQ.0) NSEQU=2
431      IF(NSEQU.LE.0.OR.NSEQU.GE.5)THEN
432         WRITE(6,*)' CURIOUS INPUT IWORK(4)=',IWORK(4)
433         ARRET=.TRUE.
434      END IF
435C -------- LAMBDA   PARAMETER FOR DENSE OUTPUT
436      LAMBDA=IWORK(5)
437      IF(LAMBDA.LT.0.OR.LAMBDA.GE.2)THEN
438         WRITE(6,*)' CURIOUS INPUT IWORK(5)=',IWORK(5)
439         ARRET=.TRUE.
440      END IF
441C -------- NRDENS   NUMBER OF DENSE OUTPUT COMPONENTS
442      NRDENS=IWORK(6)
443      IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN
444         WRITE(6,*)' CURIOUS INPUT IWORK(6)=',IWORK(6)
445         ARRET=.TRUE.
446      END IF
447C -------- PARAMETER FOR SECOND ORDER EQUATIONS
448      M1=IWORK(9)
449      M2=IWORK(10)
450      NM1=N-M1
451      IF (M1.EQ.0) M2=N
452      IF (M2.EQ.0) M2=M1
453      IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN
454       WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
455       ARRET=.TRUE.
456      END IF
457C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0
458      IF(WORK(1).EQ.0.D0)THEN
459         UROUND=1.D-16
460      ELSE
461         UROUND=WORK(1)
462         IF(UROUND.LE.0.D0.OR.UROUND.GE.1.D0)THEN
463            WRITE(6,*)'  UROUND=',WORK(1)
464            ARRET=.TRUE.
465         END IF
466      END IF
467C -------- MAXIMAL STEP SIZE
468      IF(WORK(2).EQ.0.D0)THEN
469         HMAX=XEND-X
470      ELSE
471         HMAX=WORK(2)
472      END IF
473C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
474      IF(WORK(3).EQ.0.D0)THEN
475         THET=MIN(1.0D-4,RelTol(1))
476      ELSE
477         THET=WORK(3)
478      END IF
479C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
480      IF(WORK(4).EQ.0.D0)THEN
481         FAC1=0.1D0
482      ELSE
483         FAC1=WORK(4)
484      END IF
485      IF(WORK(5).EQ.0.D0)THEN
486         FAC2=4.0D0
487      ELSE
488         FAC2=WORK(5)
489      END IF
490C -------  FAC3, FAC4   PARAMETERS FOR THE ORDER SELECTION
491      IF(WORK(6).EQ.0.D0)THEN
492         FAC3=0.7D0
493      ELSE
494         FAC3=WORK(6)
495      END IF
496      IF(WORK(7).EQ.0.D0)THEN
497         FAC4=0.9D0
498      ELSE
499         FAC4=WORK(7)
500      END IF
501C ------- SAFE1, SAFE2 SAFETY FACTORS FOR STEP SIZE PREDICTION
502      IF(WORK(8).EQ.0.D0)THEN
503         SAFE1=0.6D0
504      ELSE
505         SAFE1=WORK(8)
506      END IF
507      IF(WORK(9).EQ.0.D0)THEN
508         SAFE2=0.93D0
509      ELSE
510         SAFE2=WORK(9)
511      END IF
512C ------- WKFCN,WKJAC,WKDEC,WKSOL  ESTIMATED WORK FOR  FCN,JAC,DEC,SOL
513      IF(WORK(10).EQ.0.D0)THEN
514         WKFCN=1.D0
515      ELSE
516         WKFCN=WORK(10)
517      END IF
518      IF(WORK(11).EQ.0.D0)THEN
519         WKJAC=5.D0
520      ELSE
521         WKJAC=WORK(11)
522      END IF
523      IF(WORK(12).EQ.0.D0)THEN
524         WKDEC=1.D0
525      ELSE
526         WKDEC=WORK(12)
527      END IF
528      IF(WORK(13).EQ.0.D0)THEN
529         WKSOL=1.D0
530      ELSE
531         WKSOL=WORK(13)
532      END IF
533      WKROW=WKFCN+WKSOL
534C --------- CHECK IF TOLERANCES ARE O.K.
535      IF (ITOL.EQ.0) THEN
536          IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
537              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
538              ARRET=.TRUE.
539          END IF
540      ELSE
541        DO I=1,N
542          IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
543              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
544              ARRET=.TRUE.
545          END IF
546        END DO
547      END IF
548C *** *** *** *** *** *** *** *** *** *** *** *** ***
549C         COMPUTATION OF ARRAY ENTRIES
550C *** *** *** *** *** *** *** *** *** *** *** *** ***
551C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
552      AUTNMS=IFCN.EQ.0
553      IMPLCT=IMAS.NE.0
554      JBAND=MLJAC.LT.NM1
555C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
556C -- JACOBIAN AND MATRIX E
557      IF(JBAND)THEN
558         LDJAC=MLJAC+MUJAC+1
559         LDE=MLJAC+LDJAC
560      ELSE
561         MLJAC=NM1
562         MUJAC=NM1
563         LDJAC=NM1
564         LDE=NM1
565      END IF
566C -- MASS MATRIX
567      IF (IMPLCT) THEN
568          IF (MLMAS.NE.NM1) THEN
569              LDMAS=MLMAS+MUMAS+1
570              IF (JBAND) THEN
571                 IJOB=4
572              ELSE
573                 IJOB=3
574              END IF
575          ELSE
576              LDMAS=NM1
577              IJOB=5
578          END IF
579C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC"
580          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
581              WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF
582     & "JAC"'
583            ARRET=.TRUE.
584          END IF
585      ELSE
586          LDMAS=0
587          IF (JBAND) THEN
588             IJOB=2
589          ELSE
590             IJOB=1
591             IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7
592          END IF
593      END IF
594      LDMAS2=MAX(1,LDMAS)
595C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
596      IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN
597         WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
598     &FULL JACOBIAN'
599         ARRET=.TRUE.
600      END IF
601C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
602      KM2=(KM*(KM+1))/2
603      IEYH=21
604      IEDY=IEYH+N
605      IEFX=IEDY+N
606      IEYHH=IEFX+N
607      IEDYH=IEYHH+N
608      IEDEL=IEDYH+N
609      IEWH =IEDEL+N
610      IESCAL=IEWH+N
611      IEHH =IESCAL+N
612      IEW  =IEHH+KM
613      IEA  =IEW+KM
614      IEJAC =IEA+KM
615      IEE  =IEJAC+N*LDJAC
616      IEMAS=IEE+NM1*LDE
617      IET=IEMAS+NM1*LDMAS
618      IFAC=IET+N*KM
619      IEDE=IFAC+KM
620      IFSAFE=IEDE+(KM+2)*NRDENS
621C ------ TOTAL STORAGE REQUIREMENT -----------
622      ISTORE=IFSAFE+KM2*NRDENS-1
623      IF(ISTORE.GT.LWORK)THEN
624         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
625         ARRET=.TRUE.
626      END IF
627C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
628      IECO=21
629      IEIP=21+NRDENS
630      IENJ=IEIP+N
631      IEIPH=IENJ+KM
632C --------- TOTAL REQUIREMENT ---------------
633      ISTORE=IECO+KM-1
634      IF(ISTORE.GT.LIWORK)THEN
635         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
636         ARRET=.TRUE.
637      END IF
638C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
639      IF (ARRET) THEN
640         IDID=-1
641         RETURN
642      END IF
643      NRD=MAX(1,NRDENS)
644C -------- CALL TO CORE INTEGRATOR ------------
645      CALL SEUCOR(N,FCN,X,Y,XEND,HMAX,H,KM,RelTol,AbsTol,ITOL,JAC,IJAC,
646     &   MLJAC,MUJAC,MLMAS,MUMAS,IOUT,IDID,IJOB,M1,M2,NM1,
647     &   NMAX,UROUND,NSEQU,AUTNMS,IMPLCT,JBAND,LDJAC,LDE,LDMAS2,
648     &   WORK(IEYH),WORK(IEDY),WORK(IEFX),WORK(IEYHH),WORK(IEDYH),
649     &   WORK(IEDEL),WORK(IEWH),WORK(IESCAL),WORK(IEHH),
650     &   WORK(IEW),WORK(IEA),WORK(IEJAC),WORK(IEE),WORK(IEMAS),
651     &   WORK(IET),IWORK(IEIP),IWORK(IENJ),IWORK(IEIPH),FAC1,FAC2,FAC3,
652     &   FAC4,THET,SAFE1,SAFE2,WKJAC,WKDEC,WKROW,KM2,NRD,WORK(IFAC),
653     &   WORK(IFSAFE),LAMBDA,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,
654     &   WORK(IEDE),IWORK(IECO))
655      IWORK(14)=NFCN
656      IWORK(15)=NJAC
657      IWORK(16)=NSTEP
658      IWORK(17)=NACCPT
659      IWORK(18)=NREJCT
660      IWORK(19)=NDEC
661      IWORK(20)=NSOL
662C ----------- RETURN -----------
663      RETURN
664      END
665C
666C
667C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
668C
669      SUBROUTINE SEUCOR(N,FCN,X,Y,XEND,HMAX,H,KM,RelTol,AbsTol,ITOL,
670     &  JAC,IJAC,MLJAC,MUJAC,MLB,MUB,IOUT,IDID,IJOB,M1,M2,
671     &  NM1,NMAX,UROUND,NSEQU,AUTNMS,IMPLCT,BANDED,LFJAC,LE,
672     &  LDMAS,YH,DY,FX,YHH,DYH,DEL,WH,SCAL,HH,W,A,FJAC,E,FMAS,T,IP,
673     &  NJ,IPHES,FAC1,FAC2,FAC3,FAC4,THET,SAFE1,SAFE2,WKJAC,WKDEC,WKROW,
674     &  KM2,NRD,FACUL,FSAFE,LAMBDA,NFCN,NJAC,NSTEP,NACCPT,NREJCT,
675     &  NDEC,NSOL,DENS,ICOMP)
676C ----------------------------------------------------------
677C     CORE INTEGRATOR FOR SEULEX
678C     PARAMETERS SAME AS IN SEULEX WITH WORKSPACE ADDED
679C ----------------------------------------------------------
680C         DECLARATIONS
681C ----------------------------------------------------------
682       IMPLICIT KPP_REAL (A-H,O-Z)
683       INCLUDE 'KPP_ROOT_Parameters.h'
684       INCLUDE 'KPP_ROOT_Sparse.h'
685       INTEGER KM, KM2, K, KC, KRIGHT, KLR, KK, KRN, KOPT
686       DIMENSION Y(N),YH(N),DY(N),FX(N),YHH(N),DYH(N),DEL(N)
687       DIMENSION WH(N),SCAL(N),HH(KM),W(KM),A(KM),FJAC(LU_NONZERO)
688       DIMENSION FMAS(LDMAS,NM1),T(KM,N),IP(NM1),NJ(KM)
689       DIMENSION RelTol(*),AbsTol(*)
690       DIMENSION IPHES(N),FSAFE(KM2,NRD),FACUL(KM),E(LE,NM1)
691       DIMENSION DENS((KM+2)*NRD),ICOMP(NRD)
692       LOGICAL REJECT,LAST,ATOV,CALJAC,CALHES,AUTNMS,IMPLCT,BANDED
693       EXTERNAL FCN, JAC
694       COMMON /COSEU/XOLDD,HHH,NNRD,KRIGHT
695       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
696C --- COMPUTE COEFFICIENTS FOR DENSE OUTPUT
697       IF (IOUT.EQ.2) THEN
698          NNRD=NRD
699C --- COMPUTE THE FACTORIALS --------
700          FACUL(1)=1.D0
701          DO I=1,KM-1
702             FACUL(I+1)=I*FACUL(I)
703          END DO
704       END IF
705
706C *** *** *** *** *** *** ***
707C  INITIALISATIONS
708C *** *** *** *** *** *** ***
709       LRDE=(KM+2)*NRD
710       MLE=MLJAC
711       MUE=MUJAC
712       MBJAC=MLJAC+MUJAC+1
713       MBB=MLB+MUB+1
714       MDIAG=MLE+MUE+1
715       MDIFF=MLE+MUE-MUB
716       MBDIAG=MUB+1
717       IF (M1.GT.0) IJOB=IJOB+10
718C --- DEFINE THE STEP SIZE SEQUENCE
719       IF (NSEQU.EQ.1) THEN
720           NJ(1)=1
721           NJ(2)=2
722           NJ(3)=3
723           DO I=4,KM
724              NJ(I)=2*NJ(I-2)
725           END DO
726       END IF
727       IF (NSEQU.EQ.2) THEN
728           NJ(1)=2
729           NJ(2)=3
730           DO I=3,KM
731              NJ(I)=2*NJ(I-2)
732           END DO
733       END IF
734       DO I=1,KM
735          IF (NSEQU.EQ.3) NJ(I)=I
736          IF (NSEQU.EQ.4) NJ(I)=I+1
737       END DO
738       A(1)=WKJAC+NJ(1)*WKROW+WKDEC
739       DO I=2,KM
740          A(I)=A(I-1)+(NJ(I)-1)*WKROW+WKDEC
741       END DO
742       K=MAX0(3,MIN0(KM-2,INT(-DLOG10(RelTol(1)+AbsTol(1))*.6D0+1.5D0)))
743       HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
744       IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
745       H=MIN(ABS(H),HMAXN)
746       THETA=2*ABS(THET)
747       ERR=0.D0
748       W(1)=1.D30
749       DO I=1,N
750          IF (ITOL.EQ.0) THEN
751            SCAL(I)=AbsTol(1)+RelTol(1)*DABS(Y(I))
752          ELSE
753            SCAL(I)=AbsTol(I)+RelTol(I)*DABS(Y(I))
754          END IF
755       END DO
756       CALJAC=.FALSE.
757       REJECT=.FALSE.
758       LAST=.FALSE.
759  10   CONTINUE
760       IF (REJECT) THETA=2*ABS(THET)
761       ATOV=.FALSE.
762C *** *** *** *** *** *** ***
763C --- IS XEND REACHED IN THE NEXT STEP?
764C *** *** *** *** *** *** ***
765       H1=XEND-X
766       IF (H1.LE.UROUND) GO TO 110
767       HOPT=H
768       H=MIN(H,H1,HMAXN)
769       IF (H.GE.H1-UROUND) LAST=.TRUE.
770       IF (AUTNMS) THEN
771          CALL FCN(N,X,Y,DY)
772          NFCN=NFCN+1
773       END IF
774       IF (THETA.GT.THET.AND..NOT.CALJAC) THEN
775C *** *** *** *** *** *** ***
776C  COMPUTATION OF THE JACOBIAN
777C *** *** *** *** *** *** ***
778          NJAC=NJAC+1
779C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
780          CALL JAC(N,X,Y,FJAC)
781          CALJAC=.TRUE.
782          CALHES=.FALSE.
783       END IF
784C *** *** *** *** *** *** ***
785C --- THE FIRST AND LAST STEP
786C *** *** *** *** *** *** ***
787      IF (NSTEP.EQ.0.OR.LAST) THEN
788        IPT=0
789        NSTEP=NSTEP+1
790        DO J=1,K
791         KC=J
792         CALL SEUL(J,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
793     &             HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,FAC,
794     &             FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,MLB,
795     &             MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
796     &             ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES)
797         IF (ATOV) GOTO 10
798         IF (J.GT.1.AND.ERR.LE.1.D0) GO TO 60
799        END DO
800        GO TO 55
801      END IF
802C --- BASIC INTEGRATION STEP
803  30  CONTINUE
804      IPT=0
805      NSTEP=NSTEP+1
806      IF (NSTEP.GE.NMAX) GO TO 120
807      KC=K-1
808      DO J=1,KC
809       CALL SEUL(J,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
810     &           HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,
811     &           FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
812     &           MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
813     &           ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES)
814       IF (ATOV) GOTO 10
815      END DO
816C *** *** *** *** *** *** ***
817C --- CONVERGENCE MONITOR
818C *** *** *** *** *** *** ***
819      IF (K.EQ.2.OR.REJECT) GO TO 50
820      IF (ERR.LE.1.D0) GO TO 60
821      IF (ERR.GT.DBLE(NJ(K+1)*NJ(K))*4.D0) GO TO 100
822  50  CALL SEUL(K,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
823     &           HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,
824     &           FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
825     &           MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
826     &           ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES)
827      IF (ATOV) GOTO 10
828      KC=K
829      IF (ERR.LE.1.D0) GO TO 60
830C --- HOPE FOR CONVERGENCE IN LINE K+1
831  55  IF (ERR.GT.DBLE(NJ(K+1))*2.D0) GO TO 100
832      KC=K+1
833      CALL SEUL(KC,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
834     &           HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,
835     &           FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
836     &           MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
837     &           ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES)
838      IF (ATOV) GOTO 10
839CAdi      IF (ERR.GT.1.D0) GO TO 100
840      IF ((ERR.GT.1.D0).and.(H.gt.STEPMIN)) GO TO 100 ! Adi
841C *** *** *** *** *** *** ***
842C --- STEP IS ACCEPTED
843C *** *** *** *** *** *** ***
844  60  XOLD=X
845      X=X+H
846      IF (IOUT.EQ.2) THEN
847         KRIGHT=KC
848         DO I=1,NRD
849            DENS(I)=Y(ICOMP(I))
850         END DO
851      END IF
852      DO I=1,N
853         T1I=T(1,I)
854         IF (ITOL.EQ.0) THEN
855            SCAL(I)=AbsTol(1)+RelTol(1)*DABS(T1I)
856         ELSE
857            SCAL(I)=AbsTol(I)+RelTol(I)*DABS(T1I)
858         END IF
859         Y(I)=T1I
860      END DO
861      NACCPT=NACCPT+1
862      CALJAC=.FALSE.
863      IF (IOUT.EQ.2) THEN
864         XOLDD=XOLD
865         HHH=H
866         DO I=1,NRD
867            DENS(NRD+I)=Y(ICOMP(I))
868         END DO
869         DO KLR=1,KRIGHT-1
870C --- COMPUTE DIFFERENCES
871            IF (KLR.GE.2) THEN
872               DO KK=KLR,KC
873                  LBEG=((KK+1)*KK)/2
874                  LEND=LBEG-KK+2
875                  DO L=LBEG,LEND,-1
876                     DO I=1,NRD
877                        FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I)
878                     END DO
879                  END DO
880               END DO
881             END IF
882C --- COMPUTE DERIVATIVES AT RIGHT END ----
883             DO KK=KLR+LAMBDA,KC
884                FACNJ=NJ(KK)
885                FACNJ=FACNJ**KLR/FACUL(KLR+1)
886                IPT=((KK+1)*KK)/2
887                DO  I=1,NRD
888                   KRN=(KK-LAMBDA+1)*NRD
889                   DENS(KRN+I)=FSAFE(IPT,I)*FACNJ
890                END DO
891             END DO
892             DO J=KLR+LAMBDA+1,KC
893               DBLENJ=NJ(J)
894               DO L=J,KLR+LAMBDA+1,-1
895                 FACTOR=DBLENJ/NJ(L-1)-1.D0
896                 DO I=1,NRD
897                   KRN=(L-LAMBDA+1)*NRD+I
898                   DENS(KRN-NRD)=DENS(KRN)
899     &                          +(DENS(KRN)-DENS(KRN-NRD))/FACTOR
900                 END DO
901               END DO
902             END DO
903         END DO
904C ---  COMPUTE THE COEFFICIENTS OF THE INTERPOLATION POLYNOMIAL
905         DO IN=1,NRD
906            DO J=1,KRIGHT
907               II=NRD*J+IN
908               DENS(II)=DENS(II)-DENS(II-NRD)
909            END DO
910         END DO
911      END IF
912C --- COMPUTE OPTIMAL ORDER
913      IF (KC.EQ.2) THEN
914         KOPT=3
915         IF (REJECT) KOPT=2
916         GO TO 80
917      END IF
918      IF (KC.LE.K) THEN
919         KOPT=KC
920         IF (W(KC-1).LT.W(KC)*FAC3) KOPT=KC-1
921         IF (W(KC).LT.W(KC-1)*FAC4) KOPT=MIN0(KC+1,KM-1)
922      ELSE
923         KOPT=KC-1
924         IF (KC.GT.3.AND.W(KC-2).LT.W(KC-1)*FAC3) KOPT=KC-2
925         IF (W(KC).LT.W(KOPT)*FAC4) KOPT=MIN0(KC,KM-1)
926      END IF
927C --- AFTER A REJECTED STEP
928  80  IF (REJECT) THEN
929         K=MIN0(KOPT,KC)
930         H=DMIN1(H,HH(K))
931         REJECT=.FALSE.
932         GO TO 10
933      END IF
934C --- COMPUTE STEP SIZE FOR NEXT STEP
935      IF (KOPT.LE.KC) THEN
936         H=HH(KOPT)
937      ELSE
938         IF (KC.LT.K.AND.W(KC).LT.W(KC-1)*FAC4) THEN
939            H=HH(KC)*A(KOPT+1)/A(KC)
940         ELSE
941            H=HH(KC)*A(KOPT)/A(KC)
942         END IF
943      END IF
944      K=KOPT
945      H = DMAX1(H, STEPMIN)  ! Adi
946      GO TO 10
947C *** *** *** *** *** *** ***
948C --- STEP IS REJECTED
949C *** *** *** *** *** *** ***
950 100  K=MIN0(K,KC)
951      IF (K.GT.2.AND.W(K-1).LT.W(K)*FAC3) K=K-1
952      NREJCT=NREJCT+1
953      H=HH(K)
954      LAST=.FALSE.
955      REJECT=.TRUE.
956      IF (CALJAC) GOTO 30
957      GO TO 10
958C --- SOLUTION EXIT
959 110  CONTINUE
960      H=HOPT
961      IDID=1
962      RETURN
963C --- FAIL EXIT
964 120  WRITE (6,979) X,H
965 979  FORMAT(' EXIT OF SEULEX AT X=',D14.7,'   H=',D14.7)
966      IDID=-1
967      RETURN
968      END
969C
970C
971C *** *** *** *** *** *** ***
972C     S U B R O U T I N E    S E U L
973C *** *** *** *** *** *** ***
974C
975      SUBROUTINE SEUL(JJ,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,
976     &          H,KM,HMAXN,T,SCAL,NJ,HH,W,A,YH,DYH,DEL,WH,ERR,SAFE1,
977     &          FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
978     &          MLB,MUB,ERROLD,IPHES,ICOMP,
979     &          AUTNMS,IMPLCT,BANDED,REJECT,ATOV,FSAFE,KM2,NRD,IOUT,
980     &          IPT,M1,M2,NM1,IJOB,CALHES)
981C --- THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE
982C --- EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATE
983C --- OF THE OPTIMAL STEP SIZE
984      IMPLICIT KPP_REAL (A-H,O-Z)
985      INCLUDE 'KPP_ROOT_Parameters.h'
986      INCLUDE 'KPP_ROOT_Sparse.h'
987      INTEGER KM, KM2
988      DIMENSION Y(N),YH(N),DY(N),FX(N),DYH(N),DEL(N)
989      DIMENSION WH(N),SCAL(N),HH(KM),W(KM),A(KM)
990      DIMENSION FJAC(LU_NONZERO),E(LU_NONZERO)
991      DIMENSION FMAS(LDMAS,N),T(KM,N),IP(N),NJ(KM),IPHES(N)
992      DIMENSION FSAFE(KM2,NRD),ICOMP(NRD)
993      LOGICAL ATOV,REJECT,AUTNMS,IMPLCT,BANDED,CALHES
994      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
995      EXTERNAL FCN
996C *** *** *** *** *** *** ***
997C  COMPUTE THE MATRIX E AND ITS DECOMPOSITION
998C *** *** *** *** *** *** ***
999      HJ=H/NJ(JJ)
1000      HJI=1.D0/HJ
1001      DO  I=1,LU_NONZERO
1002         E(I)=-FJAC(I)
1003      END DO
1004      DO J=1,N
1005         E(LU_DIAG(J))=E(LU_DIAG(J))+HJI
1006      END DO
1007      CALL KppDecomp (E,IER)
1008      IF (IER.NE.0) GOTO 79
1009      NDEC=NDEC+1
1010C *** *** *** *** *** *** ***
1011C --- STARTING PROCEDURE
1012C *** *** *** *** *** *** ***
1013       IF (.NOT.AUTNMS) THEN
1014          CALL FCN(N,X+HJ,Y,DY)
1015          NFCN=NFCN+1
1016       END IF
1017       DO I=1,N
1018          YH(I)=Y(I)
1019          DEL(I)=DY(I)
1020       END DO
1021       CALL KppSolve (E,DEL)
1022       NSOL=NSOL+1
1023       M=NJ(JJ)
1024       IF (IOUT.EQ.2.AND.M.EQ.JJ) THEN
1025          IPT=IPT+1
1026          DO I=1,NRD
1027             FSAFE(IPT,I)=DEL(ICOMP(I))
1028          END DO
1029       END IF
1030C *** *** *** *** *** *** ***
1031C --- SEMI-IMPLICIT EULER METHOD
1032C *** *** *** *** *** *** ***
1033       IF (M.GT.1) THEN
1034          DO MM=1,M-1
1035             DO I=1,N
1036                YH(I)=YH(I)+DEL(I)
1037             END DO
1038             IF (AUTNMS) THEN
1039                CALL FCN(N,X+HJ*MM,YH,DYH)
1040             ELSE
1041                CALL FCN(N,X+HJ*(MM+1),YH,DYH)
1042             END IF
1043             NFCN=NFCN+1
1044             IF (MM.EQ.1.AND.JJ.LE.2) THEN
1045C --- STABILITY CHECK
1046                DEL1=0.D0
1047                DO I=1,N
1048                   DEL1=DEL1+(DEL(I)/SCAL(I))**2
1049                END DO
1050                DEL1=DSQRT(DEL1)
1051                IF (IMPLCT) THEN
1052                   DO I=1,NM1
1053                      WH(I)=DEL(I+M1)
1054                   END DO
1055                   IF (MLB.EQ.NM1) THEN
1056                      DO I=1,NM1
1057                         SUM=0.D0
1058                         DO J=1,NM1
1059                            SUM=SUM+FMAS(I,J)*WH(J)
1060                         END DO
1061                         DEL(I+M1)=SUM
1062                      END DO
1063                   ELSE
1064                      DO I=1,NM1
1065                         SUM=0.D0
1066                         DO J=MAX(1,I-MLB),MIN(NM1,I+MUB)
1067                            SUM=SUM+FMAS(I-J+MBDIAG,J)*WH(J)
1068                         END DO
1069                         DEL(I+M1)=SUM
1070                      END DO
1071                   END IF
1072                END IF
1073                IF (.NOT.AUTNMS) THEN
1074                   CALL FCN(N,X+HJ,YH,WH)
1075                   NFCN=NFCN+1
1076                   DO I=1,N
1077                      DEL(I)=WH(I)-DEL(I)*HJI
1078                   END DO
1079                ELSE
1080                   DO I=1,N
1081                      DEL(I)=DYH(I)-DEL(I)*HJI
1082                   END DO
1083                END IF
1084                CALL KppSolve (E,DEL)
1085                NSOL=NSOL+1
1086                DEL2=0.D0
1087                DO I=1,N
1088                   DEL2=DEL2+(DEL(I)/SCAL(I))**2
1089                END DO
1090                DEL2=DSQRT(DEL2)
1091                THETA=DEL2/MAX(1.D0,DEL1)
1092                IF (THETA.GT.1.D0) GOTO 79
1093             END IF
1094             CALL KppSolve (E,DYH)
1095             NSOL=NSOL+1
1096             DO I=1,N
1097                DEL(I)=DYH(I)
1098             END DO
1099             IF (IOUT.EQ.2.AND.MM.GE.M-JJ) THEN
1100                IPT=IPT+1
1101                DO I=1,NRD
1102                   FSAFE(IPT,I)=DEL(ICOMP(I))
1103                END DO
1104             END IF
1105          END DO
1106       END IF
1107       DO I=1,N
1108          T(JJ,I)=YH(I)+DEL(I)
1109       END DO
1110C *** *** *** *** *** *** ***
1111C --- POLYNOMIAL EXTRAPOLATION
1112C *** *** *** *** *** *** ***
1113       IF (JJ.EQ.1) RETURN
1114       DO L=JJ,2,-1
1115          FAC=(DBLE(NJ(JJ))/DBLE(NJ(L-1)))-1.D0
1116          DO I=1,N
1117             T(L-1,I)=T(L,I)+(T(L,I)-T(L-1,I))/FAC
1118          END DO
1119       END DO
1120       ERR=0.D0
1121       DO I=1,N
1122          ERR=ERR+MIN(ABS((T(1,I)-T(2,I)))/SCAL(I),1.D15)**2
1123       END DO
1124       IF (ERR.GE.1.D30) GOTO 79
1125       ERR=DSQRT(ERR/DBLE(N))
1126       IF (JJ.GT.2.AND.ERR.GE.ERROLD) GOTO 79
1127       ERROLD=DMAX1(4*ERR,1.D0)
1128C --- COMPUTE OPTIMAL STEP SIZES
1129       EXPO=1.D0/JJ
1130       FACMIN=FAC1**EXPO
1131       FAC=DMIN1(FAC2/FACMIN,DMAX1(FACMIN,(ERR/SAFE1)**EXPO/SAFE2))
1132       FAC=1.D0/FAC
1133       HH(JJ)=DMIN1(H*FAC,HMAXN)
1134       W(JJ)=A(JJ)/HH(JJ)
1135       RETURN
1136  79   ATOV=.TRUE.
1137       H=H*0.5D0
1138       REJECT=.TRUE.
1139       RETURN
1140       END
1141
1142
1143
1144      SUBROUTINE FUNC_CHEM(N, T, Y, P)
1145      INCLUDE 'KPP_ROOT_Parameters.h'
1146      INCLUDE 'KPP_ROOT_Global.h'
1147      INTEGER N
1148      KPP_REAL   T, Told
1149      KPP_REAL   Y(NVAR), P(NVAR)
1150      Told = TIME
1151      TIME = T
1152      CALL Update_SUN()
1153      CALL Update_RCONST()
1154      CALL Fun( Y,  FIX, RCONST, P )
1155      TIME = Told
1156      RETURN
1157      END
1158
1159
1160      SUBROUTINE JAC_CHEM(N, T, Y, J)
1161      INCLUDE 'KPP_ROOT_Parameters.h'
1162      INCLUDE 'KPP_ROOT_Global.h'
1163      INTEGER N
1164      KPP_REAL   Told, T
1165      KPP_REAL   Y(NVAR), J(LU_NONZERO)
1166      Told = TIME
1167      TIME = T
1168      CALL Update_SUN()
1169      CALL Update_RCONST()
1170      CALL Jac_SP( Y,  FIX, RCONST, J )
1171      TIME = Told
1172      RETURN
1173      END
1174
Note: See TracBrowser for help on using the repository browser.