Initial commit from Polytechnique Montreal

1 files changed, 145 insertions, 0 deletions

diff --git a/Donjon/src/DETFLU.f b/Donjon/src/DETFLU.f
new file mode 100644
index 0000000..a2938bd
--- /dev/null
+++ b/Donjon/src/DETFLU.f
@@ -0,0 +1,145 @@
+*DECK DETFLU
+      SUBROUTINE DETFLU(LHEX,NX,NY,NZ,NEL,NUN,MESHX,MESHY,MESHZ,KEYF,
+     > FLUX,NGRP,SPEC,DEVPOS,NHEX,IHEX,RESP,IPRT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute flux at detector site
+*
+*Copyright:
+* Copyright (C) 2010 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* E. Varin, M. Guyot
+*
+*Parameters:
+* LHEX   =.TRUE. if hexagonal detectors are present
+* NX     number of x mesh-splitted elements 
+* NY     number of y mesh-splitted elements 
+* NZ     number of z mesh-splitted elements
+* NEL    number of finite elements
+* NUN    number of unknowns
+* MESHX  regions coordinates according to x
+* MESHY  regions coordinates according to y
+* MESHZ  regions coordinates according to z
+* KEYF   keyflux recover from L_TRACk object
+* FLUX   flux for each mesh-splitted elements
+* NGRP   number of energy groups
+* SPEC   spectral information
+* DEVPOS detector coordinates
+* NHEX   number of hexagons in the detector
+* IHEX   index number of hexagons
+* COR    center detector coordinates
+* RESP   flux reads by the detector
+* IPRT   printing index
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NX,NY,NZ,NEL,NUN,NGRP,IPRT,NHEX,NT,KEYF(NEL),IHEX(NHEX)
+      REAL MESHX(NX+1),MESHY(NY+1),MESHZ(NZ+1),FLUX(NUN,NGRP),RESP,
+     1     DEVPOS(6),SPEC(NGRP)
+      LOGICAL LHEX
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER NXP1,NYP1,NZP1,I,J,K,I1,I2,J1,J2,K1,K2,IAM,IGR
+      REAL    X1,X2,Y1,Y2,Z1,Z2
+
+      NXP1 = NX+1
+      NYP1 = NY+1
+      NZP1 = NZ+1
+
+      X1=DEVPOS(1)
+      X2=DEVPOS(2)
+      Y1=DEVPOS(3)
+      Y2=DEVPOS(4)
+      Z1=DEVPOS(5)
+      Z2=DEVPOS(6)
+
+      IF(.NOT.LHEX) THEN
+        IF(X1.LT.MESHX(1)) X1=MESHX(1)
+        IF(X2.LT.MESHX(1)) X2=MESHX(1)
+        IF(X2.GT.MESHX(NXP1)) X2=MESHX(NXP1)
+        IF(X1.GT.MESHX(NXP1)) X1=MESHX(NXP1)
+
+        IF(Y1.LT.MESHY(1)) Y1=MESHY(1)
+        IF(Y2.LT.MESHY(1)) Y2=MESHY(1)
+        IF(Y2.GT.MESHY(NYP1)) Y2=MESHY(NYP1)
+        IF(Y1.GT.MESHY(NYP1)) Y1=MESHY(NYP1)
+      ENDIF
+
+      IF(Z1.LT.MESHZ(1)) Z1=MESHZ(1)
+      IF(Z2.LT.MESHZ(1)) Z2=MESHZ(1)
+      IF(Z2.GT.MESHZ(NZP1)) Z2=MESHZ(NZP1)
+      IF(Z1.GT.MESHZ(NZP1)) Z1=MESHZ(NZP1)
+
+      IF(.NOT.LHEX) THEN
+        I1=0
+        DO 20 I=1,NXP1
+          IF(X1.GE.MESHX(I) .AND. X1.LE.MESHX(I+1)) THEN
+            I1=I
+          ENDIF
+          IF(X2.GE.MESHX(I) .AND. X2.LE.MESHX(I+1)) THEN
+            I2=I
+            GOTO 10
+          ENDIF
+  20    CONTINUE
+
+  10    DO 30 J=1,NYP1
+          IF(Y1.GE.MESHY(J) .AND. Y1.LE.MESHY(J+1)) THEN
+              J1=J
+          ENDIF
+          IF(Y2.GE.MESHY(J) .AND. Y2.LE.MESHY(J+1)) THEN
+            J2=J
+            GOTO 40
+          ENDIF
+  30    CONTINUE
+  40    CONTINUE
+      ELSE
+        J1 = 1
+        J2 = 1
+        I1 = 1
+        I2 = NHEX
+      ENDIF
+
+      DO 50 K=1,NZP1
+          IF(Z1.GE.MESHZ(K) .AND. Z1.LE.MESHZ(K+1)) THEN
+             K1=K
+          ENDIF
+          IF(Z2.GE.MESHZ(K) .AND. Z2.LE.MESHZ(K+1)) THEN
+             K2=K
+             GOTO 60
+          ENDIF
+  50  CONTINUE
+
+  60  RESP = 0.0
+      NT = 0
+
+      IF(IPRT.GT.4) WRITE(6,*) 'POS GEOM ',I1,I2,J1,J2,K1,K2
+      DO 70 K=K1,K2
+       DO 71 J=J1,J2
+        DO 72 I=I1,I2
+          NT = NT+1
+          IF(LHEX) THEN
+            IAM = (K-1)*NX+IHEX(I)
+          ELSE
+            IAM=(K-1)*NX*NY+(J-1)*NX+I
+          ENDIF
+          DO 73 IGR=1,NGRP
+            RESP = RESP + SPEC(IGR)*FLUX(KEYF(IAM),IGR)
+  73      CONTINUE
+          IF(IPRT.GT.4) WRITE(6,*) 'DETFLU: FINITE ELEMENT NUMBER ',
+     +                             IAM
+  72  CONTINUE
+  71  CONTINUE
+  70  CONTINUE
+
+      RESP = RESP / FLOAT(NT)
+
+      RETURN
+      END