source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/ros3.f @ 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: 8.3 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 ROS3(NVAR,TIN,TOUT,STEPMIN,STEPMAX,
18     +                   STEPMIN,VAR,ATOL,RTOL,
19     +                   Info,FUNC_CHEM,JAC_CHEM)
20
21      RETURN
22      END
23
24     
25      SUBROUTINE ROS3(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       L-stable Rosenbrock 3(2), with 
33C     strongly A-stable 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, April 1996
63C     The Center for Global and Regional Environmental Research
64
65      KPP_REAL K1(NVAR), K2(NVAR), K3(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, Tplus, H, elo
71      KPP_REAL ERR, factor, facmax
72      KPP_REAL gam, c21, c31, c32, b1, b2, b3
73      KPP_REAL d1, d2, d3, a21, a31, a32
74      KPP_REAL alpha2, alpha3, g1, g2, g3
75      KPP_REAL tau, x1, x2, x3, dround, ytol
76      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
77      INTEGER    Info(5)
78      LOGICAL    IsReject,Autonomous
79      EXTERNAL   FUNC_CHEM, JAC_CHEM
80     
81      gam=   .43586652150845899941601945119356d+00
82      c21=  -.10156171083877702091975600115545d+01
83      c31=   .40759956452537699824805835358067d+01
84      c32=   .92076794298330791242156818474003d+01
85       b1=   .10000000000000000000000000000000d+01
86       b2=   .61697947043828245592553615689730d+01
87       b3=  -.42772256543218573326238373806514d+00
88       d1=   .50000000000000000000000000000000d+00
89       d2=  -.29079558716805469821718236208017d+01
90       d3=   .22354069897811569627360909276199d+00
91       a21 = 1.d0
92       a31 = 1.d0
93       a32 = 0.d0
94       alpha2 = gam
95       alpha3 = gam
96       g1=   .43586652150845899941601945119356d+00
97       g2=   .24291996454816804366592249683314d+00
98       g3=   .21851380027664058511513169485832d+01
99
100
101c     Initialization of counters, etc.
102      Autonomous = Info(1) .EQ. 1
103      uround = 1.d-15
104      dround = DSQRT(uround)
105      IF (Hmax.le.0.D0) THEN
106          Hmax = DABS(Tnext-T)
107      END IF     
108      H = DMAX1(1.d-8, Hstart)
109      Tplus = T
110      IsReject = .false.
111      Naccept  = 0
112      Nreject  = 0
113      Nfcn     = 0
114      Njac     = 0
115
116C === Starting the time loop ===     
117 10    continue 
118 
119       Tplus = T + H
120       if ( Tplus .gt. Tnext ) then
121          H = Tnext - T
122          Tplus = Tnext
123       end if
124
125       CALL JAC_CHEM(NVAR, T, y, JAC)
126       Njac = Njac+1
127       gHinv = -1.0d0/(gam*H)
128       do 15 j=1,LU_NONZERO
129         JAC(j) = -JAC(j) 
130 15    continue
131       do 20 j=1,NVAR
132         JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) - gHinv
133 20    continue
134       CALL KppDecomp (JAC, ier)
135
136       if (ier.ne.0) then
137         if ( H.gt.Hmin) then
138            H = 5.0d-1*H
139            go to 10
140         else
141            print *,'IER <> 0, H=',H
142            stop
143         end if     
144       end if 
145       
146       CALL FUNC_CHEM(NVAR, T, y, F1)
147
148C ====== NONAUTONOMOUS CASE ===============
149       IF (.not. Autonomous) THEN
150         tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T)
151         CALL FUNC_CHEM(NVAR, T+tau, y, K2)
152         nfcn=nfcn+1
153         do 30 j = 1,NVAR
154           K3(j) = ( K2(j)-F1(j) )/tau
155 30      continue
156 
157C ----- STAGE 1 (NONAUTONOMOUS) -----
158         x1 = g1*H
159         do 35 j = 1,NVAR
160           K1(j) =  F1(j) + x1*K3(j)
161 35      continue
162         CALL KppSolve (JAC, K1)
163     
164C ----- STAGE 2 (NONAUTONOMOUS) -----
165       do 40 j = 1,NVAR
166         ynew(j) = y(j) + K1(j) 
167 40    continue
168       CALL FUNC_CHEM(NVAR, T+gam*H, ynew, F1)
169       nfcn=nfcn+1
170       x1 = c21/H
171       x2 = g2*H
172       do 45 j = 1,NVAR
173         K2(j) = F1(j) + x1*K1(j) + x2*K3(j)
174 45    continue
175       CALL KppSolve (JAC, K2)
176       
177C ----- STAGE 3  (NONAUTONOMOUS) -----
178       x1 = c31/H
179       x2 = c32/H
180       x3 = g3*H
181       do 50 j = 1,NVAR
182         K3(j) = F1(j) + x1*K1(j) + x2*K2(j) + x3*K3(j)
183 50    continue
184       CALL KppSolve (JAC, K3)
185
186
187C ====== AUTONOMOUS CASE ===============
188       ELSE
189
190C ----- STAGE 1 (AUTONOMOUS) -----
191         do 60 j = 1,NVAR
192           K1(j) =  F1(j) 
193 60      continue
194         CALL KppSolve (JAC, K1)
195     
196C ----- STAGE 2 (AUTONOMOUS) -----
197       do 65 j = 1,NVAR
198         ynew(j) = y(j) + K1(j) 
199 65    continue
200       CALL FUNC_CHEM(NVAR, T + gam*H, ynew, F1)
201       nfcn=nfcn+1
202         x1 = c21/H
203         do 70 j = 1,NVAR
204           K2(j) = F1(j) + x1*K1(j) 
205 70      continue
206         CALL KppSolve (JAC, K2)
207       
208C ----- STAGE 3  (AUTONOMOUS) -----
209       x1 = c31/H
210       x2 = c32/H
211       do 90 j = 1,NVAR
212         K3(j) = F1(j) + x1*K1(j) + x2*K2(j)
213 90    continue
214       CALL KppSolve (JAC, K3)
215
216       END  IF ! Autonomousous
217
218C ---- The Solution ---
219
220       do 120 j = 1,NVAR
221         ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) 
222 120   continue
223
224
225C ====== Error estimation ========
226
227        ERR=0.d0
228        do 130 i=1,NVAR
229           ytol = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(ynew(i)))
230           ERR=ERR+((d1*K1(i)+d2*K2(i)+d3*K3(i))/ytol)**2
231 130    continue     
232        ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
233
234C ======= Choose the stepsize ===============================
235 
236        elo    = 3.0D0 ! estimator local order
237        factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
238        Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
239 
240C ======= Rejected/Accepted Step ============================
241 
242        IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
243          IsReject = .true.
244          H = DMIN1(H/10,Hnew)
245          Nreject  = Nreject+1
246        ELSE
247          DO 140 i=1,NVAR
248             y(i)  = ynew(i)
249 140      CONTINUE
250          T = Tplus
251          IF (.NOT.IsReject) THEN
252              H = Hnew   ! Do not increase stepsize if previos step was rejected
253          END IF   
254          IsReject = .false.
255          Naccept = Naccept+1
256        END IF
257
258
259C ======= End of the time loop ===============================
260      if ( T .lt. Tnext ) go to 10
261     
262     
263C ======= Output Information =================================
264      Info(2) = Nfcn
265      Info(3) = Njac
266      Info(4) = Naccept
267      Info(5) = Nreject
268      Hstart  = H
269     
270      RETURN
271      END       
272
273 
274 
275      SUBROUTINE FUNC_CHEM(N, T, Y, P)
276      INCLUDE 'KPP_ROOT_params.h'
277      INCLUDE 'KPP_ROOT_global.h'
278      INTEGER N
279      KPP_REAL   T, Told
280      KPP_REAL   Y(NVAR), P(NVAR)
281      Told = TIME
282      TIME = T
283      CALL Update_SUN()
284      CALL Update_RCONST()
285      CALL Fun( Y,  FIX, RCONST, P )
286      TIME = Told
287      RETURN
288      END
289 
290      SUBROUTINE JAC_CHEM(N, T, Y, J)
291      INCLUDE 'KPP_ROOT_params.h'
292      INCLUDE 'KPP_ROOT_global.h'
293      INTEGER N
294      KPP_REAL   Told, T
295      KPP_REAL   Y(NVAR), J(LU_NONZERO)
296      Told = TIME
297      TIME = T
298      CALL Update_SUN()
299      CALL Update_RCONST()
300      CALL Jac_SP( Y,  FIX, RCONST, J )
301      TIME = Told
302      RETURN
303      END                                                                                                                 
304 
Note: See TracBrowser for help on using the repository browser.