source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/ros2.f90 @ 3947

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

Merge of branch palm4u into trunk

File size: 7.5 KB
Line 
1      SUBROUTINE INTEGRATE( TIN, TOUT )
2 
3      USE KPP_ROOT_global
4
5! TIN - Start Time
6      KPP_REAL TIN
7! TOUT - End Time
8      KPP_REAL TOUT
9
10      INTEGER    INFO(5)
11
12      EXTERNAL FUNC_CHEM, JAC_CHEM
13
14      INFO(1) = Autonomous
15 
16      call ROS2(NVAR,TIN,TOUT,STEPMIN,STEPMAX,   &
17            STEPMIN,VAR,ATOL,RTOL, &
18            Info,FUNC_CHEM,JAC_CHEM)
19 
20
21      END SUBROUTINE INTEGRATE
22 
23
24 
25 
26      SUBROUTINE ROS2(N,T,Tnext,Hmin,Hmax,Hstart,  &
27            y,AbsTol,RelTol,      &
28            Info,FUNC_CHEM,JAC_CHEM)
29     
30      USE KPP_ROOT_params
31      USE KPP_ROOT_Jacobian_sparsity                                 
32      IMPLICIT NONE
33
34!  INPUT ARGUMENTS:
35!     y = Vector of (NVAR) concentrations, contains the
36!     initial values on input
37!     [T, Tnext] = the integration interval
38!     Hmin, Hmax = lower and upper bounds for the selected step-size.
39!      Note that for Step = Hmin the current computed
40!      solution is unconditionally accepted by the error
41!      control mechanism.
42!     AbsTol, RelTol = (NVAR) dimensional vectors of
43!      componentwise absolute and relative tolerances.
44!     FUNC_CHEM = name of routine of derivatives. KPP syntax.
45!      See the header below.
46!     JAC_CHEM = name of routine that computes the Jacobian, in
47!      sparse format. KPP syntax. See the header below.
48!     Info(1) = 1  for  autonomous   system
49!         = 0  for nonautonomous system
50!     Info(2) = 1  for third order embedded formula
51!         = 0  for first order embedded formula
52!
53!   Note: Stage 3 used to build strongly A-stable order 3 formula for error control
54!        Embed3 = (Info(2).EQ.1)
55!     if Embed3 = .true. then the third order embedded formula is used
56!         .false. then a first order embedded  formula is used
57!
58!
59!  OUTPUT ARGUMENTS:
60!     y = the values of concentrations at Tend.
61!     T = equals Tend on output.
62!     Info(2) = # of FUNC_CHEM calls.
63!     Info(3) = # of JAC_CHEM calls.
64!     Info(4) = # of accepted steps.
65!     Info(5) = # of rejected steps.
66 
67      KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR)
68      KPP_REAL F1(NVAR), JAC(LU_NONZERO)
69      KPP_REAL DFDT(NVAR)
70      KPP_REAL Hmin,Hmax,Hnew,Hstart,ghinv,uround
71      KPP_REAL y(NVAR), ynew(NVAR)
72      KPP_REAL AbsTol(NVAR), RelTol(NVAR)
73      KPP_REAL T, Tnext, H, Hold, Tplus
74      KPP_REAL ERR, factor, facmax
75      KPP_REAL tau, beta, elo, dround, a21, c31, c32
76      KPP_REAL gamma3, d1, d2, d3, gam
77      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
78      INTEGER    Info(5)
79      LOGICAL    IsReject, Autonomous, Embed3
80      EXTERNAL    FUNC_CHEM, JAC_CHEM                                 
81
82      KPP_REAL gamma, m1, m2, alpha, beta1, beta2, delta, w, e
83 
84!     Initialization of counters, etc.
85      Autonomous = Info(1) .EQ. 1
86      Embed3  = Info(2) .EQ. 1
87      uround  = 1.d-15
88      dround  = dsqrt(uround)
89      H = DMAX1(1.d-8, Hmin)
90      Tplus = T
91      IsReject = .false.
92      Naccept  = 0
93      Nreject  = 0
94      Nfcn     = 0
95      Njac     = 0
96     
97!     Method Parameters     
98      gamma = 1.d0 + 1.d0/sqrt(2.d0)
99      a21   = - 1.d0/gamma
100      m1 = -3.d0/(2.d0*gamma)
101      m2 = -1.d0/(2.d0*gamma)
102      c31 = -1.0D0/gamma**2*(1.0D0-7.0D0*gamma+9.0D0*gamma**2) &
103      /(-1.0D0+2.0D0*gamma)
104      c32 = -1.0D0/gamma**2*(1.0D0-6.0D0*gamma+6.0D0*gamma**2) &
105       /(-1.0D0+2.0D0*gamma)/2
106      gamma3 = 0.5D0 - 2*gamma
107      d1 = ((-9.0D0*gamma+8.0D0*gamma**2+2.0D0)/gamma**2/      &
108       (-1.0D0+2*gamma))/6.0D0
109      d2 = ((-1.0D0+3.0D0*gamma)/gamma**2/(-1.0D0+2.0D0*gamma))/6.0D0
110      d3 = -1.0D0/(3.0D0*gamma)
111 
112! === Starting the time loop ===
113 10    CONTINUE
114      Tplus = T + H
115      if ( Tplus .gt. Tnext ) then
116         H = Tnext - T
117         Tplus = Tnext
118      end if
119 
120      call JAC_CHEM(NVAR, T, y, JAC)
121 
122      Njac = Njac+1
123      ghinv = -1.0d0/(gamma*H)
124      DO j=1,NVAR
125         JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
126      END DO
127      CALL KppDecomp (JAC, ier)
128 
129      if (ier.ne.0) then
130        if ( H.gt.Hmin) then
131          H = 5.0d-1*H
132          go to 10
133        else
134          print *,'IER <> 0, H=',H
135          stop
136        end if
137      end if
138 
139      call FUNC_CHEM(NVAR, T, y, F1)
140 
141! ====== NONAUTONOMOUS CASE ===============
142   IF (.not. Autonomous) THEN
143     tau = dsign(dround*dmax1( 1.0d-6, dabs(T) ), T)
144     call FUNC_CHEM(NVAR, T+tau, y, K2)
145     nfcn=nfcn+1
146     DO j = 1,NVAR
147       DFDT(j) = ( K2(j)-F1(j) )/tau
148     END DO
149   END IF ! .NOT.Autonomous
150 
151! ----- STAGE 1 -----
152   DO j = 1,NVAR
153      K1(j) =  F1(j) 
154   END DO
155   IF (.NOT.Autonomous) THEN
156     delta = gamma*H
157     DO j = 1,NVAR
158       K1(j) = K1(j) + delta*DFDT(j)
159     END DO     
160   END IF ! .NOT.Autonomous
161   call KppSolve (JAC, K1)
162 
163! ----- STAGE 2  -----
164   DO j = 1,NVAR
165     ynew(j) = y(j) + a21*K1(j)
166   END DO
167   call FUNC_CHEM(NVAR, T+H, ynew, F1)
168   nfcn=nfcn+1
169   beta = 2.d0/(gamma*H)
170   DO j = 1,NVAR
171     K2(j) = F1(j) + beta*K1(j) 
172   END DO
173   IF (.NOT. Autonomous) THEN
174     delta = -gamma*H
175     DO j = 1,NVAR
176       K2(j) = K2(j) + delta*DFDT(j)
177     END DO   
178   END IF ! .NOT.Autonomous
179   call KppSolve (JAC, K2)
180 
181! ----- STAGE 3  -----
182   IF (Embed3) THEN
183   beta1 = -c31/H
184   beta2 = -c32/H
185   delta = gamma3*H
186   DO j = 1,NVAR
187     K3(j) = F1(j) + beta1*K1(j) + beta2*K2(j) 
188   END DO
189   IF (.NOT.Autonomous) THEN
190     DO j = 1,NVAR
191       K3(j) = K3(j) + delta*DFDT(j)
192     END DO   
193   END IF ! .NOT.Autonomous
194   CALL KppSolve (JAC, K3)
195    END IF ! Embed3
196 
197
198
199! ---- The Solution ---
200   DO j = 1,NVAR
201     ynew(j) = y(j) + m1*K1(j) + m2*K2(j)
202   END DO
203 
204 
205! ====== Error estimation ========
206 
207    ERR=0.d0
208    DO i=1,NVAR
209       w = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(ynew(i)))
210    IF ( Embed3 ) THEN
211     e = d1*K1(i) + d2*K2(i) + d3*K3(i)
212    ELSE
213         e = 1.d0/(2.d0*gamma)*(K1(i)+K2(i))
214    END IF  ! Embed3
215       ERR = ERR + ( e/w )**2
216    END DO
217    ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
218 
219! ======= Choose the stepsize ===============================
220 
221    IF ( Embed3 ) THEN
222        elo = 3.0D0 ! estimator local order
223    ELSE
224        elo = 2.0D0
225    END IF   
226    factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
227    Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
228 
229! ======= Rejected/Accepted Step ============================
230 
231    IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
232      IsReject = .true.
233      H = DMIN1(H/10,Hnew)
234      Nreject  = Nreject+1
235    ELSE
236      DO i=1,NVAR
237         y(i)  = ynew(i)
238      END DO
239      T = Tplus
240      IF (.NOT. IsReject) THEN
241          H = Hnew   ! Do not increase stepsize if previous step was rejected
242      END IF   
243      IsReject = .false.
244      Naccept = Naccept+1
245    END IF
246
247 
248! ======= End of the time loop ===============================
249      IF ( T .lt. Tnext ) GO TO 10
250 
251 
252 
253! ======= Output Information =================================
254      Info(2) = Nfcn
255      Info(3) = Njac
256      Info(4) = Naccept
257      Info(5) = Nreject
258 
259      END SUBROUTINE Ros2
260 
261 
262 
263      SUBROUTINE FUNC_CHEM(N, T, Y, P)     
264      USE KPP_ROOT_global
265      INTEGER N 
266      KPP_REAL   T, Told
267      KPP_REAL   Y(NVAR), P(NVAR)
268      Told = TIME
269      TIME = T
270      CALL Update_SUN()
271      CALL Update_RCONST()
272      CALL Fun( Y, FIX, RCONST, P )
273      TIME = Told
274      END SUBROUTINE FUNC_CHEM
275
276 
277      SUBROUTINE JAC_CHEM(N, T, Y, J)     
278      USE KPP_ROOT_global
279      INTEGER N 
280      KPP_REAL   Told, T
281      KPP_REAL   Y(NVAR), J(LU_NONZERO)
282      Told = TIME
283      TIME = T
284      CALL Update_SUN()
285      CALL Update_RCONST()
286      CALL Jac_SP( Y, FIX, RCONST, J )
287      TIME = Told
288      END SUBROUTINE JAC_CHEM                                           
289 
290
291
292
293
294                                                       
Note: See TracBrowser for help on using the repository browser.