source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/DVODE/example2.f90 @ 2697

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

Merge of branch palm4u into trunk

File size: 3.7 KB
Line 
1    MODULE EXAMPLE2
2
3! DEMONSTRATION PROGRAM FOR THE DVODE_F90 PACKAGE.
4
5! The following is a modification of the previous example
6! program to illustrate root finding. The problem is from
7! chemical kinetics, and consists of the following three
8! rate equations:
9!     dy1/dt = -.04d0*y1 + 1.d4*y2*y3
10!     dy2/dt = .04d0*y1 - 1.d4*y2*y3 - 3.d7*y2**2
11!     dy3/dt = 3.d7*y2**2
12! on the interval from t = 0.0d0 to t = 4.d10, with initial
13! conditions y1 = 1.0d0, y2 = y3 = 0.0d0. The problem is stiff.
14! In addition, we want to find the values of t, y1, y2,
15! and y3 at which:
16!   (1) y1 reaches the value 1.d-4, and
17!   (2) y3 reaches the value 1.d-2.
18
19! The following coding solves this problem with DVODE_F90
20! using an internally generated dense Jacobian and
21! printing results at t = .4, 4., ..., 4.d10, and at the
22! computed roots. It uses ITOL = 2 and ATOL much smaller
23! for y2 than y1 or y3 because y2 has much smaller values.
24! At the end of the run, statistical quantities of interest
25! are printed (see optional outputs in the full description
26! below). Output is written to the file example2.dat.
27
28    CONTAINS
29
30      SUBROUTINE FEX(NEQ,T,Y,YDOT)
31        IMPLICIT NONE
32        INTEGER NEQ
33        DOUBLE PRECISION T, Y, YDOT
34        DIMENSION Y(3), YDOT(3)
35
36        YDOT(1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3)
37        YDOT(3) = 3.0D7*Y(2)*Y(2)
38        YDOT(2) = -YDOT(1) - YDOT(3)
39        RETURN
40      END SUBROUTINE FEX
41
42      SUBROUTINE GEX(NEQ,T,Y,NG,GOUT)
43        IMPLICIT NONE
44        INTEGER NEQ, NG
45        DOUBLE PRECISION T, Y, GOUT
46        DIMENSION Y(3), GOUT(2)
47
48        GOUT(1) = Y(1) - 1.0D-4
49        GOUT(2) = Y(3) - 1.0D-2
50        RETURN
51      END SUBROUTINE GEX
52
53    END MODULE EXAMPLE2
54
55!******************************************************************
56
57    PROGRAM RUNEXAMPLE2
58
59      USE DVODE_F90_M
60      USE EXAMPLE2
61
62      IMPLICIT NONE
63      INTEGER ITASK, ISTATE, NG, NEQ, IOUT, JROOT, ISTATS, IERROR, I
64      DOUBLE PRECISION ATOL, RTOL, RSTATS, T, TOUT, Y
65      DIMENSION Y(3), ATOL(3), RSTATS(22), ISTATS(31), JROOT(2)
66
67      TYPE (VODE_OPTS) :: OPTIONS
68
69      OPEN (UNIT=6,FILE='example2.dat')
70      IERROR = 0
71      NEQ = 3
72      Y(1) = 1.0D0
73      Y(2) = 0.0D0
74      Y(3) = 0.0D0
75      T = 0.0D0
76      TOUT = 0.4D0
77      RTOL = 1.0D-4
78      ATOL(1) = 1.0D-8
79      ATOL(2) = 1.0D-10
80      ATOL(3) = 1.0D-8
81      ITASK = 1
82      ISTATE = 1
83      NG = 2
84!     OPTIONS = SET_OPTS(DENSE_J=.TRUE.,RELERR=RTOL,ABSERR_VECTOR=ATOL, &
85!       NEVENTS=NG)
86      OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE.,RELERR=RTOL,             &
87        ABSERR_VECTOR=ATOL,NEVENTS=NG)
88      DO 20 IOUT = 1, 12
8910      CONTINUE
90        CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,G_FCN=GEX)
91        CALL GET_STATS(RSTATS,ISTATS,NG,JROOT)
92        WRITE (6,90000) T, Y(1), Y(2), Y(3)
93        DO I = 1, NEQ
94          IF (Y(I)<0.0D0) IERROR = 1
95        END DO
9690000   FORMAT (' At t =',D12.4,'   Y =',3D14.6)
97        IF (ISTATE<0) GO TO 30
98        IF (ISTATE==2) GO TO 20
99        WRITE (6,90001) JROOT(1), JROOT(2)
10090001   FORMAT (5X,' The above line is a root, JROOT =',2I5)
101        ISTATE = 2
102        GO TO 10
10320    TOUT = TOUT*10.0D0
104      WRITE (6,90002) ISTATS(11), ISTATS(12), ISTATS(13), ISTATS(10)
105      IF (IERROR==1) THEN
106        WRITE (6,90003)
107      ELSE
108        WRITE (6,90004)
109      END IF
11090002 FORMAT (/' No. steps =',I4,'  No. f-s =',I4,'  No. J-s =',I4, &
111        '  No. g-s =',I4/)
112      STOP
11330    WRITE (6,90005) ISTATE
11490003 FORMAT (/' An error occurred.')
11590004 FORMAT (/' No errors occurred.')
11690005 FORMAT (///' Error halt.. ISTATE =',I3)
117      STOP
118    END PROGRAM RUNEXAMPLE2
Note: See TracBrowser for help on using the repository browser.