source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/ros1_ddm.f @ 2696

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

Merge of branch palm4u into trunk

File size: 9.4 KB
Line 
1      SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )
2 
3      INCLUDE 'KPP_ROOT_params.h'
4      INCLUDE 'KPP_ROOT_global.h'
5
6      INTEGER   NSENSIT
7C TIN - Start Time
8      KPP_REAL TIN
9C TOUT - End Time
10      KPP_REAL TOUT
11C TOUT - End Time
12      KPP_REAL Y( NVAR*(NSENSIT+1) )
13C ---  Note: Y contains: (1:NVAR) concentrations, followed by
14C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
15C ---                   etc.,  followed by         
16C ---                   (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
17
18      INTEGER    INFO(5)
19
20      EXTERNAL FUNC_CHEM, JAC_CHEM,  HESS_CHEM
21
22      INFO(1) = Autonomous
23 
24      CALL ROS1_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,Y,
25     +                   Info,FUNC_CHEM,JAC_CHEM,HESS_CHEM)
26 
27
28      RETURN
29      END
30 
31
32 
33 
34      SUBROUTINE ROS1_DDM(N,NSENSIT,T,Tnext,Hstart,
35     +                   y,Info,FUNC_CHEM,JAC_CHEM,HESS_CHEM)
36      IMPLICIT NONE
37      INCLUDE 'KPP_ROOT_params.h'
38      INCLUDE 'KPP_ROOT_sparse.h'                                                                                                 
39      INCLUDE 'KPP_ROOT_global.h'                                                                                                 
40C
41C  Linearly Implicit Euler with direct-decoupled calculation of sensitivities
42C  A method of theoretical interest but of no practical value
43C
44C The global variable DDMTYPE distinguishes between:
45C      DDMTYPE = 0 : sensitivities w.r.t. initial values
46C      DDMTYPE = 1 : sensitivities w.r.t. parameters
47C
48C  INPUT ARGUMENTS:
49C     y = Vector of:   (1:NVAR) concentrations, followed by
50C                      (1:NVAR) sensitivities w.r.t. first parameter, followed by
51C                       etc.,  followed by         
52C                      (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
53C         (y contains initial values at input, final values at output)
54C     [T, Tnext] = the integration interval
55C     Hmin, Hmax = lower and upper bounds for the selected step-size.
56C          Note that for Step = Hmin the current computed
57C          solution is unconditionally accepted by the error
58C          control mechanism.
59C     AbsTol, RelTol = (NVAR) dimensional vectors of
60C          componentwise absolute and relative tolerances.
61C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
62C          See the header below.
63C     JAC_CHEM = name of routine that computes the Jacobian, in
64C          sparse format. KPP syntax. See the header below.
65C     Info(1) = 1  for  Autonomous   system
66C             = 0  for nonAutonomous system
67C
68C  OUTPUT ARGUMENTS:
69C     y = the values of concentrations at Tend.
70C     T = equals TENDon output.
71C     Info(2) = # of FUNC_CHEM CALLs.
72C     Info(3) = # of JAC_CHEM CALLs.
73C     Info(4) = # of accepted steps.
74C     Info(5) = # of rejected steps.
75C     Hstart = The last accepted stepsize
76C 
77C    Adrian Sandu, December 2001
78C 
79      INTEGER   NSENSIT
80      KPP_REAL Fv(NVAR*(NSENSIT+1)), Hv(NVAR)
81      KPP_REAL DFDP(NVAR*NSENSIT)
82      KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
83      KPP_REAL HESS(NHESS)
84      KPP_REAL DJDP(NVAR*NSENSIT)
85      KPP_REAL H, Hstart
86      KPP_REAL y(NVAR*(NSENSIT+1))
87      KPP_REAL T, Tnext, Tplus
88      KPP_REAL elo,ghinv,uround
89     
90      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
91      INTEGER    Info(5)
92      LOGICAL    IsReject, Autonomous
93      EXTERNAL   FUNC_CHEM, JAC_CHEM, HESS_CHEM                                                                                               
94
95
96      H     = Hstart
97      Tplus = T
98      Nfcn  = 0
99      Njac  = 0
100 
101C === Starting the time loop ===
102 10    CONTINUE
103 
104       Tplus = T + H
105       IF ( Tplus .gt. Tnext ) THEN
106          H = Tnext - T
107          Tplus = Tnext
108       END IF
109 
110C              Initial Function and Jacobian values                                                                                                                           
111       CALL FUNC_CHEM(NVAR, T, y, Fv)
112       Nfcn = Nfcn+1
113       CALL JAC_CHEM(NVAR, T, y, JAC)
114       Njac = Njac+1
115       CALL HESS_CHEM( NVAR, T, y, HESS )
116       IF (DDMTYPE .EQ. 1) THEN
117          CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
118       END IF 
119 
120C              Form the Prediction matrix and compute its LU factorization                                                                                                                           
121       DO 40 j=1,LU_NONZERO
122         AJAC(j) = -JAC(j)
123 40    CONTINUE
124       DO 50 j=1,NVAR
125         AJAC(LU_DIAG(j)) = AJAC(LU_DIAG(j)) + 1.0d0/H
126 50    CONTINUE
127       CALL KppDecomp (AJAC, ier)
128C 
129       IF (ier.ne.0) THEN
130          PRINT *,'ROS1: Singular factorization at T=',T,'; H=',H
131          STOP
132       END IF
133 
134C ------------ STAGE 1-------------------------
135       CALL KppSolve (AJAC, Fv)
136C               --- If  derivative w.r.t. parameters
137         IF (DDMTYPE .EQ. 1) THEN
138            CALL DJACDPAR(NVAR, NSENSIT, T, y, Fv(1), DJDP)
139         END IF
140C               --- End of derivative w.r.t. parameters         
141
142       DO 100 i=1,NSENSIT
143          CALL Jac_SP_Vec (JAC, y(i*NVAR+1), Fv(i*NVAR+1))
144          CALL Hess_Vec ( HESS, y(i*NVAR+1), Fv(1), Hv )
145          IF (DDMTYPE .EQ. 0) THEN
146            DO 80 j=1,NVAR
147              Fv(i*NVAR+j) = Fv(i*NVAR+j) + Hv(j)
148 80         CONTINUE       
149          ELSE
150            DO 90 j=1,NVAR
151              Fv(i*NVAR+j) = Fv(i*NVAR+j) + Hv(j) 
152     &                + DFDP(i*NVAR+j)+ DJDP((i-1)*NVAR+j)
153 90         CONTINUE       
154          END IF 
155          CALL KppSolve (AJAC, Fv(i*NVAR+1))
156 100   CONTINUE
157 
158C ---- The Solution ---
159       DO 160 j = 1,NVAR*(NSENSIT+1)
160         y(j) = y(j) + Fv(j) 
161 160   CONTINUE
162       T = T + H
163 
164C ======= End of the time loop ===============================
165      IF ( T .lt. Tnext ) THEN
166          GO TO 10
167      END IF     
168 
169C ======= Output Information =================================
170      Info(2) = Nfcn
171      Info(3) = Njac
172      Info(4) = Naccept
173      Info(5) = Nreject
174      Hstart  = H
175 
176      RETURN
177      END
178 
179 
180 
181      SUBROUTINE FUNC_CHEM(N, T, Y, P)
182      INCLUDE 'KPP_ROOT_params.h'
183      INCLUDE 'KPP_ROOT_global.h'
184      INTEGER N
185      KPP_REAL   T, Told
186      KPP_REAL   Y(NVAR), P(NVAR)
187      Told = TIME
188      TIME = T
189      CALL Update_SUN()
190      CALL Update_RCONST()
191      CALL Fun( Y,  FIX, RCONST, P )
192      TIME = Told
193      RETURN
194      END
195
196 
197      SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
198C ---  Computes the partial derivatives of FUNC_CHEM w.r.t. parameters 
199      INCLUDE 'KPP_ROOT_params.h'
200      INCLUDE 'KPP_ROOT_global.h'
201C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients   
202      INTEGER NCOEFF, JCOEFF(NREACT)
203      COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
204     
205      INTEGER N
206      KPP_REAL   T, Told
207      KPP_REAL   Y(NVAR), P(NVAR*NSENSIT)
208      Told = TIME
209      TIME = T
210      CALL Update_SUN()
211      CALL Update_RCONST()
212C
213      IF (DDMTYPE .EQ. 0) THEN
214C ---  Note: the values below are for sensitivities w.r.t. initial values;
215C ---       they may have to be changed for other applications
216        DO j=1,NSENSIT
217          DO i=1,NVAR
218            P(i+NVAR*(j-1)) = 0.0D0
219          END DO
220        END DO
221      ELSE     
222C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
223C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
224C ---       w.r.t. which one differentiates
225        CALL dFun_dRcoeff( Y,  FIX, NCOEFF, JCOEFF, P )
226      END IF
227      TIME = Told
228      RETURN
229      END
230 
231      SUBROUTINE JAC_CHEM(N, T, Y, J)
232      INCLUDE 'KPP_ROOT_params.h'
233      INCLUDE 'KPP_ROOT_global.h'
234      INTEGER N
235      KPP_REAL   Told, T
236      KPP_REAL   Y(NVAR), J(LU_NONZERO)
237      Told = TIME
238      TIME = T
239      CALL Update_SUN()
240      CALL Update_RCONST()
241      CALL Jac_SP( Y,  FIX, RCONST, J )
242      TIME = Told
243      RETURN
244      END                                                                                                                 
245
246 
247      SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
248C ---  Computes the partial derivatives of JAC w.r.t. parameters times user vector U
249      INCLUDE 'KPP_ROOT_params.h'
250      INCLUDE 'KPP_ROOT_global.h'
251C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients   
252      INTEGER NCOEFF, JCOEFF(NREACT)
253      COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
254     
255      INTEGER N
256      KPP_REAL   T, Told
257      KPP_REAL   Y(NVAR), U(NVAR)
258      KPP_REAL   P(NVAR*NSENSIT)
259      Told = TIME
260      TIME = T
261      CALL Update_SUN()
262      CALL Update_RCONST()
263C
264      IF (DDMTYPE .EQ. 0) THEN
265C ---  Note: the values below are for sensitivities w.r.t. initial values;
266C ---       they may have to be changed for other applications
267        DO j=1,NSENSIT
268          DO i=1,NVAR
269            P(i+NVAR*(j-1)) = 0.0D0
270          END DO
271        END DO
272      ELSE     
273C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
274C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
275C ---       w.r.t. which one differentiates
276        CALL dJac_dRcoeff( Y,  FIX, U, NCOEFF, JCOEFF, P )
277      END IF
278      TIME = Told
279      RETURN
280      END
281
282 
283      SUBROUTINE HESS_CHEM(N, T, Y, HESS)
284      INCLUDE 'KPP_ROOT_params.h'
285      INCLUDE 'KPP_ROOT_global.h'
286      INTEGER N
287      KPP_REAL   Told, T
288      KPP_REAL   Y(NVAR), HESS(NHESS)
289      Told = TIME
290      TIME = T
291      CALL Update_SUN()
292      CALL Update_RCONST()
293      CALL Hessian( Y, FIX, RCONST, HESS )
294      TIME = Told
295      RETURN
296      END                                                                                                                 
297
298
299
300                                                                                                                           
Note: See TracBrowser for help on using the repository browser.