Initial commit from Polytechnique Montreal

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