source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/DVODE/example1.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.6 KB
Line 
1    MODULE EXAMPLE1
2
3! DEMONSTRATION PROGRAM FOR THE DVODE_F90 PACKAGE.
4
5! The following is a simple example problem, with the coding
6! needed for its solution by DVODE_F90. The problem is from
7! chemical kinetics, and consists of the following three rate
8! 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
15! The following coding solves this problem with DVODE_F90,
16! using a user supplied Jacobian and printing results at
17! t = .4, 4.,...,4.d10. It uses ITOL = 2 and ATOL much smaller
18! for y2 than y1 or y3 because y2 has much smaller values. At
19! the end of the run, statistical quantities of interest are
20! printed. (See optional output in the full DVODE description
21! below.) Output is written to the file example1.dat.
22
23    CONTAINS
24
25      SUBROUTINE FEX(NEQ,T,Y,YDOT)
26        IMPLICIT NONE
27        INTEGER NEQ
28        DOUBLE PRECISION T, Y, YDOT
29        DIMENSION Y(NEQ), YDOT(NEQ)
30
31        YDOT(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)
32        YDOT(3) = 3.E7*Y(2)*Y(2)
33        YDOT(2) = -YDOT(1) - YDOT(3)
34        RETURN
35      END SUBROUTINE FEX
36
37      SUBROUTINE JEX(NEQ,T,Y,ML,MU,PD,NRPD)
38        IMPLICIT NONE
39        INTEGER NEQ, ML, MU, NRPD
40        DOUBLE PRECISION PD, T, Y
41        DIMENSION Y(NEQ), PD(NRPD,NEQ)
42
43        PD(1,1) = -.04D0
44        PD(1,2) = 1.D4*Y(3)
45        PD(1,3) = 1.D4*Y(2)
46        PD(2,1) = .04D0
47        PD(2,3) = -PD(1,3)
48        PD(3,2) = 6.E7*Y(2)
49        PD(2,2) = -PD(1,2) - PD(3,2)
50        RETURN
51      END SUBROUTINE JEX
52
53    END MODULE EXAMPLE1
54
55!******************************************************************
56
57    PROGRAM RUNEXAMPLE1
58
59      USE DVODE_F90_M
60      USE EXAMPLE1
61
62      IMPLICIT NONE
63      DOUBLE PRECISION ATOL, RTOL, T, TOUT, Y, RSTATS
64      INTEGER NEQ, ITASK, ISTATE, ISTATS, IOUT, IERROR, I
65      DIMENSION Y(3), ATOL(3), RSTATS(22), ISTATS(31)
66
67      TYPE (VODE_OPTS) :: OPTIONS
68
69      OPEN (UNIT=6,FILE='example1.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.D-4
78      ATOL(1) = 1.D-8
79      ATOL(2) = 1.D-14
80      ATOL(3) = 1.D-6
81      ITASK = 1
82      ISTATE = 1
83!     OPTIONS = SET_OPTS(DENSE_J=.TRUE.,ABSERR_VECTOR=ATOL,RELERR=RTOL, &
84!       USER_SUPPLIED_JACOBIAN=.TRUE.)
85      OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE.,ABSERR_VECTOR=ATOL,      &
86        RELERR=RTOL,USER_SUPPLIED_JACOBIAN=.TRUE.)
87      DO IOUT = 1, 12
88        CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JEX)
89        CALL GET_STATS(RSTATS,ISTATS)
90        WRITE (6,90003) T, Y(1), Y(2), Y(3)
91        DO I = 1, NEQ
92          IF (Y(I)<0.0D0) IERROR = 1
93        END DO
94        IF (ISTATE<0) THEN
95          WRITE (6,90004) ISTATE
96          STOP
97        END IF
98        TOUT = TOUT*10.0D0
99      END DO
100      WRITE (6,90000) ISTATS(11), ISTATS(12), ISTATS(13), ISTATS(19), &
101        ISTATS(20), ISTATS(21), ISTATS(22)
102      IF (IERROR==1) THEN
103        WRITE (6,90001)
104      ELSE
105        WRITE (6,90002)
106      END IF
10790000 FORMAT (/'  No. steps =',I4,'   No. f-s =',I4,'  No. J-s =',I4, &
108        '   No. LU-s =',I4/'  No. nonlinear iterations =', &
109        I4/'  No. nonlinear convergence failures =', &
110        I4/'  No. error test failures =',I4/)
11190001 FORMAT (/' An error occurred.')
11290002 FORMAT (/' No errors occurred.')
11390003 FORMAT (' At t =',D12.4,'   y =',3D14.6)
11490004 FORMAT (///' Error halt: ISTATE =',I3)
115      STOP
116    END PROGRAM RUNEXAMPLE1
Note: See TracBrowser for help on using the repository browser.