source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/ros4.f @ 3288

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

Merge of branch palm4u into trunk

File size: 10.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 ROS4(NVAR,TIN,TOUT,STEPMIN,STEPMAX,
18     +                   STEPMIN,VAR,ATOL,RTOL,
19     +                   Info,FUNC_CHEM,JAC_CHEM)
20
21      RETURN
22      END
23 
24
25 
26 
27      SUBROUTINE ROS4(N,T,Tnext,Hmin,Hmax,Hstart,
28     +                   y,AbsTol,RelTol,
29     +                   Info,FUNC_CHEM,JAC_CHEM)
30      IMPLICIT NONE
31      INCLUDE 'KPP_ROOT_params.h'
32      INCLUDE 'KPP_ROOT_sparse.h'                                                                                                 
33C
34C  Four Stages, Fourth Order L-stable Rosenbrock Method, 
35C     with embedded L-stable, third order method for error control
36C     Simplified version of E. Hairer's atmros4; the coefficients are slightly
37C     different
38C
39C
40C  INPUT ARGUMENTS:
41C     y = Vector of (NVAR) concentrations, contains the
42C         initial values on input
43C     [T, Tnext] = the integration interval
44C     Hmin, Hmax = lower and upper bounds for the selected step-size.
45C          Note that for Step = Hmin the current computed
46C          solution is unconditionally accepted by the error
47C          control mechanism.
48C     AbsTol, RelTol = (NVAR) dimensional vectors of
49C          componentwise absolute and relative tolerances.
50C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
51C          See the header below.
52C     JAC_CHEM = name of routine that computes the Jacobian, in
53C          sparse format. KPP syntax. See the header below.
54C     Info(1) = 1  for  Autonomous   system
55C             = 0  for nonAutonomous system
56C
57C  OUTPUT ARGUMENTS:
58C     y = the values of concentrations at Tend.
59C     T = equals TENDon output.
60C     Info(2) = # of FUNC_CHEM CALLs.
61C     Info(3) = # of JAC_CHEM CALLs.
62C     Info(4) = # of accepted steps.
63C     Info(5) = # of rejected steps.
64C     Hstart = The last accepted stepsize
65C
66C    Adrian Sandu, December 2001
67C
68      KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR), K4(NVAR)
69      KPP_REAL F1(NVAR)
70      KPP_REAL DFDT(NVAR)
71      KPP_REAL JAC(LU_NONZERO)
72      KPP_REAL Hmin,Hmax,Hstart
73      KPP_REAL y(NVAR), ynew(NVAR)
74      KPP_REAL AbsTol(NVAR), RelTol(NVAR)
75      KPP_REAL T, Tnext, H, Hnew, Tplus
76      KPP_REAL elo,ghinv,uround
77      KPP_REAL ERR, factor, facmax
78      KPP_REAL w, e, dround, tau
79      KPP_REAL hgam1, hgam2, hgam3, hgam4
80      KPP_REAL hc21, hc31, hc32, hc41, hc42, hc43
81     
82      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
83      INTEGER    Info(5)
84      LOGICAL    IsReject, Autonomous
85      EXTERNAL   FUNC_CHEM, JAC_CHEM                                                                                               
86
87
88C   The method coefficients
89      DOUBLE PRECISION gamma, gamma2, gamma3, gamma4
90      PARAMETER ( gamma  =  0.5728200000000000D+00 )
91      PARAMETER ( gamma2 = -0.1769193891319233D+01 )
92      PARAMETER ( gamma3 =  0.7592633437920482D+00 )
93      PARAMETER ( gamma4 = -0.1049021087100450D+00 )
94      DOUBLE PRECISION a21, a31, a32, a41, a42, a43
95      PARAMETER (  a21 = 0.2000000000000000D+01  )   
96      PARAMETER (  a31 = 0.1867943637803922D+01  )   
97      PARAMETER (  a32 = 0.2344449711399156D+00  )   
98      DOUBLE PRECISION alpha2, alpha3
99      PARAMETER (  alpha2  = 0.1145640000000000D+01  )   
100      PARAMETER (  alpha3  = 0.6552168638155900D+00  )   
101      DOUBLE PRECISION c21, c31, c32, c41, c42, c43
102      PARAMETER (  c21  = -0.7137615036412310D+01  )
103      PARAMETER (  c31  =  0.2580708087951457D+01  )
104      PARAMETER (  c32  =  0.6515950076447975D+00  )
105      PARAMETER (  c41  = -0.2137148994382534D+01  )
106      PARAMETER (  c42  = -0.3214669691237626D+00  )
107      PARAMETER (  c43  = -0.6949742501781779D+00  )
108      DOUBLE PRECISION b1, b2, b3, b4
109      PARAMETER (  b1 = 0.2255570073418735D+01  )
110      PARAMETER (  b2 = 0.2870493262186792D+00  )
111      PARAMETER (  b3 = 0.4353179431840180D+00  )
112      PARAMETER (  b4 = 0.1093502252409163D+01 ) 
113      DOUBLE PRECISION d1, d2, d3, d4
114      PARAMETER (  d1 = -0.2815431932141155D+00  )
115      PARAMETER (  d2 = -0.7276199124938920D-01  )
116      PARAMETER (  d3 = -0.1082196201495311D+00  )
117      PARAMETER (  d4 = -0.1093502252409163D+01  ) 
118 
119c     Initialization of counters, etc.
120      Autonomous = Info(1) .EQ. 1
121      uround = 1.d-15
122      dround = DSQRT(uround)
123      IF (Hmax.le.0.D0) THEN
124          Hmax = DABS(Tnext-T)
125      END IF     
126      H = DMAX1(1.d-8, Hstart)
127      Tplus = T
128      IsReject = .false.
129      Naccept  = 0
130      Nreject  = 0
131      Nfcn     = 0
132      Njac     = 0
133 
134C === Starting the time loop ===
135 10    CONTINUE
136 
137       Tplus = T + H
138       IF ( Tplus .gt. Tnext ) THEN
139          H = Tnext - T
140          Tplus = Tnext
141       END IF
142 
143C              Initial Function and Jacobian values                                                                                                                           
144       CALL FUNC_CHEM( T, y, F1 )
145       CALL JAC_CHEM( T, y, JAC )
146 
147C              The time derivative for non-Autonomous case                                                                                                                           
148       IF (.not. Autonomous) THEN
149         tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T)
150         CALL FUNC_CHEM( T+tau, y, K2 )
151         nfcn=nfcn+1
152         DO 20 j = 1,NVAR
153           DFDT(j) = ( K2(j)-F1(j) )/tau
154 20      CONTINUE
155       END IF
156 
157C              Form the Prediction matrix and compute its LU factorization                                                                                                                           
158       Njac = Njac+1
159       ghinv = 1.0d0/(gamma*H)
160       DO 30 j=1,LU_NONZERO
161         JAC(j) = -JAC(j)
162 30    CONTINUE
163       DO 40 j=1,NVAR
164         JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
165 40    CONTINUE
166       CALL KppDecomp (JAC, ier)
167C
168       IF (ier.ne.0) THEN
169         IF ( H.gt.Hmin) THEN
170            H = 5.0d-1*H
171            GO TO 10
172         ELSE
173            PRINT *,'ROS4: Singular factorization at T=',T,'; H=',H
174            STOP
175         END IF
176       END IF
177                                                                                                                           
178 
179C ------------ STAGE 1-------------------------
180       DO 50 j = 1,NVAR
181         K1(j) =  F1(j)
182 50    CONTINUE
183       IF (.NOT. Autonomous) THEN
184          hgam1 = H*gamma
185          DO 60 j=1,NVAR
186            K1(j) = K1(j) + hgam1*DFDT(j)           
187 60       CONTINUE
188       END IF
189       CALL KppSolve (JAC, K1)
190 
191C ----------- STAGE 2 -------------------------
192       DO 70 j = 1,NVAR
193         ynew(j) = y(j) + a21*K1(j)
194 70    CONTINUE
195       CALL FUNC_CHEM( T+alpha2*H, ynew, F1)
196       nfcn=nfcn+1
197       hc21 = c21/H
198       DO 80 j = 1,NVAR
199         K2(j) = F1(j) + hc21*K1(j)
200 80    CONTINUE
201       IF (.NOT. Autonomous) THEN
202         hgam2 = H*gamma2
203         DO 90 j=1,NVAR
204            K2(j) = K2(j) + hgam2*DFDT(j)           
205 90      CONTINUE
206       END IF
207       CALL KppSolve (JAC, K2)
208 
209 
210C ------------ STAGE 3 -------------------------
211       DO 100 j = 1,NVAR
212         ynew(j) = y(j) + a31*K1(j) + a32*K2(j)
213 100   CONTINUE
214       CALL FUNC_CHEM( T+alpha3*H, ynew, F1)
215       nfcn=nfcn+1
216       hc31 = c31/H
217       hc32 = c32/H
218       DO 110 j = 1,NVAR
219         K3(j) = F1(j) + hc31*K1(j) + hc32*K2(j)
220 110   CONTINUE
221       IF (.NOT. Autonomous) THEN
222         hgam3 = H*gamma3
223         DO 120 j=1,NVAR
224            K3(j) = K3(j) + hgam3*DFDT(j)           
225 120     CONTINUE
226       END IF
227       CALL KppSolve (JAC, K3)
228 
229C ------------ STAGE 4 -------------------------
230C              Note: uses the same function value as stage 3
231       hc41 = c41/H
232       hc42 = c42/H
233       hc43 = c43/H
234       DO 140 j = 1,NVAR
235         K4(j) = F1(j) + hc41*K1(j) + hc42*K2(j) + hc43*K3(j)
236 140   CONTINUE
237       IF (.NOT. Autonomous) THEN
238          hgam4 = H*gamma4
239          DO 150 j=1,NVAR
240            K4(j) = K4(j) + hgam4*DFDT(j)
241 150      CONTINUE
242       END IF
243       CALL KppSolve (JAC, K4)
244 
245
246
247C ---- The Solution ---
248       DO 160 j = 1,NVAR
249         ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) + b4*K4(j)
250 160   CONTINUE
251 
252 
253C ====== Error estimation ========
254 
255        ERR=0.d0
256        DO 170 j = 1,NVAR
257           w = AbsTol(j) + RelTol(j)*DMAX1(DABS(y(j)),DABS(ynew(j)))
258           e = d1*K1(j) + d2*K2(j) + d3*K3(j) + d4*K4(j)
259           ERR = ERR + ( e/w )**2
260 170    CONTINUE
261        ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
262 
263C ======= Choose the stepsize ===============================
264     
265        elo    = 4.0D0 ! estimator local order
266        factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
267        Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
268 
269C ======= Rejected/Accepted Step ============================
270 
271        IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
272          IsReject = .true.
273          H = DMIN1(H/10,Hnew)
274          Nreject  = Nreject+1
275        ELSE
276          DO 180 i=1,NVAR
277             y(i)  = ynew(i)
278 180      CONTINUE
279          T = Tplus
280          IF (.NOT.IsReject) THEN
281              H = Hnew   ! Do not increase stepsize if previos step was rejected
282          END IF   
283          IsReject = .false.
284          Naccept = Naccept+1
285        END IF
286 
287C ======= End of the time loop ===============================
288      IF ( T .lt. Tnext ) GO TO 10
289 
290C ======= Output Information =================================
291      Info(2) = Nfcn
292      Info(3) = Njac
293      Info(4) = Naccept
294      Info(5) = Nreject
295      Hstart  = H
296 
297      RETURN
298      END
299 
300 
301      SUBROUTINE FUNC_CHEM( T, Y, P )
302      INCLUDE 'KPP_ROOT_params.h'
303      INCLUDE 'KPP_ROOT_global.h'
304      KPP_REAL   T, Told
305      KPP_REAL   Y(NVAR), P(NVAR)
306      Told = TIME
307      TIME = T
308      CALL Update_SUN()
309      CALL Update_RCONST()
310      CALL Fun( Y,  FIX, RCONST, P )
311      TIME = Told
312      RETURN
313      END
314
315 
316      SUBROUTINE JAC_CHEM( T, Y, J )
317      INCLUDE 'KPP_ROOT_params.h'
318      INCLUDE 'KPP_ROOT_global.h'
319      KPP_REAL   Told, T
320      KPP_REAL   Y(NVAR), J(LU_NONZERO)
321      Told = TIME
322      TIME = T
323      CALL Update_SUN()
324      CALL Update_RCONST()
325      CALL Jac_SP( Y,  FIX, RCONST, J )
326      TIME = Told
327      RETURN
328      END                                                                                                                 
329
330
331
332
333
334                                                                                                                           
Note: See TracBrowser for help on using the repository browser.