summaryrefslogtreecommitdiff
path: root/Dragon/src/EVORK.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/EVORK.f')
-rw-r--r--Dragon/src/EVORK.f193
1 files changed, 193 insertions, 0 deletions
diff --git a/Dragon/src/EVORK.f b/Dragon/src/EVORK.f
new file mode 100644
index 0000000..93758ee
--- /dev/null
+++ b/Dragon/src/EVORK.f
@@ -0,0 +1,193 @@
+*DECK EVORK
+ SUBROUTINE EVORK(Y,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,MU1,IMA,MAXA,
+ 1 NSUPF,NFISS,KFISS,YSF,ADPL,BDPL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Fifth-order Runge-Kutta Cash-Karp step with monitoring of local
+* truncation error to ensure accuracy and adjust stepsize.
+* 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
+* Y dependent variable vector.
+* N size of the dependent variable vector.
+* X independent variable.
+* HTRY stepsize to be attempted.
+* EPS required accuracy.
+* YSCAL vector against which the error is scaled.
+* HDID stepsize that was actually accomplished.
+* HNEXT estimated next stepsize.
+* MU1 position of each diagonal element in vectors ADPL and ASS.
+* IMA position of the first non-zero column element in vectors
+* ADPL and ASS.
+* 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, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
+* "Numerical recipes in Fortran, Second edition, Chapter 16,
+* Cambridge, 1992.
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER N,MU1(N),IMA(N),MAXA,NSUPF,NFISS,KFISS(NFISS)
+ REAL Y(N),X,HTRY,EPS,YSCAL(N),HDID,HNEXT,YSF(NFISS,NSUPF,2),
+ 1 ADPL(MAXA,2),BDPL(N,2)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (A2=.2,A3=.3,A4=.6,A5=1.,A6=.875,B21=.2,B31=3./40.,
+ *B32=9./40.,B41=.3,B42=-.9,B43=1.2,B51=-11./54.,B52=2.5,
+ *B53=-70./27.,B54=35./27.,B61=1631./55296.,B62=175./512.,
+ *B63=575./13824.,B64=44275./110592.,B65=253./4096.,C1=37./378.,
+ *C3=250./621.,C4=125./594.,C6=512./1771.,DC1=C1-2825./27648.,
+ *DC3=C3-18575./48384.,DC4=C4-13525./55296.,DC5=-277./14336.,
+ *DC6=C6-.25)
+ PARAMETER (SAFETY=0.85,PGROW=-.2,PSHRNK=-.25,GROW=1.5,SHRNK=0.5)
+ CHARACTER HSMG*131
+ REAL, ALLOCATABLE, DIMENSION(:) :: YTEMP,YGAR
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: AK
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(YTEMP(N),YGAR(N),AK(N,6))
+*
+ NSUPL=N-NSUPF
+ H=HTRY
+ 10 CALL ALLUM(N,ADPL(1,1),Y(1),AK(1,1),MU1,IMA,1)
+ CALL ALLUM(N,ADPL(1,2),Y(1),YGAR(1),MU1,IMA,1)
+ DO 25 I=1,NSUPF
+ DO 20 J=1,NFISS
+ AK(NSUPL+I,1)=AK(NSUPL+I,1)+YSF(J,I,1)*Y(KFISS(J))
+ YGAR(NSUPL+I)=YGAR(NSUPL+I)+YSF(J,I,2)*Y(KFISS(J))
+ 20 CONTINUE
+ 25 CONTINUE
+ DO 30 I=1,N
+ AK(I,1)=AK(I,1)+BDPL(I,1)+X*(YGAR(I)+BDPL(I,2))
+ YTEMP(I)=Y(I)+H*B21*AK(I,1)
+ 30 CONTINUE
+*
+ CALL ALLUM(N,ADPL(1,1),YTEMP(1),AK(1,2),MU1,IMA,1)
+ CALL ALLUM(N,ADPL(1,2),YTEMP(1),YGAR(1),MU1,IMA,1)
+ DO 45 I=1,NSUPF
+ DO 40 J=1,NFISS
+ AK(NSUPL+I,2)=AK(NSUPL+I,2)+YSF(J,I,1)*YTEMP(KFISS(J))
+ YGAR(NSUPL+I)=YGAR(NSUPL+I)+YSF(J,I,2)*YTEMP(KFISS(J))
+ 40 CONTINUE
+ 45 CONTINUE
+ DO 50 I=1,N
+ AK(I,2)=AK(I,2)+BDPL(I,1)+(X+A2*H)*(YGAR(I)+BDPL(I,2))
+ YTEMP(I)=Y(I)+H*(B31*AK(I,1)+B32*AK(I,2))
+ 50 CONTINUE
+*
+ CALL ALLUM(N,ADPL(1,1),YTEMP(1),AK(1,3),MU1,IMA,1)
+ CALL ALLUM(N,ADPL(1,2),YTEMP(1),YGAR(1),MU1,IMA,1)
+ DO 65 I=1,NSUPF
+ DO 60 J=1,NFISS
+ AK(NSUPL+I,3)=AK(NSUPL+I,3)+YSF(J,I,1)*YTEMP(KFISS(J))
+ YGAR(NSUPL+I)=YGAR(NSUPL+I)+YSF(J,I,2)*YTEMP(KFISS(J))
+ 60 CONTINUE
+ 65 CONTINUE
+ DO 70 I=1,N
+ AK(I,3)=AK(I,3)+BDPL(I,1)+(X+A3*H)*(YGAR(I)+BDPL(I,2))
+ YTEMP(I)=Y(I)+H*(B41*AK(I,1)+B42*AK(I,2)+B43*AK(I,3))
+ 70 CONTINUE
+*
+ CALL ALLUM(N,ADPL(1,1),YTEMP(1),AK(1,4),MU1,IMA,1)
+ CALL ALLUM(N,ADPL(1,2),YTEMP(1),YGAR(1),MU1,IMA,1)
+ DO 85 I=1,NSUPF
+ DO 80 J=1,NFISS
+ AK(NSUPL+I,4)=AK(NSUPL+I,4)+YSF(J,I,1)*YTEMP(KFISS(J))
+ YGAR(NSUPL+I)=YGAR(NSUPL+I)+YSF(J,I,2)*YTEMP(KFISS(J))
+ 80 CONTINUE
+ 85 CONTINUE
+ DO 90 I=1,N
+ AK(I,4)=AK(I,4)+BDPL(I,1)+(X+A4*H)*(YGAR(I)+BDPL(I,2))
+ YTEMP(I)=Y(I)+H*(B51*AK(I,1)+B52*AK(I,2)+B53*AK(I,3)+B54*AK(I,4))
+ 90 CONTINUE
+*
+ CALL ALLUM(N,ADPL(1,1),YTEMP(1),AK(1,5),MU1,IMA,1)
+ CALL ALLUM(N,ADPL(1,2),YTEMP(1),YGAR(1),MU1,IMA,1)
+ DO 105 I=1,NSUPF
+ DO 100 J=1,NFISS
+ AK(NSUPL+I,5)=AK(NSUPL+I,5)+YSF(J,I,1)*YTEMP(KFISS(J))
+ YGAR(NSUPL+I)=YGAR(NSUPL+I)+YSF(J,I,2)*YTEMP(KFISS(J))
+ 100 CONTINUE
+ 105 CONTINUE
+ DO 110 I=1,N
+ AK(I,5)=AK(I,5)+BDPL(I,1)+(X+A5*H)*(YGAR(I)+BDPL(I,2))
+ YTEMP(I)=Y(I)+H*(B61*AK(I,1)+B62*AK(I,2)+B63*AK(I,3)+B64*AK(I,4)+
+ 1 B65*AK(I,5))
+ 110 CONTINUE
+*
+ CALL ALLUM(N,ADPL(1,1),YTEMP(1),AK(1,6),MU1,IMA,1)
+ CALL ALLUM(N,ADPL(1,2),YTEMP(1),YGAR(1),MU1,IMA,1)
+ DO 125 I=1,NSUPF
+ DO 120 J=1,NFISS
+ AK(NSUPL+I,6)=AK(NSUPL+I,6)+YSF(J,I,1)*YTEMP(KFISS(J))
+ YGAR(NSUPL+I)=YGAR(NSUPL+I)+YSF(J,I,2)*YTEMP(KFISS(J))
+ 120 CONTINUE
+ 125 CONTINUE
+ DO 130 I=1,N
+ AK(I,6)=AK(I,6)+BDPL(I,1)+(X+A6*H)*(YGAR(I)+BDPL(I,2))
+ YTEMP(I)=Y(I)+H*(C1*AK(I,1)+C3*AK(I,3)+C4*AK(I,4)+C6*AK(I,6))
+ YGAR(I)=H*(DC1*AK(I,1)+DC3*AK(I,3)+DC4*AK(I,4)+DC5*AK(I,5)+
+ 1 DC6*AK(I,6))
+ 130 CONTINUE
+*
+ ERRMAX=0.0
+ DO 140 I=1,N
+ ERRMAX=MAX(ERRMAX,ABS(YGAR(I)/YSCAL(I)))
+ 140 CONTINUE
+ ERRMAX=ERRMAX/EPS
+ IF (ERRMAX.EQ.0.0) THEN
+ HDID=H
+ HNEXT=GROW*H
+ X=X+H
+ DO 150 I=1,N
+ Y(I)=YTEMP(I)
+ 150 CONTINUE
+ GO TO 170
+ ELSE IF (ERRMAX.LE.1.0) THEN
+ HDID=H
+ HNEXT=MIN(GROW,SAFETY*(ERRMAX**PGROW))*H
+ X=X+H
+ DO 160 I=1,N
+ Y(I)=YTEMP(I)
+ 160 CONTINUE
+ GO TO 170
+ ELSE
+ H=MAX(SHRNK,SAFETY*(ERRMAX**PSHRNK))*H
+ XNEW=X+H
+ IF (X.EQ.XNEW) THEN
+ WRITE(HSMG,'(35HEVORK: STEPSIZE NOT SIGNIFICANT (H=,1P,E11.4,
+ 1 6H HTRY=,E11.4,2H).)') H,HTRY
+ CALL XABORT(HSMG)
+ ENDIF
+ GO TO 10
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ 170 DEALLOCATE(AK,YGAR,YTEMP)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:33:22 +0000