summaryrefslogtreecommitdiff
path: root/Dragon/src/SPHDRV.f
blob: d6fee88cce4cf32558494a3c8f5104ef9ea84ca4 (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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
*DECK SPHDRV
      SUBROUTINE SPHDRV(IPTRK,IFTRK,IPMACR,IPFLX,IPRINT,IMC,NGCOND,
     1 NMERGE,NALBP,IGRMIN,IGRMAX,SPH)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Calculation of the SPH factors. These factors are used to multiply the
* cross sections and to divide the fluxes. The SPH factors calculation
* is generally application dependent. New SPH algorithms should be
* implemented in this driver.
*
*Copyright:
* Copyright (C) 2011 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): A. Hebert
*
*Parameters: input
* IPTRK   pointer to the macro-tracking LCM object.
* IFTRK   unit of the macro-tracking binary sequential file.
* IPMACR  pointer to the Macrolib (L_MACROLIB signature).
* IPFLX   pointer towards an initialization flux (L_FLUX signature).
* IPRINT  print flag (equal to 0 for no print).
* IMC     type of macro-calculation (=1 diffusion or SPN; 
*         =2 other options;
*         =3 type PIJ with Bell acceleration).
* NGCOND  number of condensed groups.
* NMERGE  number of merged regions.
* NALBP   number of physical albedos.
* IGRMIN  first group to process.
* IGRMAX  last group to process.
*
*Parameters: output
* SPH     SPH homogenization factors.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPTRK,IPMACR,IPFLX
      INTEGER   IFTRK,IPRINT,IMC,NGCOND,NMERGE,NALBP,IGRMIN,IGRMAX
      REAL      SPH(NMERGE+NALBP,NGCOND)
*----
*  LOCAL VARIABLES
*----
      PARAMETER (NSTATE=40)
      CHARACTER TEXT12*12,CNDOOR*12,CTITRE*72,SUFF(2)*2,HSMG*131
      INTEGER   ISTATE(NSTATE)
      LOGICAL   ILK
      TYPE(C_PTR) IPSPH
*----
*  ALLOCATABLE ARRAYS
*----
      INTEGER, ALLOCATABLE, DIMENSION(:) :: MAT2,KEY2,MERG2
      REAL, ALLOCATABLE, DIMENSION(:) :: VOL2,VREF,VMAC
      DATA      SUFF/'00','01'/
*----
*  RECOVER SPH-RELATED INFORMATION
*----
      CALL LCMLEN(IPMACR,'SPH',ILONG,ITYLCM)
      IF(ILONG.EQ.0) CALL XABORT('SPHDRV: MISSING SPH DIRECTORY.')
      IPSPH=LCMDID(IPMACR,'SPH')
      CALL LCMGET(IPSPH,'STATE-VECTOR',ISTATE)
      NSPH=ISTATE(1)
      KSPH=ISTATE(2)
      MAXIT=ISTATE(3)
      MAXNBI=ISTATE(4)
      IF((NSPH.EQ.0).OR.(NSPH.EQ.1)) CALL XABORT('SPHDRV: INVALID VALU'
     > //'E OF NSPH.')
*----
*  RECOVER AND USE AN EXISTING MACRO-TRACKING.
*----
      IF(C_ASSOCIATED(IPTRK)) THEN
         IF(NSPH.GE.2) THEN
            CALL LCMGTC(IPSPH,'SPH$TRK',12,CNDOOR)
            CALL LCMGET(IPSPH,'SPH-EPSILON',EPSPH)
         ENDIF
         CALL LCMGTC(IPTRK,'SIGNATURE',12,TEXT12)
         IF(TEXT12.NE.'L_TRACK') THEN
            CALL XABORT('SPHDRV: TRACKING DATA STRUCTURE EXPECTED.')
         ENDIF
         CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
         NREG2=ISTATE(1)
         NUN2=ISTATE(2)
         ILK=ISTATE(3).EQ.0
         CALL LCMGTC(IPTRK,'TRACK-TYPE',12,CNDOOR)
         IF(CNDOOR.EQ.'MCCG') THEN
            CALL LCMLEN(IPTRK,'KEYFLX',LKFL,ITYLCM)
            NFUNL=LKFL/NREG2
         ELSE
            NFUNL=1
         ENDIF
*----
*  THE NUMBER OF UNKNOWNS OF A CURRENT-BASED MULTICELL ITERATION IS
*  INCREASED TO HOLD INTERFACE CURRENT COMPONENTS.
*----    
         IF(NSPH.EQ.4) THEN
           IF(CNDOOR.EQ.'SYBIL') NUN2=NUN2+ISTATE(9)
           IF((CNDOOR.EQ.'EXCELL').AND.(ISTATE(7).EQ.5))
     >       NUN2=NUN2+ISTATE(28)
         ENDIF
*
         IF((CNDOOR.EQ.'EXCELL').OR.(CNDOOR.EQ.'MCCG')) THEN
            ISCAT=ISTATE(6)
         ELSE IF(CNDOOR.EQ.'SN') THEN
            ISCAT=ISTATE(16)
         ELSE IF(CNDOOR.EQ.'BIVAC') THEN
            IF(ISTATE(14).EQ.0) THEN
               ISCAT=1
            ELSE
               ISCAT=ISTATE(16)
            ENDIF
         ELSE IF(CNDOOR.EQ.'TRIVAC') THEN
            IF(ISTATE(30).EQ.0) THEN
               ISCAT=1
            ELSE
               ISCAT=ISTATE(32)
            ENDIF
         ELSE
            ISCAT=1
         ENDIF
         ISCAT=ABS(ISCAT)
         ALLOCATE(VOL2(NREG2),MAT2(NREG2),KEY2(NREG2*NFUNL))
         CALL LCMGET(IPTRK,'VOLUME',VOL2)
         CALL LCMGET(IPTRK,'MATCOD',MAT2)
         CALL LCMGET(IPTRK,'KEYFLX',KEY2)
         CALL LCMLEN(IPTRK,'TITLE',LENGT,ITYLCM)
         IF(LENGT.GT.0) THEN
            CALL LCMGTC(IPTRK,'TITLE',72,CTITRE)
         ELSE
            CTITRE='*** NO TITLE PROVIDED ***'
         ENDIF
         NBMIX2=0
         IF(KSPH.EQ.5) THEN
*          HEBERT-BENOIST ALBS TECHNIQUE.
           DO 20 IREG=1,NREG2
           NBMIX2=MAX(NBMIX2,MAT2(IREG))
   20      CONTINUE
           ALLOCATE(MERG2(NBMIX2))
           DO 25 IBM=1,NBMIX2
           MERG2(IBM)=1
   25      CONTINUE
           ILK=.FALSE.
         ELSE
           DO 30 IREG=1,NREG2
           NBMIX2=MAX(NBMIX2,MAT2(IREG))
   30      CONTINUE
           IF(NBMIX2.NE.NMERGE) THEN
              WRITE(HSMG,'(41HSPHDRV: INVALID NUMBER OF MACRO-REGIONS (,
     >        2I6,2H).)') NBMIX2,NMERGE
              CALL XABORT(HSMG)
           ENDIF
           ALLOCATE(MERG2(NBMIX2))
           DO 35 IBM=1,NBMIX2
           MERG2(IBM)=IBM
   35      CONTINUE
         ENDIF
