Initial commit from Polytechnique Montreal

1 files changed, 114 insertions, 0 deletions

diff --git a/Dragon/src/EVOODE.f b/Dragon/src/EVOODE.f
new file mode 100644
index 0000000..dbdf0de
--- /dev/null
+++ b/Dragon/src/EVOODE.f
@@ -0,0 +1,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