source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/rodas3.f @ 3423

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

Merge of branch palm4u into trunk

File size: 7.2 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        CALL RODAS3(NVAR,TIN,TOUT,STEPMIN,STEPMAX,STEPMIN,
17     +                   VAR,ATOL,RTOL,
18     +                   Info,FUNC_CHEM,JAC_CHEM)
19
20      RETURN
21      END
22
23
24
25      SUBROUTINE RODAS3(N,T,Tnext,Hmin,Hmax,Hstart,
26     +                   y,AbsTol,RelTol,
27     +                   Info,FUNC_CHEM,JAC_CHEM)
28      IMPLICIT NONE
29      INCLUDE 'KPP_ROOT_params.h'
30      INCLUDE 'KPP_ROOT_sparse.h'
31
32C     Stiffly accurate Rosenbrock 3(2), with 
33C     stiffly accurate embedded formula for error control. 
34C
35C     All the arguments aggree with the KPP syntax.
36C
37C  INPUT ARGUMENTS:
38C     y = Vector of (NVAR) concentrations, contains the
39C         initial values on input
40C     [T, Tnext] = the integration interval
41C     Hmin, Hmax = lower and upper bounds for the selected step-size.
42C          Note that for Step = Hmin the current computed
43C          solution is unconditionally accepted by the error
44C          control mechanism.
45C     AbsTol, RelTol = (NVAR) dimensional vectors of 
46C          componentwise absolute and relative tolerances.
47C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
48C          See the header below.
49C     JAC_CHEM = name of routine that computes the Jacobian, in
50C          sparse format. KPP syntax. See the header below.
51C     Info(1) = 1  for  autonomous   system
52C             = 0  for nonautonomous system
53C
54C  OUTPUT ARGUMENTS:
55C     y = the values of concentrations at Tend.
56C     T = equals Tend on output.
57C     Info(2) = # of FUNC_CHEM calls.
58C     Info(3) = # of JAC_CHEM calls.
59C     Info(4) = # of accepted steps.
60C     Info(5) = # of rejected steps.
61C   
62C     Adrian Sandu, March 1996
63C     The Center for Global and Regional Environmental Research
64
65      KPP_REAL   K1(NVAR), K2(NVAR), K3(NVAR), K4(NVAR)
66      KPP_REAL   F1(NVAR), JAC(LU_NONZERO)
67      KPP_REAL   Hmin,Hmax,Hnew,Hstart,ghinv,uround
68      KPP_REAL   y(NVAR), ynew(NVAR)
69      KPP_REAL   AbsTol(NVAR), RelTol(NVAR)
70      KPP_REAL   T, Tnext, H, Hold, Tplus
71      KPP_REAL   ERR, factor, facmax
72      KPP_REAL   c43, tau, x1, x2, ytol, elo
73     
74      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
75      INTEGER    Info(5)
76      LOGICAL    IsReject,Autonomous
77      EXTERNAL   FUNC_CHEM, JAC_CHEM
78
79c     Initialization of counters, etc.
80      Autonomous = Info(1) .EQ. 1
81      uround = 1.d-15
82      c43 = - 8.d0/3.d0 
83      H = DMAX1(1.d-8, Hstart)
84      Hmin = DMAX1(Hmin,uround*(Tnext-T))
85      Hmax = DMIN1(Hmax,Tnext-T)
86      Tplus = T
87      IsReject = .false.
88      Naccept  = 0
89      Nreject  = 0
90      Nfcn     = 0
91      Njac     = 0
92
93
94C === Starting the time loop ===     
95 10    continue 
96       Tplus = T + H
97       if ( Tplus .gt. Tnext ) then
98          H = Tnext - T
99          Tplus = Tnext
100       end if
101
102       CALL JAC_CHEM(NVAR, T, y, JAC)
103       Njac = Njac+1
104       gHinv = -2.0d0/H
105       do 20 j=1,NVAR
106         JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + gHinv
107 20    continue
108       CALL KppDecomp (JAC, ier)
109
110       if (ier.ne.0) then
111         if ( H.gt.Hmin) then
112            H = 5.0d-1*H
113            go to 10
114         else
115            print *,'IER <> 0, H=',H
116            stop
117         end if     
118       end if 
119       
120       CALL FUNC_CHEM(NVAR, T, y, F1)
121
122C ====== NONAUTONOMOUS CASE ===============
123       IF (.not. Autonomous) THEN
124         tau = DSQRT( uround*DMAX1( 1.0d-5, DABS(T) ) )
125         CALL FUNC_CHEM(NVAR, T+tau, y, K2)
126         nfcn=nfcn+1
127         do 30 j = 1,NVAR
128           K3(j) = ( K2(j)-F1(j) )/tau
129 30      continue
130 
131C ----- STAGE 1 (NONAUTONOMOUS) -----
132         x1 = 0.5*H
133         do 40 j = 1,NVAR
134           K1(j) =  F1(j) + x1*K3(j) 
135 40      continue
136         CALL KppSolve (JAC, K1)
137
138C ----- STAGE 2 (NONAUTONOMOUS) -----
139         x1 = 4.d0/H
140         x2 = 1.5d0*H
141         do 50 j = 1,NVAR
142           K2(j) = F1(j) - x1*K1(j) + x2*K3(j)
143 50      continue
144         CALL KppSolve (JAC, K2)
145
146C ====== AUTONOMOUS CASE ===============
147       ELSE
148C ----- STAGE 1 (AUTONOMOUS) -----
149         do 60 j = 1,NVAR
150           K1(j) =  F1(j) 
151 60      continue
152         CALL KppSolve (JAC, K1)
153     
154C ----- STAGE 2 (AUTONOMOUS) -----
155         x1 = 4.d0/H
156         do 70 j = 1,NVAR
157           K2(j) = F1(j) - x1*K1(j) 
158 70      continue
159         CALL KppSolve (JAC, K2)
160       END IF
161       
162C ----- STAGE 3 -----
163       do 80 j = 1,NVAR
164         ynew(j) = y(j) - 2.0d0*K1(j) 
165 80    continue
166       CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
167       nfcn=nfcn+1
168       do 90 j = 1,NVAR
169         K3(j) = F1(j) + ( -K1(j) + K2(j) )/H
170 90    continue
171       CALL KppSolve (JAC, K3)
172
173C ----- STAGE 4 -----
174       do 100 j = 1,NVAR
175         ynew(j) = y(j) - 2.0d0*K1(j) - K3(j)
176 100   continue
177       CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
178       nfcn=nfcn+1
179       do 110 j = 1,NVAR
180         K4(j) = F1(j) + ( -K1(j) + K2(j) - C43*K3(j)  )/H
181 110   continue
182       CALL KppSolve (JAC, K4)
183
184C ---- The Solution ---
185
186       do 120 j = 1,NVAR
187         ynew(j) = y(j) - 2.0d0*K1(j) - K3(j) - K4(j) 
188 120   continue
189
190
191C ====== Error estimation ========
192
193        ERR=0.d0
194        do 130 i=1,NVAR
195           ytol = AbsTol(i) + RelTol(i)*DABS(ynew(i))
196           ERR = ERR + ( K4(i)/ytol )**2
197 130    continue     
198        ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
199
200C ======= Choose the stepsize ===============================
201        elo    = 3.0D0 ! estimator local order
202        factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
203        Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
204 
205C ======= Rejected/Accepted Step ============================
206 
207        IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
208          IsReject = .true.
209          H = DMIN1(H/10,Hnew)
210          Nreject  = Nreject+1
211        ELSE
212          DO 140 i=1,NVAR
213             y(i)  = ynew(i)
214 140      CONTINUE
215          T = Tplus
216          IF (.NOT.IsReject) THEN
217              H = Hnew   ! Do not increase stepsize if previos step was rejected
218          END IF   
219          IsReject = .false.
220          Naccept = Naccept+1
221        END IF
222
223
224C ======= End of the time loop ===============================
225      if ( T .lt. Tnext ) go to 10
226 
227     
228     
229C ======= Output Information =================================
230      Info(2) = Nfcn
231      Info(3) = Njac
232      Info(4) = Naccept
233      Info(5) = Nreject
234      Hstart  = H
235     
236      RETURN
237      END       
238
239 
240 
241      SUBROUTINE FUNC_CHEM(N, T, Y, P)
242      INCLUDE 'KPP_ROOT_params.h'
243      INCLUDE 'KPP_ROOT_global.h'
244      INTEGER N
245      KPP_REAL   T, Told
246      KPP_REAL   Y(NVAR), P(NVAR)
247      Told = TIME
248      TIME = T
249      CALL Update_SUN()
250      CALL Update_RCONST()
251      CALL Fun( Y, FIX, RCONST, P )
252      TIME = Told
253      RETURN
254      END
255
256 
257      SUBROUTINE JAC_CHEM(N, T, Y, J)
258      INCLUDE 'KPP_ROOT_params.h'
259      INCLUDE 'KPP_ROOT_global.h'
260      INTEGER N
261      KPP_REAL   Told, T
262      KPP_REAL   Y(NVAR), J(LU_NONZERO)
263      Told = TIME
264      TIME = T
265      CALL Update_SUN()
266      CALL Update_RCONST()
267      CALL Jac_SP( Y, FIX, RCONST, J )
268      TIME = Told
269      RETURN
270      END                                                                                                                 
271 
Note: See TracBrowser for help on using the repository browser.