From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Donjon/src/AFMTAV.f | 173 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 Donjon/src/AFMTAV.f (limited to 'Donjon/src/AFMTAV.f') diff --git a/Donjon/src/AFMTAV.f b/Donjon/src/AFMTAV.f new file mode 100644 index 0000000..481280f --- /dev/null +++ b/Donjon/src/AFMTAV.f @@ -0,0 +1,173 @@ +*DECK AFMTAV + SUBROUTINE AFMTAV (NBURN,ITM,XBMAX,XBMIN,YS,NBMIN,NBMAX,XB,SIGAV) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Time average calculation using different approximation. +* +*Copyright: +* Copyright (C) 1996 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): +* M.T. Sissaoui +* +*Parameters: input +* NBURN total number of steps. +* ITM type of the approximation (1-Lagrange; 2-spline; 3-Hermite) +* XBMAX highest value. +* XBMIN lower value. +* YS parameter to be integrated +* NBMIN +* NBMAX +* XB steps +* +*Parameters: output +* SIGAV average value of YS +* +*----------------------------------------------------------------------- +* + REAL YS(NBURN),XB(NBURN),SIGAV + REAL UU(2) + DOUBLE PRECISION DD + REAL, ALLOCATABLE, DIMENSION(:) :: Y,U +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(Y(NBURN),U(NBURN)) +* + IF(NBMAX.GT.0) THEN + INMAX=NBMAX-2 + INMIN=NBMIN + ELSE + INMAX=ABS(NBMAX)-1 + INMIN=NBMIN + IF(ITM.EQ.1) + 1 CALL XABORT('AFMTAV: MORE BURNUP STEPS ARE REQUIRED TO USE ' + 1 //' LAGRANGE METHOD, CHOOSE HERMIT OR SPLINE METHOD') + ENDIF + IF(ABS(NBMAX).GT.NBMIN) THEN + SIGAV=0.0 +* + IF(ITM.EQ.1) THEN +* TIME AVERAGE CALCULATION USING LAGRANGE APPROXIMATION. +* + DO 113 IR=INMIN,INMAX + I1=IR + I2=IR+1 + I3=IR+2 + XBI=MAX(XBMIN,XB(IR)) + XBF=MIN(XBMAX,XB(IR+1)) + TX=XBF-XBI + TX2=XBF**2-XBI**2 + TX3=XBF**3-XBI**3 + X12=XB(I1)-XB(I2) + X13=XB(I1)-XB(I3) + X23=XB(I2)-XB(I3) + XA12=XB(I1)+XB(I2) + XA13=XB(I1)+XB(I3) + XA23=XB(I2)+XB(I3) + XM12=XB(I1)*XB(I2) + XM13=XB(I1)*XB(I3) + XM23=XB(I2)*XB(I3) + Y1=YS(I1)/(X12*X13) + Y2=-YS(I2)/(X12*X23) + Y3=YS(I3)/(X13*X23) +* + SIGAV=SIGAV + + 1 Y1*(TX3/3.0-XA23*TX2/2.0+XM23*TX)+ + 1 Y2*(TX3/3.0-XA13*TX2/2.0+XM13*TX)+ + 1 Y3*(TX3/3.0-XA12*TX2/2.0+XM12*TX) + 113 CONTINUE +* + ELSE IF(ITM.EQ.2) THEN +* TIME AVERAGE CALCULATION USING SPLINE APPROXIMATION. +* THE LOWER BOUNDARY CONDITION IS SET TO BE NATURAL + Y(1)=0.0 + U(1)=0.0 +* THE UPPER BOUNDARY CONDITION IS SET EITHER TO BE NATURAL + QN=0.0 + UN=0.0 +* + DO 103 IR=2,NBURN-1 + SIG=(XB(IR)-XB(IR-1))/(XB(IR+1)-XB(IR-1)) + P=SIG*Y(IR-1)+2.0 + Y(IR)=(SIG-1.0)/P + U(IR)=(6.*((YS(IR+1)-YS(IR))/(XB(IR+1)-XB(IR))- + 1 (YS(IR)-YS(IR-1))/(XB(IR)-XB(IR-1)))/(XB(IR+1)- + 1 XB(IR-1))-SIG*U(IR-1))/P + 103 CONTINUE +* + Y(NBURN)=(UN-QN*U(NBURN-1))/(QN*Y(NBURN-1)+1.0) +* + DO 104 K=NBURN-1,1,-1 + Y(K)=Y(K)*Y(K+1)+U(K) + 104 CONTINUE +* +* COMPUTE THE INTEGRAL OF THE X-SECTION + INMAX=NBMAX-2 + INMIN=NBMIN + DO 300 IR=INMIN,INMAX + H=XB(IR+1)-XB(IR) + XBI=MAX(XBMIN,XB(IR)) + XBF=MIN(XBMAX,XB(IR+1)) +* + DB=XBF-XBI + HF=XB(IR+1)-XBF + HI=XB(IR+1)-XBI +* + AI=-0.5*(HF**2-HI**2)/H + BI=DB-AI + CI=-(AI/6)*H**2-(HF**4-HI**4)/(24*H) + DI=-(BI/6)*H**2-(HF**4-HI**4)/(24*H) +* + SIGAV=SIGAV+AI*YS(IR)+BI*YS(IR+1)+ + 1 CI*Y(IR)+DI*Y(IR+1) + 300 CONTINUE + ELSE IF(ITM.EQ.3) THEN +* TIME AVERAGE CALCULATION USING HERMIT APPROXIMATION. + DO 101 I=1,NBURN + Y(I)=YS(I) + 101 CONTINUE +* TAKE THE DERIVATIVE WITH RESPECT TO BURNUP OR NEUTRON EXPOSURE AT +* TABULATION POINTS. + CALL ALDERV(NBURN,XB,Y) +* +* COMPUTE THE INTEGRAL OF THE X-SECTION + DD=0.0D0 + DO 200 IR=1,NBURN-1 + IF((XBMIN.LT.XB(IR+1)).AND.(XBMAX.GT.XB(IR))) THEN + DX=XB(IR+1)-XB(IR) + XBI=MAX(XBMIN,XB(IR)) + XBF=MIN(XBMAX,XB(IR+1)) + CC=0.5*(XBF-XBI) + U1=(XBI-0.5*(XB(IR)+XB(IR+1)))/DX + U2=(XBF-0.5*(XB(IR)+XB(IR+1)))/DX + UU(1)=0.5*(-(U2-U1)*0.577350269189626+U1+U2) + UU(2)=0.5*((U2-U1)*0.577350269189626+U1+U2) + DO 190 J=1,2 + H1=3.0*(0.5-UU(J))**2-2.0*(0.5-UU(J))**3 + H2=(0.5-UU(J))**2-(0.5-UU(J))**3 + H3=3.0*(0.5+UU(J))**2-2.0*(0.5+UU(J))**3 + H4=-(0.5+UU(J))**2+(0.5+UU(J))**3 + DD=DD+(H1*YS(IR)+H2*Y(IR)*DX+H3*YS(IR+1)+ + 1 H4*Y(IR+1)*DX)*CC + 190 CONTINUE + ENDIF + 200 CONTINUE + SIGAV=REAL(DD) + ENDIF + SIGAV=SIGAV/(XBMAX-XBMIN) + ELSE + SIGAV=YS(NBMIN) + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(U,Y) + RETURN + END -- cgit v1.2.3