Initial commit from Polytechnique Montreal

1 files changed, 224 insertions, 0 deletions

diff --git a/Dragon/src/EVOKAP.f b/Dragon/src/EVOKAP.f
new file mode 100644
index 0000000..a280232
--- /dev/null
+++ b/Dragon/src/EVOKAP.f
@@ -0,0 +1,224 @@
+*DECK EVOKAP
+      SUBROUTINE EVOKAP(Y,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,MU1,IMA,MAXA,
+     1 NSUPF,NFISS,KFISS,YSF,ADPL,BDPL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Fourth-order Kaps-Rentrop step for integrating stiff O.D.E.'s, with
+* monitoring of local truncation error to 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 and S.A. Teukolsky, 'Integrating stiff ordinary differen-
+* tial equations', Computers in physics, 3 (3), 88 (May/June 1989).
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  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 (MAXTRY=40,SAFETY=0.85,GROW=1.5,PGROW=-.25,SHRNK=0.5,
+     1 PSHRNK=-1./3.)
+      PARAMETER (GAM=.231,GAM21=-.270629667752/GAM,
+     1 GAM31=.311254483294/GAM,GAM32=.852445628482E-2/GAM,
+     2 GAM41=.282816832044/GAM,GAM42=-.457959483281/GAM,
+     3 GAM43=-.111208333333/GAM,ALF21=.462,ALF31=-.815668168327E-1,
+     4 ALF32=.961775150166,C1=.217487371653,C2=.486229037990,C3=0.,
+     5 C4=.296283590357,CC1=-.717088504499,CC2=1.77617912176,
+     6 CC3=-.590906172617E-1,GAM2X=GAM*(1.+GAM21),
+     7 GAM3X=GAM*(1.+GAM31+GAM32),GAM4X=GAM*(1.+GAM41+GAM42+GAM43))
+      CHARACTER HSMG*131
+      REAL, ALLOCATABLE, DIMENSION(:) :: DYDX,TEMP,YSAV,DYSAV,DFDX,ASS
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: AK
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(DYDX(N),TEMP(N),YSAV(N),DYSAV(N),DFDX(N),AK(N,4),
+     1 ASS(IMA(N)))
+*
+      XSAV=X
+      NSUPL=N-NSUPF
+      CALL ALLUM(N,ADPL(1,1),Y(1),DYDX(1),MU1,IMA,1)
+      CALL ALLUM(N,ADPL(1,2),Y(1),TEMP(1),MU1,IMA,1)
+      DO 15 I=1,NSUPF
+      DO 10 J=1,NFISS
+      DYDX(NSUPL+I)=DYDX(NSUPL+I)+YSF(J,I,1)*Y(KFISS(J))
+      TEMP(NSUPL+I)=TEMP(NSUPL+I)+YSF(J,I,2)*Y(KFISS(J))
+   10 CONTINUE
+   15 CONTINUE
+      DO 20 I=1,N
+      DYDX(I)=DYDX(I)+X*TEMP(I)+BDPL(I,1)+X*BDPL(I,2)
+      YSAV(I)=Y(I)
+      DYSAV(I)=DYDX(I)
+   20 CONTINUE
+      H=HTRY
+      DO 200 JTRY=1,MAXTRY
+      HSQ=H*H
+      CALL ALLUM(N,ADPL(1,2),YSAV(1),DFDX(1),MU1,IMA,1)
+      DO 35 I=1,NSUPF
+      DO 30 J=1,NFISS
+      DFDX(NSUPL+I)=DFDX(NSUPL+I)+YSF(J,I,2)*YSAV(KFISS(J))
+   30 CONTINUE
+   35 CONTINUE
+      DO 40 I=1,IMA(N)
+      ASS(I)=-H*GAM*(ADPL(I,1)+XSAV*ADPL(I,2))
+   40 CONTINUE
+      DO 50 I=1,N
+      DFDX(I)=DFDX(I)+BDPL(I,2)
+      ASS(MU1(I))=1.+ASS(MU1(I))
+   50 CONTINUE
+      CALL ALLUF(N,ASS,MU1,IMA)
+      DO 60 I=1,N
+      AK(I,1)=H*DYSAV(I)+HSQ*GAM*DFDX(I)
+   60 CONTINUE
+      CALL ALLUS(NSUPL,MU1(1),IMA(1),ASS(1),AK(1,1))
+      IF(NSUPF.GT.0) THEN
+        DO 75 I=1,NSUPF
+        DO 70 J=1,NFISS
+        AK(NSUPL+I,1)=AK(NSUPL+I,1)+H*GAM*(YSF(J,I,1)+XSAV*YSF(J,I,2))
+     1  *AK(KFISS(J),1)
+   70   CONTINUE
+   75   CONTINUE
+        CALL ALLUS(NSUPF,MU1(NSUPL+1),IMA(NSUPL+1),ASS(1),AK(NSUPL+1,1))
+      ENDIF
+      DO 80 I=1,N
+      Y(I)=YSAV(I)+ALF21*AK(I,1)
+   80 CONTINUE
+      X=XSAV+ALF21*H
+      CALL ALLUM(N,ADPL(1,1),Y(1),DYDX(1),MU1,IMA,1)
+      CALL ALLUM(N,ADPL(1,2),Y(1),TEMP(1),MU1,IMA,1)
+      DO 95 I=1,NSUPF
+      DO 90 J=1,NFISS
+      DYDX(NSUPL+I)=DYDX(NSUPL+I)+YSF(J,I,1)*Y(KFISS(J))
+      TEMP(NSUPL+I)=TEMP(NSUPL+I)+YSF(J,I,2)*Y(KFISS(J))
+   90 CONTINUE
+   95 CONTINUE
+      DO 100 I=1,N
+      DYDX(I)=DYDX(I)+X*TEMP(I)+BDPL(I,1)+X*BDPL(I,2)
+      AK(I,2)=H*DYDX(I)+HSQ*GAM2X*DFDX(I)+GAM21*AK(I,1)
+  100 CONTINUE
+      CALL ALLUS(NSUPL,MU1(1),IMA(1),ASS(1),AK(1,2))
+      IF(NSUPF.GT.0) THEN
+        DO 106 I=1,NSUPF
+        DO 105 J=1,NFISS
+        AK(NSUPL+I,2)=AK(NSUPL+I,2)+H*GAM*(YSF(J,I,1)+XSAV*YSF(J,I,2))
+     1  *AK(KFISS(J),2)
+  105   CONTINUE
+  106   CONTINUE
+        CALL ALLUS(NSUPF,MU1(NSUPL+1),IMA(NSUPL+1),ASS(1),AK(NSUPL+1,2))
+      ENDIF
+      DO 110 I=1,N
+      AK(I,2)=AK(I,2)-GAM21*AK(I,1)
+      Y(I)=YSAV(I)+ALF31*AK(I,1)+ALF32*AK(I,2)
+  110 CONTINUE
+      X=XSAV+(ALF31+ALF32)*H
+      CALL ALLUM(N,ADPL(1,1),Y(1),DYDX(1),MU1,IMA,1)
+      CALL ALLUM(N,ADPL(1,2),Y(1),TEMP(1),MU1,IMA,1)
+      DO 125 I=1,NSUPF
+      DO 120 J=1,NFISS
+      DYDX(NSUPL+I)=DYDX(NSUPL+I)+YSF(J,I,1)*Y(KFISS(J))
+      TEMP(NSUPL+I)=TEMP(NSUPL+I)+YSF(J,I,2)*Y(KFISS(J))
+  120 CONTINUE
+  125 CONTINUE
+      DO 130 I=1,N
+      DYDX(I)=DYDX(I)+X*TEMP(I)+BDPL(I,1)+X*BDPL(I,2)
+      TEMP(I)=GAM31*AK(I,1)+GAM32*AK(I,2)
+      AK(I,3)=H*DYDX(I)+GAM3X*HSQ*DFDX(I)+TEMP(I)
+  130 CONTINUE
+      CALL ALLUS(NSUPL,MU1(1),IMA(1),ASS(1),AK(1,3))
+      IF(NSUPF.GT.0) THEN
+        DO 136 I=1,NSUPF
+        DO 135 J=1,NFISS
+        AK(NSUPL+I,3)=AK(NSUPL+I,3)+H*GAM*(YSF(J,I,1)+XSAV*YSF(J,I,2))
+     1  *AK(KFISS(J),3)
+  135   CONTINUE
+  136   CONTINUE
+        CALL ALLUS(NSUPF,MU1(NSUPL+1),IMA(NSUPL+1),ASS(1),AK(NSUPL+1,3))
+      ENDIF
+      DO 140 I=1,N
+      AK(I,3)=AK(I,3)-TEMP(I)
+      TEMP(I)=GAM41*AK(I,1)+GAM42*AK(I,2)+GAM43*AK(I,3)
+      AK(I,4)=H*DYDX(I)+HSQ*GAM4X*DFDX(I)+TEMP(I)
+  140 CONTINUE
+      CALL ALLUS(NSUPL,MU1(1),IMA(1),ASS(1),AK(1,4))
+      IF(NSUPF.GT.0) THEN
+        DO 146 I=1,NSUPF
+        DO 145 J=1,NFISS
+        AK(NSUPL+I,4)=AK(NSUPL+I,4)+H*GAM*(YSF(J,I,1)+XSAV*YSF(J,I,2))
+     1  *AK(KFISS(J),4)
+  145   CONTINUE
+  146   CONTINUE
+        CALL ALLUS(NSUPF,MU1(NSUPL+1),IMA(NSUPL+1),ASS(1),AK(NSUPL+1,4))
+      ENDIF
+      DO 150 I=1,N
+      AK(I,4)=AK(I,4)-TEMP(I)
+      Y(I)=YSAV(I)+C1*AK(I,1)+C2*AK(I,2)+C3*AK(I,3)+C4*AK(I,4)
+      TEMP(I)=YSAV(I)+CC1*AK(I,1)+CC2*AK(I,2)+CC3*AK(I,3)
+  150 CONTINUE
+      X=XSAV+H
+      IF (X.EQ.XSAV) THEN
+         WRITE(HSMG,'(36HEVOKAP: STEPSIZE NOT SIGNIFICANT (H=,1P,E11.4,
+     1   6H HTRY=,E11.4,2H).)') H,HTRY
+         CALL XABORT(HSMG)
+      ENDIF
+      ERRMAX=0.
+      DO 160 I=1,N
+      ERRMAX=MAX(ERRMAX,ABS((Y(I)-TEMP(I))/YSCAL(I)))
+  160 CONTINUE
+      ERRMAX=ERRMAX/EPS
+      IF (ERRMAX.EQ.0.) THEN
+         HDID=H
+         HNEXT=GROW*H
+         GO TO 210
+      ELSE IF (ERRMAX.LE.1.) THEN
+         HDID=H
+         HNEXT=MIN(GROW,SAFETY*(ERRMAX**PGROW))*H
+         GO TO 210
+      ELSE
+         H=MAX(SHRNK,SAFETY*(ERRMAX**PSHRNK))*H
+      ENDIF
+  200 CONTINUE
+      CALL XABORT('EVOKAP: EXCEEDED MAXTRY.')
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+  210 DEALLOCATE(ASS,AK,DFDX,DYSAV,YSAV,TEMP,DYDX)
+      RETURN
+      END