diff options
Diffstat (limited to 'Donjon/src/AFMCPT.f')
| -rw-r--r-- | Donjon/src/AFMCPT.f | 467 |
1 files changed, 467 insertions, 0 deletions
diff --git a/Donjon/src/AFMCPT.f b/Donjon/src/AFMCPT.f new file mode 100644 index 0000000..5b697a7 --- /dev/null +++ b/Donjon/src/AFMCPT.f @@ -0,0 +1,467 @@ +*DECK AFMCPT + SUBROUTINE AFMCPT (KENTRY,NBURN,NGRP,NISO,NL,IMPX, + 1 SMACB,XBORB,XPURB,XXENB,XT1FB,XT2FB,XT1CB, + 1 XT2CB,XT1MB,XT2MB,XD1CB,XD2CB,XD1MB,XD2MB, + 1 XSMB,XNP9B,XMFDB,XMMDB,XPF1B,XPF2B,XPF1LB,XPF2LB, + 1 DENSITB,CPW1B,CPW2B,FLUXB,OVERVB,CHIB, + 1 IJJ,NJJ,HISO,CTITRE,NMIX,SIGMA,NTYP,TF,TC, + 1 TM,DC,DM,BOR,XEN,SM,RNP9,XI,TFR,TCR,TMR,XIR, + 1 OVERV,FLUX,CHI,SCAT,MMIX,NPS,PW,XBRH,XBURN, + 1 LTAV,IRAV,IDF,JTAB,IXYZ,ILIN) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the cross sections +* +*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 +* KENTRY address of the linked list or handle to the xsm file. +* NGRP number of energy groups. +* NISO number of extracted isotopes. +* NL number of Legendre orders (=1 for isotropic scattering). +* NBURN number of burnup steps. +* NTYP total number of cross sections type. +* TF fuel temperature. +* TC coolant temperature. +* TM moderator temperature. +* DC coolant density. +* DM moderator density. +* BOR Boron concentration. +* XEN Xenon concentration. +* SM Samarium concentration. +* RNP9 Neptunium concentration. +* XI moderator purity. +* NMIX mixture number. +* NPS NPS-2 power shift. +* IXYZ type of diffusion coefficient (=0: isotropic; =1: directional) +* +*Parameters: +* IMPX +* SMACB +* XBORB +* XPURB +* XXENB +* XT1FB +* XT2FB +* XT1CB +* XT2CB +* XT1MB +* XT2MB +* XD1CB +* XD2CB +* XD1MB +* XD2MB +* XSMB +* XNP9B +* XMFDB +* XMMDB +* XPF1B +* XPF2B +* XPF1LB +* XPF2LB +* DENSITB +* CPW1B +* CPW2B +* FLUXB +* OVERVB +* CHIB +* IJJ +* NJJ +* HISO +* CTITRE +* SIGMA +* TFR +* TCR +* TMR +* XIR +* OVERV +* FLUX +* CHI +* SCAT +* MMIX +* PW +* XBRH +* XBURN +* LTAV +* IRAV +* IDF +* JTAB +* ILIN +* +*----------------------------------------------------------------------- +* + CHARACTER HMICRO*12,CTITRE*72 + LOGICAL LTAV + DOUBLE PRECISION XCOF(3) + REAL CPF1(3) + DIMENSION KENTRY(*), + 1 SMACB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XBORB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XXENB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XT1FB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XT2FB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XT1CB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XT2CB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XT1MB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XT2MB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XD1CB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XD2CB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XD1MB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XD2MB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XSMB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XNP9B(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XMFDB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XMMDB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XPF1B(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XPF2B(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XPF1LB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XPF2LB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 XPURB(NGRP*NGRP,NTYP,NISO,NBURN,*), + 1 DENSITB(NISO,NBURN,*),CPW1B(2,NBURN,*), + 1 CPW2B(2,NBURN,*),FLUXB(NGRP,NBURN,*), + 1 CHIB(NGRP,NBURN,*),OVERVB(NGRP,NBURN,*), + 2 HISO(*),JTAB(*), + 6 ELMT(3),PW(*), + 7 SIGMA(MMIX,NGRP,NTYP),XBRH(*), + 8 OVERV(MMIX,*),FLUX(MMIX,*), + 7 CHI(MMIX,*),SCAT(MMIX,NL,NGRP,*), + 1 XBURN(NBURN,*),NJJ(*),IJJ(*) + REAL, ALLOCATABLE, DIMENSION(:) :: DEL + REAL, ALLOCATABLE, DIMENSION(:,:) :: CPW1,CPW2 + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: SMAC +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(SMAC(NGRP*NGRP,NTYP,NISO),DEL(NISO),CPW1(3,NPS), + 1 CPW2(3,NPS)) +* +* Remove warning for uninitialized variable + ISOB=0 + D2MB=0.0 + D2CB=0.0 + D1M=0.0 + D2C=0.0 + CX=0.0 + CSM=0.0 + D1C=0.0 + D2M=0.0 + CNP9=0.0 + CB=0.0 + PUR=0.0 +*---------------------------------------------------------------* + IF(NBURN.EQ.0) CALL XABORT('AFMCPT: ZERO NUMBER OF MIXTURES.') + IF(NGRP.EQ.0) CALL XABORT('AFMCPT: ZERO NUMBER OF GROUPS.') +*---------------------------------------------------------------* + NG2=NGRP*NGRP + NTM=4+2*IXYZ + IF(IMPX.GT.5) THEN + WRITE(6,'(9H AFMCPT: ,A)') CTITRE + WRITE(6,*) ' NTYP ',NTYP + ENDIF + XECON=0.0 + TFAV=0.0 + DCAV=0.0 + CPF2=0.0 + CF=0.0 + DO 81 ISO=1,NISO + DO 80 ITY=1,NTYP + DO 79 IGR=1,NG2 + SMAC(IGR,ITY,ISO)=0.0 + 79 CONTINUE + 80 CONTINUE + 81 CONTINUE + DO 82 IGR=1,NGRP + FLUX(NMIX,IGR)=0.0 + OVERV(NMIX,IGR)=0.0 + CHI(NMIX,IGR)=0.0 + 82 CONTINUE + DO 161 ISO=1,NISO + DEL(ISO)=0.0 + 161 CONTINUE + DO 73 II=1,3 + XCOF(II)=0.0D0 + ELMT(II)=0.0 + CPF1(II)=0.0 + 73 CONTINUE + DO 74 K=1,NPS + DO II=1,3 + CPW1(II,K)=0.0 + CPW2(II,K)=0.0 + ENDDO + 74 CONTINUE +* + PREF=PW(1) + IPFBM=KENTRY(2) +*---------------------------------------------------------------* +* COMPUTE THE HISTORY COEFFICIENTS + NPSX=NPS+NPS-1 + DO 101 K=2,NPSX +* XSECTION FOR SNAP-SHOT OR (TIME AVERAGE-HOMOG.) + IF(LTAV) THEN + IRMAX=IRAV + IRMIN=IRAV + XCOF(1)=1.0D0 + ELSE +* + IF(K.GT.NPS) THEN + XIRAD=XBRH(K-NPS+1) + ELSE + XIRAD=XBRH(NPS)-XBRH(K-1) + ENDIF + CALL AFMLOC(NBURN,NTP,XIRAD,XIRAD,XBURN(1,IDF), + 1 IRMAX,IRMIN,XCOF,ILIN) + ENDIF +* +* INTERPOLATE THE HISTORY COEFFICIENTS + IH=0 + DO 910 I = IRMIN,IRMAX + IH=IH+1 + IF(K.LE.NPS) THEN + IF(PW(K).GT.PREF) THEN + CPW1(IH,K)=CPW1B(1,I,IDF) + CPW1(IH,K) + CPW2(IH,K)=CPW2B(1,I,IDF) + CPW2(IH,K) + ELSE + CPW1(IH,K)=CPW1B(2,I,IDF) + CPW1(IH,K) + CPW2(IH,K)=CPW2B(2,I,IDF) + CPW2(IH,K) + ENDIF +* + IF(K.GT.2) THEN + IF(PW(K-1).GT.PREF) THEN + CPW1(IH,K-1)=-CPW1B(1,I,IDF) + CPW1(IH,K-1) + CPW2(IH,K-1)=-CPW2B(1,I,IDF) + CPW2(IH,K-1) + ELSE + CPW1(IH,K-1)=-CPW1B(2,I,IDF) + CPW1(IH,K-1) + CPW2(IH,K-1)=-CPW2B(2,I,IDF) + CPW2(IH,K-1) + ENDIF + ENDIF + ENDIF + 910 CONTINUE + 101 CONTINUE +* + YF=0.0 + DO 111 K=2,NPS + IF(K.EQ.2) THEN +* CORRECTE THE STURATING PSEUDO-FISSILE ISOTOPE + IF(PW(K).GT.PREF) THEN + XPW=ALOG(PW(K)/PW(1)) + XPWM=1.0/PW(K)-1.0/PW(1) + YF=0.0 + ELSE + XPW=PW(K)-PW(1) + XPWM=(PW(K)-PW(1))**2 + YF=1.0 + ENDIF + DO IH=1,3 + CPF1(IH)=CPF1(IH) +CPW1(IH,K)*XPW + CPW2(IH,K)*XPWM + ENDDO +C + ELSE + IF(PW(K).GT.PREF) THEN + XPW=ALOG(PW(K)/PW(1)) + XPWM=1.0/PW(K)-1.0/PW(1) + ELSE + XPW=PW(K)-PW(1) + XPWM=(PW(K)-PW(1))**2 +C + ENDIF + DO IH=1,3 + CPF1(IH)=CPF1(IH) +CPW1(IH,K)*XPW + CPW2(IH,K)*XPWM + ENDDO + ENDIF + 111 CONTINUE +*---------------------------------------------------------------* +* COMPUTE THE DEVIATION OF THE PSEUDO-ISOTOPE +* CPF1 AND CPF2 +* CPF2=CPF1*CPF1 +*---------------------------------------------------------------* +* APPLY THE FEEDBACK MODEL + T1F=SQRT(TF)-SQRT(TFR) + T2F=TF-TFR + T1C=ALOG(TC/TCR) + T2C=1.0/TC - 1.0/TCR + T1M=ALOG(TM/TMR) + T2M=1.0/TM - 1.0/TMR +* +* RECOVER LOCAL PARAMETER COEFFICIENT AND X-SECTIONS + II=0 + DMOD=0.0 + DO 900 I = IRMIN,IRMAX + II=II+1 + RXCOEF=REAL(XCOF(II)) + DO 249 IGR=1,NGRP + FLUX(NMIX,IGR)=FLUX(NMIX,IGR)+RXCOEF*FLUXB(IGR,I,IDF) + OVERV(NMIX,IGR)=OVERV(NMIX,IGR)+RXCOEF*OVERVB(IGR,I,IDF) + 249 CONTINUE +* COMPUTE DELTA-CONCENTRATION + DO 49 ISO=1,NISO + IF(DENSITB(ISO,I,IDF).EQ.0.0) GO TO 49 + WRITE(HMICRO,'(3A4)') (HISO((ISO-1)*3+IH),IH=1,3) + IF(HMICRO.EQ.'XE135') THEN + DEL(ISO)=XEN/(IRMAX-IRMIN+1)-DENSITB(ISO,I,IDF)*RXCOEF + 1 +DEL(ISO) + CX=DEL(ISO) + ELSE IF(HMICRO.EQ.'BMOD') THEN + ISOB=ISO +* +* ERROR IN CALCULATING BORON CONCENTRATION, +* BUT ADEQUATELY APPROXIMATED IF THE BORON REF EQ 0.0ppm +* +* DEL(ISO)=BOR/(IRMAX-IRMIN+1)- +* 1 DENSITB(ISO,I,IDF)*RXCOEF+DEL(ISO) + DEL(ISO)=-DENSITB(ISO,I,IDF)*RXCOEF+DEL(ISO) + CB=DEL(ISO) + ELSE IF(HMICRO.EQ.'CWAT') THEN + DEL(ISO)=(DC-1.0)*DENSITB(ISO,I,IDF) + D1C=DEL(ISO) + D2C=D1C*D1C + D2CB=D1C*D1C + ELSE IF(HMICRO.EQ.'MWAT') THEN + DMOD=DENSITB(ISO,I,IDF)/(IRMAX-IRMIN+1)+DMOD + DEL(ISO)=(DM-1.0)*DENSITB(ISO,I,IDF) + D1M=ALOG(DM) + D2M=1.0/(DM*DENSITB(ISO,I,IDF)) - 1.0/DENSITB(ISO,I,IDF) + D2MB=DEL(ISO) +* PURITY + PUR=(XI-XIR)*DM*DENSITB(ISO,I,IDF) + ELSE IF(HMICRO.EQ.'SM149') THEN + DEL(ISO)=SM/(IRMAX-IRMIN+1)-DENSITB(ISO,I,IDF)*RXCOEF + 1 +DEL(ISO) + CSM=DEL(ISO) + ELSE IF(HMICRO.EQ.'NP239') THEN + DEL(ISO)=RNP9/(IRMAX-IRMIN+1)-DENSITB(ISO,I,IDF)*RXCOEF + 1 +DEL(ISO) + CNP9=DEL(ISO) + ELSE IF(HMICRO.EQ.'FPC') THEN +* + CF= CF+DENSITB(ISO,I,IDF)*RXCOEF + DEL(ISO)=CPF1(II)*RXCOEF+DEL(ISO) + ELSE IF(HMICRO.EQ.'MACR') THEN + DEL(ISO)=DENSITB(ISO,I,IDF) + DO 271 IGR=1,NGRP + CHI(NMIX,IGR)=CHIB(IGR,I,IDF)*RXCOEF+CHI(NMIX,IGR) + 271 CONTINUE + ENDIF + 49 CONTINUE + 900 CONTINUE +* R.C. 24/05/2011 +* boron unit correction +* Bnat = 10.811 g/mol +* O16 = 15.9949 g/mol +* H1 = 1.0078 g/mol +* D2 = 2.0141 g/mol + CB=DMOD*BOR/10.811*((1.0078*(1-XI)+2.0141*XI)*2+15.9949)/3+CB + CB=CB*DM + DEL(ISOB)=CB +* R.C. +* CORRECT THE FUEL TEMPERATURE + CQ=REAL((CPF1(1)*XCOF(1)+CPF1(2)*XCOF(2)+CPF1(3)*XCOF(3)+CF)/CF) + RCQ=1.0-CQ + CQ2=CQ*CQ + RCQ2=1.0-CQ2 + IF(JTAB(1).EQ.0) THEN + CQ=0.0 + RCQ=0.0 + CQ2=0.0 + RCQ2=0.0 + CX=0.0 + CSM=0.0 + CNP9=0.0 + D1C=0.0 + D2C=0.0 + D2CB=0.0 + T1F=0.0 + T2F=0.0 + T1C=0.0 + T2C=0.0 + ENDIF +* +* RECOVER MACROSCOPIC X-SECTIONS + II=0 + DO 901 I = IRMIN,IRMAX + II=II+1 + RXCOEF=REAL(XCOF(II)) + CPF2=CPF1(II)*CPF1(II) + DO 98 ISO=1,NISO + DO 99 ITY=1,NTYP + IF(ISO.EQ.1) THEN + ZZ=1.0 + IF(ITY.EQ.1) THEN + VD2M=D2MB + VD2C=D2CB + ELSE IF(ITY.GE.5) THEN + VD2M=D2MB + VD2C=D2C + ELSE + VD2M=D2M + VD2C=D2C + ENDIF + ELSE + ZZ=0.0 + VD2M=D2M + VD2C=D2C + ENDIF + DO 100 IGR=1,NG2 + SMAC(IGR,ITY,1)=SMAC(IGR,ITY,1)+ + 1 SMACB(IGR,ITY,1,I,IDF)*RXCOEF*ZZ+ + 1 (XBORB(IGR,ITY,ISO,I,IDF)*RXCOEF*CB + + 1 XPURB(IGR,ITY,ISO,I,IDF)*RXCOEF*PUR*ZZ + + 1 XXENB(IGR,ITY,ISO,I,IDF)*RXCOEF*CX + + 1 XT1FB(IGR,ITY,ISO,I,IDF)*RXCOEF*T1F*CQ + + 1 XT1FB(IGR,ITY,ISO,1,IDF)*RXCOEF*T1F*RCQ + + 1 XT2FB(IGR,ITY,ISO,I,IDF)*RXCOEF*T2F*CQ2 + + 1 XT2FB(IGR,ITY,ISO,1,IDF)*RXCOEF*T2F*RCQ2 + + 1 XT1CB(IGR,ITY,ISO,I,IDF)*RXCOEF*T1C + + 1 XT2CB(IGR,ITY,ISO,I,IDF)*RXCOEF*T2C + + 1 XT1MB(IGR,ITY,ISO,I,IDF)*RXCOEF*T1M + + 1 XT2MB(IGR,ITY,ISO,I,IDF)*RXCOEF*T2M + + 1 XD1CB(IGR,ITY,ISO,I,IDF)*RXCOEF*D1C + + 1 XD2CB(IGR,ITY,ISO,I,IDF)*RXCOEF*VD2C + + 1 XD1MB(IGR,ITY,ISO,I,IDF)*RXCOEF*D1M + + 1 XD2MB(IGR,ITY,ISO,I,IDF)*RXCOEF*VD2M + + 1 XSMB(IGR,ITY,ISO,I,IDF)*RXCOEF*CSM + + 1 XNP9B(IGR,ITY,ISO,I,IDF)*RXCOEF*CNP9 + + 1 XMFDB(IGR,ITY,ISO,I,IDF)*RXCOEF*T1F*D1C + + 1 XMMDB(IGR,ITY,ISO,I,IDF)*RXCOEF*T1C*D1C + + 1 XPF1B(IGR,ITY,ISO,I,IDF)*RXCOEF*CPF1(II)*(1.-YF) + + 1 XPF2B(IGR,ITY,ISO,I,IDF)*RXCOEF*CPF2*(1.-YF)+ + 1 XPF1LB(IGR,ITY,ISO,I,IDF)*RXCOEF*CPF1(II)*YF + + 1 XPF2LB(IGR,ITY,ISO,I,IDF)*RXCOEF*CPF2*YF)*DEL(ISO) +* + 100 CONTINUE + 99 CONTINUE + 98 CONTINUE + 901 CONTINUE +* STORE SCATTERING + IL= 1 + ITY=5+2*IXYZ+IL + IGAR=0 + DO 130 JGR=1,NGRP + DO 120 IGR=IJJ(JGR),IJJ(JGR)-NJJ(JGR)+1,-1 + IGAR=IGAR+1 + SCAT(NMIX,IL,IGR,JGR)=SMAC(IGAR,ITY,1) +* TOTAL OR ABS + SMAC(IGR,2,1)=SMAC(IGR,2,1)+SCAT(NMIX,IL,IGR,JGR) + 120 CONTINUE + 130 CONTINUE +* STORE X-SECTIONS + DO 261 ITY=1,NTYP + DO 260 IGR=1,NGRP + SIGMA(NMIX,IGR,ITY)=SMAC(IGR,ITY,1) + 260 CONTINUE + 261 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(CPW2,CPW1,DEL,SMAC) + RETURN + END |