source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/oldies/ros1.f @ 3698

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

Merge of branch palm4u into trunk

File size: 4.6 KB
Line 
1      SUBROUTINE INTEGRATE( TIN, TOUT )
2 
3      INCLUDE 'KPP_ROOT_params.h'
4      INCLUDE 'KPP_ROOT_global.h'
5
6C TIN - Start Time
7      KPP_REAL TIN
8C TOUT - End Time
9      KPP_REAL TOUT
10
11      INTEGER    INFO(5)
12
13      EXTERNAL FUNC_CHEM, JAC_CHEM
14
15      INFO(1) = Autonomous
16 
17      CALL ROS1(NVAR,TIN,TOUT,STEPMIN,VAR,
18     +                   Info,FUNC_CHEM,JAC_CHEM)
19 
20
21      RETURN
22      END
23 
24
25 
26 
27      SUBROUTINE ROS1(N,T,Tnext,Hstart,
28     +                   y,Info,FUNC_CHEM,JAC_CHEM)
29     
30      INCLUDE 'KPP_ROOT_params.h'
31      INCLUDE 'KPP_ROOT_sparse.h'                                                                                                 
32C
33C  Linearly Implicit Euler
34C  A method of theoretical interest but of no practical value
35C
36C  INPUT ARGUMENTS:
37C     y = Vector of (NVAR) concentrations, contains the
38C         initial values on input
39C     [T, Tnext] = the integration interval
40C     Hmin, Hmax = lower and upper bounds for the selected step-size.
41C          Note that for Step = Hmin the current computed
42C          solution is unconditionally accepted by the error
43C          control mechanism.
44C     AbsTol, RelTol = (NVAR) dimensional vectors of
45C          componentwise absolute and relative tolerances.
46C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
47C          See the header below.
48C     JAC_CHEM = name of routine that computes the Jacobian, in
49C          sparse format. KPP syntax. See the header below.
50C     Info(1) = 1  for  Autonomous   system
51C             = 0  for nonAutonomous system
52C
53C  OUTPUT ARGUMENTS:
54C     y = the values of concentrations at Tend.
55C     T = equals TENDon output.
56C     Info(2) = # of FUNC_CHEM CALLs.
57C     Info(3) = # of JAC_CHEM CALLs.
58C     Info(4) = # of accepted steps.
59C     Info(5) = # of rejected steps.
60C     Hstart = The last accepted stepsize
61C 
62C    Adrian Sandu, December 2001
63C 
64      KPP_REAL Fv(NVAR)
65      KPP_REAL JAC(LU_NONZERO)
66      KPP_REAL H, Hstart
67      KPP_REAL y(NVAR)
68      KPP_REAL T, Tnext, Tplus
69      KPP_REAL elo,ghinv,uround
70     
71      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j
72      INTEGER    Info(5)
73      LOGICAL    IsReject, Autonomous
74      EXTERNAL    FUNC_CHEM, JAC_CHEM                                                                                               
75
76
77      H     = Hstart
78      Tplus = T
79      Nfcn  = 0
80      Njac  = 0
81 
82C === Starting the time loop ===
83 10    CONTINUE
84 
85       Tplus = T + H
86       IF ( Tplus .gt. Tnext ) THEN
87          H = Tnext - T
88          Tplus = Tnext
89       END IF
90 
91C              Initial Function and Jacobian values                                                                                                                           
92       CALL FUNC_CHEM(NVAR, T, y, Fv)
93       Nfcn = Nfcn+1
94       CALL JAC_CHEM(NVAR, T, y, JAC)
95       Njac = Njac+1
96 
97C              Form the Prediction matrix and compute its LU factorization                                                                                                                           
98       DO 40 j=1,NVAR
99         JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) - 1.0d0/H
100 40    CONTINUE
101       CALL KppDecomp (JAC, ier)
102C 
103       IF (ier.ne.0) THEN
104          PRINT *,'ROS1: Singular factorization at T=',T,'; H=',H
105          STOP
106       END IF
107 
108C ------------ STAGE 1-------------------------
109       CALL KppSolve (JAC, Fv)
110 
111C ---- The Solution ---
112       DO 160 j = 1,NVAR
113         y(j) = y(j) - Fv(j) 
114 160   CONTINUE
115       T = T + H
116 
117C ======= End of the time loop ===============================
118      IF ( T .lt. Tnext ) GO TO 10
119 
120C ======= Output Information =================================
121      Info(2) = Nfcn
122      Info(3) = Njac
123      Info(4) = Njac
124      Info(5) = 0
125      Hstart  = H
126 
127      RETURN
128      END
129 
130 
131      SUBROUTINE FUNC_CHEM(N, T, Y, P)
132      INCLUDE 'KPP_ROOT_params.h'
133      INCLUDE 'KPP_ROOT_global.h'
134      INTEGER N
135      KPP_REAL   T, Told
136      KPP_REAL   Y(NVAR), P(NVAR)
137      Told = TIME
138      TIME = T
139      CALL Update_SUN()
140      CALL Update_RCONST()
141      CALL Fun( Y,  FIX, RCONST, P )
142      TIME = Told
143      RETURN
144      END
145
146 
147      SUBROUTINE JAC_CHEM(N, T, Y, J)
148      INCLUDE 'KPP_ROOT_params.h'
149      INCLUDE 'KPP_ROOT_global.h'
150      INTEGER N
151      KPP_REAL   Told, T
152      KPP_REAL   Y(NVAR), J(LU_NONZERO)
153      Told = TIME
154      TIME = T
155      CALL Update_SUN()
156      CALL Update_RCONST()
157      CALL Jac_SP( Y,  FIX, RCONST, J )
158      TIME = Told
159      RETURN
160      END                                                                                                                 
161
162
163
164
165
166                                                                                                                           
Note: See TracBrowser for help on using the repository browser.