*
*        RECOVER TABULATED FUNCTIONS.
         CALL XDRTA2
      ELSE
         ISCAT=1
         ILK=.FALSE.
         NBMIX2=0
         NREG2=0
         NUN2=0
      ENDIF
*----
*  CHECK VOLUME CONSISTENCY
*----
      ALLOCATE(VREF(NMERGE),VMAC(NMERGE))
      VMAC(:)=0.0
      CALL LCMGET(IPMACR,'VOLUME',VREF)
      DO IREG=1,NREG2
         IBM=MAT2(IREG)
         IF(IBM.GT.0) VMAC(IBM)=VMAC(IBM)+VOL2(IREG)
      ENDDO
      VREFT=SUM(VREF(:NMERGE))
      VMACT=SUM(VMAC(:NMERGE))
      DO IBM=1,NMERGE
         ERR=ABS(VREF(IBM)/VREFT-VMAC(IBM)/VMACT)
         IF(ERR.GT.1.0E-4*ABS(VREF(IBM)/VREFT)) THEN
            WRITE(HSMG,'(38HSPHDRV: INCONSISTENT VOLUME IN MIXTURE,I5,
     >      3H BY,F7.2,2H %)') IBM,ERR*100.0
            CALL XABORT(HSMG)
         ENDIF
      ENDDO
      DEALLOCATE(VMAC,VREF)
*----
*  GENERAL PROCEDURE FOR COMPUTING THE SPH FACTORS
*----
      CALL SPHEQU(NBMIX2,IPTRK,IFTRK,IPMACR,IPFLX,CNDOOR,NSPH,KSPH,
     1 MAXIT,MAXNBI,EPSPH,IPRINT,IMC,NGCOND,NMERGE,NALBP,ISCAT,NREG2,
     2 NUN2,MAT2,VOL2,KEY2,MERG2,ILK,CTITRE,IGRMIN,IGRMAX,SPH)
      IF(C_ASSOCIATED(IPTRK)) DEALLOCATE(MERG2,KEY2,MAT2,VOL2)
*----
*  PRINT SPH FACTORS
*----
      IF(IPRINT.GT.1) THEN
         WRITE(6,'(/21H SPHDRV: SPH FACTORS:)')
         WRITE(6,200) ((IKK,IGR,SPH(IKK,IGR),IKK=1,NMERGE+NALBP),IGR=1,
     >   NGCOND)
      ENDIF
      RETURN
*
  200 FORMAT(4X,4HSPH(,I5,1H,,I3,2H)=,F9.5,:,4X,4HSPH(,I5,1H,,I3,2H)=,
     > F9.5,:,4X,4HSPH(,I5,1H,,I3,2H)=,F9.5,:,4X,4HSPH(,I5,1H,,I3,2H)=,
     > F9.5,:,4X,4HSPH(,I5,1H,,I3,2H)=,F9.5)
      END