summaryrefslogtreecommitdiff
path: root/Donjon/src/DETFLU.f
blob: a2938bdeee3f65e1ca92e537139e4fa3d0ea6e77 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
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