1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
*DECK EVOODE
SUBROUTINE EVOODE(YSTART,NVAR,X1,X2,EPS,H1,NOK,NBAD,ITYPE,MU1,
1 IMA,MAXA,NSUPF,NFISS,KFISS,YSF,ADPL,BDPL)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Runge-Kutta or Kaps-Rentrop driver with adaptive stepsize control
* special version for isotopic depletion calculations.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input/output
* YSTART dependent variable vector.
* NVAR dimension of the dependent variable vector (number of
* depleting isotopes).
* X1 initial value of the independent variable.
* X2 final value of the independent variable.
* EPS required accuracy.
* H1 guessed first stepsize.
* NOK number of good steps taken.
* NBAD number of bad steps taken.
* ITYPE type of ODE solution:
* =1 fifth-order Runge-Kutta method;
* =2 fourth-order Kaps-Rentrop method.
* MU1 position of each diagonal element in matrix ADPL.
* IMA position of the first non-zero column element in matrix ADPL.
* MAXA first dimension of matrix ADPL.
* NSUPF number of depleting fission products.
* NFISS number of fissile isotopes producing fission products.
* KFISS position in chain of the fissile isotopes.
* YSF components of the product of the fission yields and fission
* rates.
* ADPL depletion matrix components.
* BDPL depletion source components.
*
*-----------------------------------------------------------------------
*
* REFERENCE:
* W.H. PRESS, B.P. FLANNERY, S.A. TEUKOLSKY AND W.T. VETTERLING,
* 'NUMERICAL RECIPIES (FORTRAN VERSION)', CAMBRIDGE UNIVERSITY PRESS,
* CHAPTER 15, CAMBRIDGE (1990).
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NVAR,NOK,NBAD,ITYPE,MU1(NVAR),IMA(NVAR),MAXA,NSUPF,NFISS,
1 KFISS(NFISS)
REAL YSTART(NVAR),X1,X2,EPS,H1,YSF(NFISS,NSUPF,2),ADPL(MAXA,2),
1 BDPL(NVAR,2)
*----
* LOCAL VARIABLES
*----
PARAMETER (MAXSTP=50000,TINY=1.E-8)
REAL, ALLOCATABLE, DIMENSION(:) :: YSCAL,Y
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(YSCAL(NVAR),Y(NVAR))
*
NSUPL=NVAR-NSUPF
X=X1
H=SIGN(H1,X2-X1)
NOK=0
NBAD=0
DO 10 I=1,NVAR
Y(I)=YSTART(I)
10 CONTINUE
DO 50 NSTP=1,MAXSTP
IF((X+H-X2)*(X+H-X1).GT.0.0) H=X2-X
CALL ALLUM(NVAR,ADPL(1,1),Y(1),YSCAL(1),MU1,IMA,1)
CALL ALLUM(NVAR,ADPL(1,2),Y(1),YSTART(1),MU1,IMA,1)
DO 25 I=1,NSUPF
DO 20 J=1,NFISS
YSCAL(NSUPL+I)=YSCAL(NSUPL+I)+YSF(J,I,1)*Y(KFISS(J))
YSTART(NSUPL+I)=YSTART(NSUPL+I)+YSF(J,I,2)*Y(KFISS(J))
20 CONTINUE
25 CONTINUE
DO 30 I=1,NVAR
YSCAL(I)=YSCAL(I)+BDPL(I,1)+X*(YSTART(I)+BDPL(I,2))
YSCAL(I)=MAX(ABS(Y(I))+ABS(H*YSCAL(I)),TINY)
30 CONTINUE
IF(ITYPE.EQ.1) THEN
CALL EVORK(Y,NVAR,X,H,EPS,YSCAL,HDID,HNEXT,MU1,IMA,MAXA,
1 NSUPF,NFISS,KFISS,YSF,ADPL,BDPL)
ELSE IF(ITYPE.EQ.2) THEN
CALL EVOKAP(Y,NVAR,X,H,EPS,YSCAL,HDID,HNEXT,MU1,IMA,MAXA,
1 NSUPF,NFISS,KFISS,YSF,ADPL,BDPL)
ENDIF
IF(HDID.EQ.H) THEN
NOK=NOK+1
ELSE
NBAD=NBAD+1
ENDIF
IF((X-X2)*(X2-X1).GE.0.0) THEN
DO 40 I=1,NVAR
YSTART(I)=Y(I)
40 CONTINUE
GO TO 60
ENDIF
H=HNEXT
50 CONTINUE
CALL XABORT('EVOODE: TOO MANY STEPS.')
*----
* SCRATCH STORAGE DEALLOCATION
*----
60 DEALLOCATE(Y,YSCAL)
END
|
