|
3983
|
1 SUBROUTINE ODESSA_STESA (NEQ, Y, NROW, NCOL, YH, WM, IWM, EWT, |
|
|
2 1 SAVF, ACOR, PAR, NRS, F, JAC, DF, PJAC, PDF, SOLVE) |
|
|
3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
|
|
4 EXTERNAL F, JAC, DF, PJAC, PDF, SOLVE |
|
|
5 DIMENSION NEQ(*), Y(NROW,*), YH(NROW,NCOL,*), WM(*), IWM(*), |
|
|
6 1 EWT(NROW,*), SAVF(*), ACOR(NROW,*), PAR(*), NRS(*) |
|
|
7 PARAMETER (ONE=1.0D0,ZERO=0.0D0) |
|
|
8 COMMON /ODE001/ ROWND, ROWNS(173), |
|
|
9 1 TESCO(3,12), RDUM1, EL0, H, RDUM2(4), TN, RDUM3, |
|
|
10 2 IOWND1(14), IOWNS(4), |
|
|
11 3 IALTH, LMAX, IDUM1, IERPJ, IERSL, JCUR, IDUM2, KFLAG, L, IDUM3, |
|
|
12 4 MITER, IDUM4(4), N, NQ, IDUM5, NFE, IDUM6(2) |
|
|
13 COMMON /ODE002/ DUPS, DSMS, DDNS, |
|
|
14 1 IOWND2(3), IDUM7, NSV, IDUM8(2), IDF, IDUM9, JOPT, KFLAGS |
|
|
15 C----------------------------------------------------------------------- |
|
|
16 C ODESSA_STESA IS CALLED BY ODESSA_STODE TO PERFORM AN EXPLICIT |
|
|
17 C CALCULATION FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS DY(I)/DP(J), |
|
|
18 C I = 1,N; J = 1,NPAR. |
|
|
19 C |
|
|
20 C IN ADDITION TO THE VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION |
|
|
21 C WITH ODESSA_STESA USES THE FOLLOWING.. |
|
|
22 C Y = AN NROW (=N) BY NCOL (=NSV) REAL ARRAY CONTAINING THE |
|
|
23 C CORRECTED DEPENDENT VARIABLES ON OUTPUT.. |
|
|
24 C Y(I,1) , I = 1,N = STATE VARIABLES (INPUT); |
|
|
25 C Y(I,J) , I = 1,N , J = 2,NSV , |
|
|
26 C = SENSITIVITY COEFFICIENTS, DY(I)/DP(J). |
|
|
27 C YH = AN N BY NSV BY LMAX REAL ARRAY CONTAINING THE PREDICTED |
|
|
28 C DEPENDENT VARIABLES AND THEIR APPROXIMATE SCALED DERIVATIVES. |
|
|
29 C SAVF = A REAL ARRAY OF LENGTH N USED TO STORE FIRST DERIVATIVES |
|
|
30 C OF DEPENDENT VARIABLES IF MITER = 2 OR 5. |
|
|
31 C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING THE EQUATION |
|
|
32 C PARAMETERS OF INTEREST. |
|
|
33 C NRS = AN INTEGER ARRAY OF LENGTH NPAR + 1 CONTAINING THE NUMBER |
|
|
34 C OF REPEATED STEPS (KFLAGS .LT. 0) DUE TO THE SENSITIVITY |
|
|
35 C CALCULATIONS.. |
|
|
36 C NRS(1) = TOTAL NUMBER OF REPEATED STEPS |
|
|
37 C NRS(I) , I = 2,NPAR = NUMBER OF REPEATED STEPS DUE |
|
|
38 C TO PARAMETER I. |
|
|
39 C NSV = NUMBER OF SOLUTION VECTORS = NPAR + 1. |
|
|
40 C KFLAGS = LOCAL ERROR TEST FLAG, = 0 IF TEST PASSES, .LT. 0 IF TEST |
|
|
41 C FAILS, AND STEP NEEDS TO BE REPEATED. ERROR TEST IS APPLIED |
|
|
42 C TO EACH SOLUTION VECTOR INDEPENDENTLY. |
|
|
43 C DUPS, DSMS, DDNS = REAL SCALARS USED FOR COMPUTING RHUP, RHSM, RHDN, |
|
|
44 C ON RETURN TO ODESSA_STODE (IALTH .EQ. 1). |
|
|
45 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, IALTH, LMAX, |
|
|
46 C IERPJ, IERSL, JCUR, KFLAG, L, MITER, N, NQ, NFE, AND JOPT. |
|
|
47 C----------------------------------------------------------------------- |
|
|
48 DUPS = ZERO |
|
|
49 DSMS = ZERO |
|
|
50 DDNS = ZERO |
|
|
51 HL0 = H*EL0 |
|
|
52 EL0I = ONE/EL0 |
|
|
53 TI2 = ONE/TESCO(2,NQ) |
|
|
54 TI3 = ONE/TESCO(3,NQ) |
|
|
55 C IF MITER = 2 OR 5 (OR IDF = 0), SUPPLY DERIVATIVES AT CORRECTED |
|
|
56 C Y(*,1) VALUES FOR NUMERICAL DIFFERENTIATION IN PJAC AND/OR PDF. |
|
|
57 IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. IDF .EQ. 0) GO TO 10 |
|
|
58 GO TO 15 |
|
|
59 10 CALL F (NEQ, TN, Y, PAR, SAVF) |
|
|
60 NFE = NFE + 1 |
|
|
61 C IF JCUR = 0, UPDATE THE JACOBIAN MATRIX. |
|
|
62 C IF MITER = 5, LOAD CORRECTED Y(*,1) VALUES INTO Y(*,2). |
|
|
63 15 IF (JCUR .EQ. 1) GO TO 30 |
|
|
64 IF (MITER .NE. 5) GO TO 25 |
|
|
65 DO 20 I = 1,N |
|
|
66 20 Y(I,2) = Y(I,1) |
|
|
67 25 CALL PJAC (NEQ, Y, Y(1,2), N, WM, IWM, EWT, SAVF, ACOR(1,2), |
|
|
68 1 PAR, F, JAC, JOPT) |
|
|
69 IF (IERPJ .NE. 0) RETURN |
|
|
70 C----------------------------------------------------------------------- |
|
|
71 C THIS IS A LOOPING POINT FOR THE SENSITIVITY CALCULATIONS. |
|
|
72 C----------------------------------------------------------------------- |
|
|
73 C FOR EACH PARAMETER PAR(*), A SENSITIVITY SOLUTION VECTOR IS COMPUTED |
|
|
74 C USING THE SAME STEP SIZE (H) AND ORDER (NQ) AS IN ODESSA_STODE. |
|
|
75 C A LOCAL ERROR TEST IS APPLIED INDEPENDENTLY TO EACH SOLUTION VECTOR. |
|
|
76 C----------------------------------------------------------------------- |
|
|
77 30 DO 100 J = 2,NSV |
|
|
78 JPAR = J - 1 |
|
|
79 C EVALUATE INHOMOGENEITY TERM, TEMPORARILY LOAD INTO Y(*,JPAR+1). ------ |
|
|
80 CALL PDF(NEQ, Y, WM, SAVF, ACOR(1,J), Y(1,J), PAR, |
|
|
81 1 F, DF, JPAR) |
|
|
82 C----------------------------------------------------------------------- |
|
|
83 C LOAD RHS OF SENSITIVITY SOLUTION (CORRECTOR) EQUATION.. |
|
|
84 C |
|
|
85 C RHS = DY/DP - EL(1)*H*D(DY/DP)/DT + EL(1)*H*DF/DP |
|
|
86 C |
|
|
87 C----------------------------------------------------------------------- |
|
|
88 DO 40 I = 1,N |
|
|
89 40 Y(I,J) = YH(I,J,1) - EL0*YH(I,J,2) + HL0*Y(I,J) |
|
|
90 C----------------------------------------------------------------------- |
|
|
91 C SOLVE CORRECTOR EQUATION: THE SOLUTIONS ARE LOCATED IN Y(*,JPAR+1). |
|
|
92 C THE EXPLICIT FORMULA IS.. |
|
|
93 C |
|
|
94 C (I - EL(1)*H*JAC) * DY/DP(CORRECTED) = RHS |
|
|
95 C |
|
|
96 C----------------------------------------------------------------------- |
|
|
97 CALL SOLVE (WM, IWM, Y(1,J), DUM) |
|
|
98 IF (IERSL .NE. 0) RETURN |
|
|
99 C ESTIMATE LOCAL TRUNCATION ERROR. ------------------------------------- |
|
|
100 DO 50 I = 1,N |
|
|
101 50 ACOR(I,J) = (Y(I,J) - YH(I,J,1))*EL0I |
|
|
102 ERR = ODESSA_VNORM(N, ACOR(1,J), EWT(1,J))*TI2 |
|
|
103 IF (ERR .GT. ONE) GO TO 200 |
|
|
104 C----------------------------------------------------------------------- |
|
|
105 C LOCAL ERROR TEST PASSED. SET KFLAGS TO 0 TO INDICATE THIS. |
|
|
106 C IF IALTH = 1, COMPUTE DSMS, DDNS, AND DUPS (IF L .LT. LMAX). |
|
|
107 C----------------------------------------------------------------------- |
|
|
108 KFLAGS = 0 |
|
|
109 IF (IALTH .GT. 1) GO TO 100 |
|
|
110 IF (L .EQ. LMAX) GO TO 70 |
|
|
111 DO 60 I= 1,N |
|
|
112 60 Y(I,J) = ACOR(I,J) - YH(I,J,LMAX) |
|
3987
|
113 DUPS = DMAX1(DUPS,ODESSA_VNORM(N,Y(1,J),EWT(1,J))*TI3) |
|
|
114 70 DSMS = DMAX1(DSMS,ERR) |
|
3983
|
115 100 CONTINUE |
|
|
116 RETURN |
|
|
117 C----------------------------------------------------------------------- |
|
|
118 C THIS SECTION IS REACHED IF THE ERROR TOLERANCE FOR SENSITIVITY |
|
|
119 C SOLUTION VECTOR JPAR HAS BEEN VIOLATED. KFLAGS IS MADE NEGATIVE TO |
|
|
120 C INDICATE THIS. IF KFLAGS = -1, SET KFLAG EQUAL TO ZERO SO THAT KFLAG |
|
|
121 C IS SET TO -1 ON RETURN TO ODESSA_STODE BEFORE REPEATING THE STEP. |
|
|
122 C INCREMENT NRS(1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO ALL |
|
|
123 C SENSITIVITY SOLUTION VECTORS) BY ONE. |
|
|
124 C INCREMENT NRS(JPAR+1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO |
|
|
125 C SOLUTION VECTOR JPAR+1) BY ONE. |
|
|
126 C LOAD DSMS FOR RH CALCULATION IN ODESSA_STODE. |
|
|
127 C----------------------------------------------------------------------- |
|
|
128 200 KFLAGS = KFLAGS - 1 |
|
|
129 IF (KFLAGS .EQ. -1) KFLAG = 0 |
|
|
130 NRS(1) = NRS(1) + 1 |
|
|
131 NRS(J) = NRS(J) + 1 |
|
|
132 DSMS = ERR |
|
|
133 RETURN |
|
|
134 C-------------------- END OF SUBROUTINE ODESSA_STESA ---------------------- |
|
|
135 END |
