source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/sdirk.f @ 3821

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

Merge of branch palm4u into trunk

File size: 19.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 (LWORK=2*NVAR*NVAR+12*NVAR+7,LIWORK=2*NVAR+7)
14      PARAMETER (LRCONT=5*NVAR+2)
15
16      KPP_REAL WORK(LWORK)
17      INTEGER IWORK(LIWORK)   
18      COMMON /CONT/ ICONT(4),RCONT(LRCONT)
19      EXTERNAL FUNC_CHEM,JAC_CHEM
20
21      ITOL=1     ! --- VECTOR TOLERANCES
22      IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
23      IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
24
25      DO i=1,20
26        IWORK(i) = 0
27        WORK(i) = 0.D0
28      ENDDO
29     
30      IWORK(3) = 8
31
32      CALL ATMSDIRK(NVAR,FUNC_CHEM,TIN,VAR,TOUT,STEPMIN,     
33     &                  RTOL,ATOL,ITOL,     
34     &                  JAC_CHEM ,IJAC, FUNC_CHEM ,IMAS,                   
35     &                  WORK,LWORK,IWORK,LIWORK,LRCONT,IDID)       
36
37      IF (IDID.LT.0) THEN
38        print *,'ATMSDIRK: Unsucessfull exit at T=',
39     &          TIN,' (IDID=',IDID,')'
40      ENDIF
41
42      RETURN
43      END
44
45
46      SUBROUTINE ATMSDIRK(N,FCN,X,Y,XEND,H,
47     &                  RelTol,AbsTol,ITOL,
48     &                  JAC ,IJAC, MAS ,IMAS,
49     &                  WORK,LWORK,IWORK,LIWORK,LRCONT,IDID)
50C ----------------------------------------------------------
51C *** *** *** *** *** *** *** *** *** *** *** *** ***
52C          DECLARATIONS 
53C *** *** *** *** *** *** *** *** *** *** *** *** *** 
54      IMPLICIT KPP_REAL (A-H,O-Z)
55      DIMENSION Y(N),AbsTol(1),RelTol(1),WORK(LWORK),IWORK(LIWORK)
56      LOGICAL IMPLCT,JBAND,ARRET
57      EXTERNAL FCN,JAC,MAS
58      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
59
60C *** *** *** *** *** *** ***
61C        SETTING THE PARAMETERS 
62C *** *** *** *** *** *** ***
63       NFCN=0
64       NJAC=0
65       NSTEP=0
66       NACCPT=0
67       NREJCT=0
68       NDEC=0
69       NSOL=0
70       ARRET=.FALSE.
71C -------- SWITCH FOR TRANSFORMATION OF JACOBIAN TO HESS_CHEM FORM ---
72      NHESS1 = 0         ! ADRIAN
73C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
74      IF(IWORK(2).EQ.0)THEN
75         NMAX=100000
76      ELSE
77         NMAX=IWORK(2)
78         IF(NMAX.LE.0)THEN
79            WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
80            ARRET=.TRUE.
81         END IF
82      END IF
83C -------- NIT    MAXIMAL NUMBER OF NEWTON ITERATIONS
84      IF(IWORK(3).EQ.0)THEN
85         NIT=8
86      ELSE
87         NIT=IWORK(3)
88         IF(NIT.LE.0)THEN
89            WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
90            ARRET=.TRUE.
91         END IF
92      END IF
93C -------- METH    SWITCH FOR THE COEFFICIENTS OF THE METHOD 
94      METH = 2
95C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 
96      IF(WORK(1).EQ.0.D0)THEN
97         UROUND=1.D-16
98      ELSE
99         UROUND=WORK(1)
100         IF(UROUND.LE.1.D-19.OR.UROUND.GE.1.D0)THEN
101            WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
102            ARRET=.TRUE.
103         END IF
104      END IF
105C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
106      IF(WORK(2).EQ.0.D0)THEN
107         SAFE=0.9D0
108      ELSE
109         SAFE=WORK(2)
110         IF(SAFE.LE..001D0.OR.SAFE.GE.1.D0)THEN
111            WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
112            ARRET=.TRUE.
113         END IF
114      END IF
115C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
116      IF(WORK(3).EQ.0.D0)THEN
117         THET=0.001D0
118      ELSE
119         THET=WORK(3)
120      END IF
121C --- FNEWT   STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
122      IF(WORK(4).EQ.0.D0)THEN
123         FNEWT=0.03D0
124      ELSE
125         FNEWT=WORK(4)
126      END IF
127C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
128      IF(WORK(5).EQ.0.D0)THEN
129         QUOT1=1.D0
130      ELSE
131         QUOT1=WORK(5)
132      END IF
133      IF(WORK(6).EQ.0.D0)THEN
134         QUOT2=1.2D0
135      ELSE
136         QUOT2=WORK(6)
137      END IF
138C -------- MAXIMAL STEP SIZE
139      IF(WORK(7).EQ.0.D0)THEN
140         HMAX=XEND-X
141      ELSE
142         HMAX=WORK(7)
143      END IF
144C --------- CHECK IF TOLERANCES ARE O.K.
145      IF (ITOL.EQ.0) THEN
146          IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
147              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
148              ARRET=.TRUE.
149          END IF
150      ELSE
151          DO 15 I=1,N
152          IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
153              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
154              ARRET=.TRUE.
155          END IF
156  15      CONTINUE
157      END IF
158
159C *** *** *** *** *** *** *** *** *** *** *** *** ***
160C         COMPUTATION OF ARRAY ENTRIES
161C *** *** *** *** *** *** *** *** *** *** *** *** ***
162C ---- IMPLICIT, BANDED OR NOT ?
163      IMPLCT=IMAS.NE.0
164      ARRET=.FALSE.
165C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
166C -- JACOBIAN
167         LDJAC=N
168C -- MATRIX E FOR LINEAR ALGEBRA
169         LDE=N
170C -- MASS MATRIX
171      IF (IMPLCT) THEN
172          print *,'IMPLCT 1'
173      ELSE
174          LDMAS=0
175      END IF
176      LDMAS2=MAX(1,LDMAS)
177
178C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
179      IEYHAT=8
180      IEZ=IEYHAT+N
181      IEY0=IEZ+N
182      IEZ1=IEY0+N
183      IEZ2=IEZ1+N
184      IEZ3=IEZ2+N
185      IEZ4=IEZ3+N
186      IEZ5=IEZ4+N
187      IESCAL=IEZ5+N
188      IEF1=IESCAL+N
189      IEG1=IEF1+N
190      IEH1=IEG1+N
191      IEJAC=IEH1+N
192      IEMAS=IEJAC+N*LDJAC
193      IEE=IEMAS+N*LDMAS
194
195C ------ TOTAL STORAGE REQUIREMENT -----------
196      ISTORE=IEE+N*LDE-1
197      IF(ISTORE.GT.LWORK)THEN
198         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
199         ARRET=.TRUE.
200      END IF
201C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
202      IEIP=5
203      IEHES=IEIP+N
204C --------- TOTAL REQUIREMENT ---------------
205      ISTORE=IEHES+N-1
206      IF(ISTORE.GT.LIWORK)THEN
207         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
208         ARRET=.TRUE.
209      END IF
210C --------- CONTROL OF LENGTH OF COMMON BLOCK "CONT" -------
211      IF(LRCONT.LT.(5*N+2))THEN
212         WRITE(6,*)' INSUFF. STORAGE FOR RCONT, MIN. LRCONT=',5*N+2
213         ARRET=.TRUE.
214      END IF
215C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
216      IF (ARRET) THEN
217         IDID=-1
218         RETURN
219      END IF
220C -------- CALL TO CORE INTEGRATOR ------------
221      CALL SDICOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
222     &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,IOUT,IDID,
223     &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,METH,NHESS1,
224     &   IMPLCT,JBAND,LDJAC,LDE,LDMAS2,
225     &   WORK(IEYHAT),WORK(IEZ),WORK(IEY0),WORK(IEZ1),WORK(IEZ2),
226     &   WORK(IEZ3),WORK(IEZ4),WORK(IEZ5),WORK(IESCAL),WORK(IEF1),
227     &   WORK(IEG1),WORK(IEH1),WORK(IEJAC),WORK(IEE),
228     &   WORK(IEMAS),IWORK(IEIP),IWORK(IEHES))
229C ----------- RETURN -----------
230      RETURN
231      END
232C
233C
234C
235C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
236C
237      SUBROUTINE SDICOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
238     &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,IOUT,IDID,
239     &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,METH,NHESS1,
240     &   IMPLCT,BANDED,LDJAC,LE,LDMAS,
241     &   YHAT,Z,Y0,Z1,Z2,Z3,Z4,Z5,SCAL,F1,G1,H1,FJAC,E,FMAS,IP,IPHES)
242C ----------------------------------------------------------
243C     CORE INTEGRATOR FOR SDIRK4
244C     PARAMETERS SAME AS IN SDIRK4 WITH WORKSPACE ADDED
245C ----------------------------------------------------------
246C         DECLARATIONS
247C ----------------------------------------------------------
248      IMPLICIT KPP_REAL (A-H,O-Z)
249      INCLUDE 'KPP_ROOT_Parameters.h'
250      INCLUDE 'KPP_ROOT_Sparse.h'
251      KPP_REAL Y(N),YHAT(N),Z(N),Y0(N),Z1(N),Z2(N),Z3(N),Z4(N),Z5(N)
252      KPP_REAL SCAL(N),F1(N),G1(N),H1(N)
253      KPP_REAL FJAC(LU_NONZERO),E(LU_NONZERO),FMAS(LDMAS,N)
254      KPP_REAL AbsTol(1),RelTol(1)
255      INTEGER IP(N),IPHES(N)
256      LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,NEWTRE
257      COMMON /CONT/NN,NN2,NN3,NN4,XOLD,HSOL,CONT(5*NVAR)
258      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
259      EXTERNAL MAS, FCN, JAC
260
261C *** *** *** *** *** *** ***
262C  INITIALISATIONS
263C *** *** *** *** *** *** ***
264
265C --------- DUPLIFY N FOR COMMON BLOCK CONT -----
266      NN=N
267      NN2=2*N
268      NN3=3*N
269      NN4=4*N
270
271C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
272      IF(IMPLCT) CALL MAS(N,FMAS,LDMAS)
273
274C ---------- CONSTANTS ---------
275      MBDIAG=MUMAS+1
276      IF (METH.EQ.2) THEN
277C ---------- METHOD WITH GAMMA = 4/15 ---------------
278          GAMMA=4.0D0/15.0D0
279          C2=23.0D0/30.0D0
280          C3=17.0D0/30.0D0
281          C4=2881.0D0/28965.0D0+GAMMA
282          ALPH21=15.0D0/8.0D0
283          ALPH31=1577061.0D0/922880.0D0
284          ALPH32=-23427.0D0/115360.0D0
285          ALPH41=647163682356923881.0D0/2414496535205978880.0D0
286          ALPH42=-593512117011179.0D0/3245291041943520.0D0
287          ALPH43=559907973726451.0D0/1886325418129671.0D0
288          ALPH51=724545451.0D0/796538880.0D0
289          ALPH52=-830832077.0D0/267298560.0D0
290          ALPH53=30957577.0D0/2509272.0D0
291          ALPH54=-69863904375173.0D0/6212571137048.0D0
292          E1=7752107607.0D0/11393456128.0D0
293          E2=-17881415427.0D0/11470078208.0D0
294          E3=2433277665.0D0/179459416.0D0
295          E4=-96203066666797.0D0/6212571137048.0D0
296          D11= 24.74416644927758D0
297          D12= -4.325375951824688D0
298          D13= 41.39683763286316D0
299          D14= -61.04144619901784D0
300          D15= -3.391332232917013D0
301          D21= -51.98245719616925D0
302          D22= 10.52501981094525D0
303          D23= -154.2067922191855D0
304          D24= 214.3082125319825D0
305          D25= 14.71166018088679D0
306          D31= 33.14347947522142D0
307          D32= -19.72986789558523D0
308          D33= 230.4878502285804D0
309          D34= -287.6629744338197D0
310          D35= -18.99932366302254D0
311          D41= -5.905188728329743D0
312          D42= 13.53022403646467D0
313          D43= -117.6778956422581D0
314          D44= 134.3962081008550D0
315          D45= 8.678995715052762D0
316         ETA1=23.D0/8.D0
317         ANU1= 0.9838473040915402D0
318         ANU2= 0.3969226768377252D0
319         AMU1= 0.6563374010466914D0
320         AMU3= 0.3372498196189311D0
321      ELSE
322         PRINT *, 'WRONG  CHOICE OF <METH>'
323      END IF
324      POSNEG=SIGN(1.D0,XEND-X)
325      HMAX1=MIN(ABS(HMAX),ABS(XEND-X))
326      IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
327      H=MIN(ABS(H),HMAX1)
328      H=SIGN(H,POSNEG)
329      HOLD=H
330      CFAC=SAFE*(1+2*NIT)
331      NEWTRE=.FALSE.
332      REJECT=.FALSE.
333      FIRST=.TRUE.
334      FACCO1=1.D0
335      FACCO2=1.D0
336      FACCO3=1.D0
337      FACCO4=1.D0
338      FACCO5=1.D0
339      NSING=0
340      XOLD=X
341      IF (ITOL.EQ.0) THEN
342          DO 8 I=1,N
343   8      SCAL(I)=1.D0 / ( AbsTol(1)+RelTol(1)*DABS(Y(I)) )
344      ELSE
345          DO 9 I=1,N
346   9      SCAL(I)=1.D0 / ( AbsTol(I)+RelTol(I)*DABS(Y(I)) )
347      END IF
348
349C --- BASIC INTEGRATION STEP 
350  10  CONTINUE
351
352C *** *** *** *** *** *** ***
353C  COMPUTATION OF THE JACOBIAN
354C *** *** *** *** *** *** ***
355      NJAC=NJAC+1
356      CALL JAC(N,X,Y,FJAC)
357      CALJAC=.TRUE.
358  20  CONTINUE
359
360C *** *** *** *** *** *** ***
361C  COMPUTE THE MATRIX E AND ITS DECOMPOSITION
362C *** *** *** *** *** *** ***
363      FAC1=1.D0/(H*GAMMA)
364      IF (IMPLCT) THEN
365         print *, 'IMPLCT 4'
366      ELSE  ! EXPLICIT SYSTEM
367C --- THE MATRIX E (MAS=IDENTITY, JACOBIAN A FULL MATRIX)
368c         DO 526 J=1,N
369c           DO 525 I=1,N
370c 525         E(I,J)=-FJAC(I,J)
371c 526       E(J,J)=E(J,J)+FAC1
372c         CALL DEC(N,LE,E,IP,IER)
373          DO K=1,LU_NONZERO
374             E(K) = -FJAC(K)
375          END DO
376          DO I=1,N
377             IDG = LU_DIAG(I)
378             E(IDG) = E(IDG) + FAC1
379          END DO
380          CALL KppDecomp ( E, IER)
381         
382         IF (IER.NE.0) GOTO 79
383      END IF
384      NDEC=NDEC+1
385  30  CONTINUE
386
387      IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
388      XPH=X+H
389C --- LOOP FOR THE 5 STAGES
390      FACCO1=DMAX1(FACCO1,UROUND)**0.8D0
391      FACCO2=DMAX1(FACCO2,UROUND)**0.8D0
392      FACCO3=DMAX1(FACCO3,UROUND)**0.8D0
393      FACCO4=DMAX1(FACCO4,UROUND)**0.8D0
394      FACCO5=DMAX1(FACCO5,UROUND)**0.8D0
395
396C *** *** *** *** *** *** ***
397C  STARTING VALUES FOR NEWTON ITERATION
398C *** *** *** *** *** *** ***
399      DO 59 ISTAGE=1,5
400      IF (ISTAGE.EQ.1) THEN
401          XCH=X+GAMMA*H
402          IF (FIRST.OR.NEWTRE) THEN
403              DO 132 I=1,N
404 132          Z(I)=0.D0
405          ELSE
406              S=1.D0+GAMMA*H/HOLD
407              DO 232 I=1,N
408c  232          Z(I) = 0.D0
409 232          Z(I)=S*(CONT(I+NN)+S*(CONT(I+NN2)+S*(CONT(I+NN3)
410     &             +S*CONT(I+NN4))))-YHAT(I)
411
412          END IF
413          DO 31 I=1,N
414  31      G1(I)=0.D0
415          FACCON=FACCO1
416      END IF
417      IF (ISTAGE.EQ.2) THEN
418          XCH=X+C2*H
419          DO 131 I=1,N
420          Z1I=Z1(I)
421          Z(I)=ETA1*Z1I
422 131      G1(I)=ALPH21*Z1I
423          FACCON=FACCO2
424      END IF
425      IF (ISTAGE.EQ.3) THEN
426          XCH=X+C3*H
427          DO 231 I=1,N
428          Z1I=Z1(I)
429          Z2I=Z2(I)
430          Z(I)=ANU1*Z1I+ANU2*Z2I
431 231      G1(I)=ALPH31*Z1I+ALPH32*Z2I
432          FACCON=FACCO3
433      END IF
434      IF (ISTAGE.EQ.4) THEN
435          XCH=X+C4*H
436          DO 331 I=1,N
437          Z1I=Z1(I)
438          Z3I=Z3(I)
439          Z(I)=AMU1*Z1I+AMU3*Z3I
440 331      G1(I)=ALPH41*Z1I+ALPH42*Z2(I)+ALPH43*Z3I
441          FACCON=FACCO4
442      END IF
443      IF (ISTAGE.EQ.5) THEN
444          XCH=XPH
445          DO 431 I=1,N
446          Z1I=Z1(I)
447          Z2I=Z2(I)
448          Z3I=Z3(I)
449          Z4I=Z4(I)
450          Z(I)=E1*Z1I+E2*Z2I+E3*Z3I+E4*Z4I
451          YHAT(I)=Z(I)
452 431      G1(I)=ALPH51*Z1I+ALPH52*Z2I+ALPH53*Z3I+ALPH54*Z4I
453          FACCON=FACCO5
454      END IF
455
456
457
458C *** *** *** *** *** *** *** *** *** *** ***
459C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
460C *** *** *** *** *** *** *** *** *** *** ***
461            NEWT=0
462            THETA=ABS(THET)
463            IF (REJECT) THETA=2*ABS(THET)
464  40        CONTINUE
465            IF (NEWT.GE.NIT) THEN
466                H=H/2.D0
467                REJECT=.TRUE.
468                NEWTRE=.TRUE.
469                IF (CALJAC) GOTO 20
470                GOTO 10
471            END IF
472
473C ---     COMPUTE THE RIGHT-HAND SIDE
474            DO 41 I=1,N
475            H1(I)=G1(I)-Z(I)
476  41        CONT(I)=Y(I)+Z(I)
477            CALL FCN(N,XCH,CONT,F1)
478            NFCN=NFCN+1
479
480C ---     KppSolve THE LINEAR SYSTEMS
481            IF (IMPLCT) THEN
482                print *, 'IMPLCT 2'
483            ELSE
484                DO 345 I=1,N
485 345            F1(I)=H1(I)*FAC1+F1(I)
486C                CALL SOL(N,LE,E,F1,IP)
487                CALL KppSolve(E, F1)
488            END IF
489            NEWT=NEWT+1
490            DYNO=0.D0
491C --- NORM 2 ---
492            DO 57 I=1,N
493  57        DYNO=DYNO+(F1(I)*SCAL(I))**2
494            DYNO=DSQRT(DYNO/N)
495C --- NORM INF ---
496C            DO 57 I=1,N
497C  57        DYNO=DMAX1( DYNO, DABS(F1(I)*SCAL(I)) )
498
499
500C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
501            IF (NEWT.GE.2.AND.NEWT.LT.NIT) THEN
502                THETA=DYNO/DYNOLD
503                IF (THETA.LT.0.99D0) THEN
504                    FACCON=THETA/(1.0D0-THETA)
505                    DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)
506                    QNEWT=DMAX1(1.0D-4,DMIN1(16.0D0,DYTH/FNEWT))
507                    IF (QNEWT.GE.1.0D0) THEN
508                         H=.8D0*H*QNEWT**(-1.0D0/(NIT-NEWT))
509                         REJECT=.TRUE.
510                         NEWTRE=.TRUE.
511                         IF (CALJAC) GOTO 20
512                         GOTO 10
513                    END IF
514                ELSE
515                    NEWTRE=.TRUE.
516                    GOTO 78
517                END IF
518            END IF
519            DYNOLD=DYNO
520            DO 58 I=1,N
521  58        Z(I)=Z(I)+F1(I)
522            NSOL=NSOL+1
523            IF (FACCON*DYNO.GT.FNEWT) GOTO 40
524
525C --- END OF SIMPILFIED NEWTON
526      IF (ISTAGE.EQ.1) THEN
527          DO I=1,N
528            Z1(I) = Z(I)
529          END DO 
530          FACCO1=FACCON
531      END IF
532      IF (ISTAGE.EQ.2) THEN
533          DO I=1,N
534            Z2(I) = Z(I)
535          END DO 
536          FACCO2=FACCON
537      END IF
538      IF (ISTAGE.EQ.3) THEN
539          DO I=1,N
540            Z3(I) = Z(I)
541          END DO 
542          FACCO3=FACCON
543      END IF
544      IF (ISTAGE.EQ.4) THEN
545          DO I=1,N
546            Z4(I) = Z(I)
547          END DO 
548          FACCO4=FACCON
549      END IF
550      IF (ISTAGE.EQ.5) THEN
551          DO I=1,N
552            Z5(I) = Z(I)
553          END DO 
554          FACCO5=FACCON
555      END IF
556  59  CONTINUE
557
558
559C *** *** *** *** *** *** ***
560C  ERROR ESTIMATION 
561C *** *** *** *** *** *** ***
562      NSTEP=NSTEP+1
563      IF (IMPLCT) THEN
564          print *,'IMPLCT 3'
565      ELSE
566          DO 461 I=1,N
567 461      CONT(I)=FAC1*(Z5(I)-YHAT(I))
568      END IF
569
570      CALL KppSolve(E, CONT)
571
572      ERR=0.D0
573C ---- NORM 2 ---
574      DO 64 I=1,N
575  64  ERR=ERR+(CONT(I)*SCAL(I))**2
576      ERR=DMAX1(DSQRT(ERR/N),1.D-10)
577
578C ---- NORM INF ---
579C      DO 64 I=1,N
580c  64  ERR=DMAX1( ERR, DABS( CONT(I)*SCAL(I) ) )
581
582C --- COMPUTATION OF HNEW
583C --- WE REQUIRE .25<=HNEW/H<=10.
584      FAC=DMIN1(SAFE,CFAC/(NEWT+2*NIT))
585      QUOT=DMAX1(.25D0,DMIN1(10.D0,(ERR)**.25D0/FAC))
586      HNEW= H/QUOT
587
588C *** *** *** *** *** *** ***
589C  IS THE ERROR SMALL ENOUGH ?
590C *** *** *** *** *** *** ***
591      IF (ERR.LT.1.D0) THEN
592C --- STEP IS ACCEPTED 
593         FIRST=.FALSE.
594         NACCPT=NACCPT+1
595         HOLD=H
596         XOLD=X
597C --- COEFFICIENTS FOR CONTINUOUS SOLUTION
598         DO 74 I=1,N
599          Z1I=Z1(I)
600          Z2I=Z2(I)
601          Z3I=Z3(I)
602          Z4I=Z4(I)
603          Z5I=Z5(I)
604         CONT(I)=Y(I)
605         Y(I)=Y(I)+Z5I 
606         CONT(I+NN) =D11*Z1I+D12*Z2I+D13*Z3I+D14*Z4I+D15*Z5I
607         CONT(I+NN2)=D21*Z1I+D22*Z2I+D23*Z3I+D24*Z4I+D25*Z5I
608         CONT(I+NN3)=D31*Z1I+D32*Z2I+D33*Z3I+D34*Z4I+D35*Z5I
609         CONT(I+NN4)=D41*Z1I+D42*Z2I+D43*Z3I+D44*Z4I+D45*Z5I
610         YHAT(I)=Z5I
611         IF (ITOL.EQ.0) THEN
612           SCAL(I)=1.D0/( AbsTol(1)+RelTol(1)*DABS(Y(I)) )
613         ELSE
614           SCAL(I)=1.D0/( AbsTol(I)+RelTol(I)*DABS(Y(I)) )
615         END IF
616  74     CONTINUE
617         X=XPH
618         CALJAC=.FALSE.
619         IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
620            H=HOPT
621            IDID=1
622            RETURN
623         END IF
624         IF (IJAC.EQ.0) CALL FCN(N,X,Y,Y0)
625         NFCN=NFCN+1
626         HNEW=POSNEG*DMIN1(DABS(HNEW),HMAX1)
627         HOPT=HNEW
628         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
629         REJECT=.FALSE.
630         NEWTRE=.FALSE.
631         IF ((X+HNEW/QUOT1-XEND)*POSNEG.GT.0.D0) THEN
632            H=XEND-X
633         ELSE
634            QT=HNEW/H
635            IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30
636            H = HNEW
637         END IF
638         IF (THETA.LE.THET) GOTO 20
639         GOTO 10
640
641      ELSE
642C --- STEP IS REJECTED 
643         REJECT=.TRUE.
644         IF (FIRST) THEN
645             H=H/10.D0
646         ELSE
647             H=HNEW
648         END IF
649         IF (NACCPT.GE.1) NREJCT=NREJCT+1
650         IF (CALJAC) GOTO 20
651         GOTO 10
652      END IF
653
654C --- UNEXPECTED STEP-REJECTION
655  78  CONTINUE
656      IF (IER.NE.0) THEN
657          WRITE (6,*) ' MATRIX IS SINGULAR, IER=',IER,' X=',X,' H=',H
658          NSING=NSING+1
659          IF (NSING.GE.6) GOTO 79
660      END IF
661      H=H*0.5D0
662      REJECT=.TRUE.
663      IF (CALJAC) GOTO 20
664      GOTO 10
665
666C --- FAIL EXIT
667  79  WRITE (6,979) X,H,IER
668 979  FORMAT(' EXIT OF SDIRK4 AT X=',D14.7,'   H=',D14.7,'   IER=',I4)
669      IDID=-1
670      RETURN
671      END
672C
673
674 
675      SUBROUTINE FUNC_CHEM(N, T, Y, P)
676      INCLUDE 'KPP_ROOT_Parameters.h'
677      INCLUDE 'KPP_ROOT_Global.h'
678      INTEGER N
679      KPP_REAL   T, Told
680      KPP_REAL   Y(NVAR), P(NVAR)
681      Told = TIME
682      TIME = T
683      CALL Update_SUN()
684      CALL Update_RCONST()
685      CALL Fun( Y,  FIX, RCONST, P )
686      TIME = Told
687      RETURN
688      END
689
690 
691      SUBROUTINE JAC_CHEM(N, T, Y, J)
692      INCLUDE 'KPP_ROOT_Parameters.h'
693      INCLUDE 'KPP_ROOT_Global.h'
694      INTEGER N
695      KPP_REAL   Told, T
696      KPP_REAL   Y(NVAR), J(LU_NONZERO)
697      Told = TIME
698      TIME = T
699      CALL Update_SUN()
700      CALL Update_RCONST()
701      CALL Jac_SP( Y,  FIX, RCONST, J )
702      TIME = Told
703      RETURN
704      END                                                                                                                 
705
Note: See TracBrowser for help on using the repository browser.