summaryrefslogtreecommitdiff
path: root/Dragon/src/LIBCOR.f
blob: 2fe872df31421bb5ea5963b3f34efa2bff77603e (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
*DECK LIBCOR
      SUBROUTINE LIBCOR (IPLIB,NGRO,ISOT,JSOT,HNAMIS1,HNAMIS2)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the correlation information between a pair of resonant
* isotopes for the CALENDF method.
*
*Copyright:
* Copyright (C) 2003 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
* IPLIB   pointer to the internal library (L_LIBRARY signature).
* NGRO    number of energy groups.
* ISOT    position in list of the first isotope.
* JSOT    position in list of the second isotope.
* HNAMIS1 local name of the first isotope:
*         HNAMIS1(1:8)  is the local isotope name;
*         HNAMIS1(9:12) is a suffix function of the mixture index.
* HNAMIS2 local name of the second isotope
*         HNAMIS2(1:8)  is the local isotope name;
*         HNAMIS2(9:12) is a suffix function of the mixture index.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPLIB
      INTEGER NGRO,ISOT,JSOT
      CHARACTER HNAMIS1*12,HNAMIS2*12
*----
*  LOCAL VARIABLES
*----
      PARAMETER(MAXNPT=12)
      TYPE(C_PTR) JPLIB,IP1,IP2,JP1,JP2,KP1,KP2
      REAL SIGQT1(MAXNPT),SIGQT2(MAXNPT),WSLD1(MAXNPT**2),
     1 WSLD2(MAXNPT**2)
      DOUBLE PRECISION SUMA1,SUMB1,SUMA2,SUMB2
      INTEGER, ALLOCATABLE, DIMENSION(:) :: NFS1,NFS2
      REAL, ALLOCATABLE, DIMENSION(:) :: TBIN1,TBIN2,EBIN,DBIN,PROB1,
     1 PROB2
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: COMOM
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(NFS1(NGRO),NFS2(NGRO))
*
      JPLIB=LCMGID(IPLIB,'ISOTOPESLIST')
      IP1=LCMGIL(JPLIB,ISOT) ! set ISOT-th isotope
      IP2=LCMGIL(JPLIB,JSOT) ! set JSOT-th isotope
      CALL LCMLEN(IP1,'BIN-NFS',LENGT1,ITYLCM)
      CALL LCMLEN(IP2,'BIN-NFS',LENGT2,ITYLCM)
      IF((LENGT1.EQ.0).OR.(LENGT1.NE.LENGT2)) CALL XABORT('LIBCOR: UNA'
     1 //'BLE TO FIND CONSISTENT BIN TYPE INFORMATION.')
      CALL LCMGET(IP1,'BIN-NFS',NFS1)
      CALL LCMGET(IP2,'BIN-NFS',NFS2)
      LBIN=0
      IGRMIN=1
      IGRMAX=NGRO
      DO 10 IGRP=NGRO,1,-1
      IF(NFS1(IGRP).NE.NFS2(IGRP)) CALL XABORT('INVALID BIN INFO.')
      IF((IGRMAX.EQ.IGRP).AND.(NFS1(IGRP).EQ.0)) IGRMAX=IGRP-1
      LBIN=LBIN+NFS1(IGRP)
   10 CONTINUE
      DO 20 IGRP=1,NGRO
      IF((IGRMIN.EQ.IGRP).AND.(NFS1(IGRP).EQ.0)) IGRMIN=IGRP+1
   20 CONTINUE
      ALLOCATE(TBIN1(LBIN),TBIN2(LBIN),EBIN(LBIN+1),DBIN(LBIN))
      CALL LCMGET(IP1,'BIN-ENERGY',EBIN)
      CALL LCMGET(IP1,'BIN-NTOT0',TBIN1)
      CALL LCMGET(IP2,'BIN-NTOT0',TBIN2)
      CALL LCMSIX(IP1,'PT-TABLE',1)
      CALL LCMSIX(IP2,'PT-TABLE',1)
*---
*  LOOP OVER THE RESONANT ENERGY GROUPS.
*---
      LBIN=0
      JP1=LCMGID(IP1,'GROUP-PT')
      JP2=LCMGID(IP2,'GROUP-PT')
      DO 130 IGRP=IGRMIN,IGRMAX
      SUMA1=0.0D0
      SUMB1=0.0D0
      SUMA2=0.0D0
      SUMB2=0.0D0
      DO 30 IGF=1,NFS1(IGRP)
      SIGTA=MAX(0.002,TBIN1(LBIN+IGF))
      SIGTB=MAX(0.002,TBIN2(LBIN+IGF))
      DELM=LOG(EBIN(LBIN+IGF)/EBIN(LBIN+IGF+1))
      SUMA1=SUMA1+TBIN1(LBIN+IGF)*DELM
      SUMB1=SUMB1+SIGTA*DELM
      SUMA2=SUMA2+TBIN2(LBIN+IGF)*DELM
      SUMB2=SUMB2+SIGTB*DELM
      TBIN1(LBIN+IGF)=SIGTA
      TBIN2(LBIN+IGF)=SIGTB
      DBIN(LBIN+IGF)=DELM
   30 CONTINUE
      DO 40 IGF=1,NFS1(IGRP)
      TBIN1(LBIN+IGF)=TBIN1(LBIN+IGF)*REAL(SUMA1/SUMB1)
      TBIN2(LBIN+IGF)=TBIN2(LBIN+IGF)*REAL(SUMA2/SUMB2)
   40 CONTINUE
*
      CALL LCMLEL(JP1,IGRP,N1,ITYLCM)
      CALL LCMLEL(JP2,IGRP,N2,ITYLCM)
      IF((N1.EQ.0).OR.(N2.EQ.0)) GO TO 120
      KP1=LCMGIL(JP1,IGRP)
      KP2=LCMGIL(JP2,IGRP)
      CALL LCMLEN(KP1,'SIGQT-SIGS',NQT1,ITYLCM)
      CALL LCMLEN(KP2,'SIGQT-SIGS',NQT2,ITYLCM)
      CALL LCMLEN(KP1,'PROB-TABLE',NQT10,ITYLCM)
      CALL LCMLEN(KP2,'PROB-TABLE',NQT20,ITYLCM)
      ALLOCATE(PROB1(NQT10),PROB2(NQT20))
      CALL LCMGET(KP1,'PROB-TABLE',PROB1)
      CALL LCMGET(KP2,'PROB-TABLE',PROB2)
      DO 50 I=1,NQT1
      SIGQT1(I)=PROB1(MAXNPT+I)
   50 CONTINUE
      DO 60 I=1,NQT2
      SIGQT2(I)=PROB2(MAXNPT+I)
   60 CONTINUE
*
      ALLOCATE(COMOM(NQT1*NQT2))
      CALL LIBCOM(NFS1(IGRP),DBIN(LBIN+1),TBIN1(LBIN+1),
     1 TBIN2(LBIN+1),NQT1,NQT2,COMOM)
      CALL LIBOMG(NQT1,1-NQT1/2,SIGQT1,NQT2,1-NQT2/2,SIGQT2,
     1 COMOM,WSLD1)
      DEALLOCATE(COMOM)
*---
*  CHECK NORMALIZATION OF THE CORRELATED WEIGHT MATRIX.
*---
      DO 80 I=1,NQT1
      SUM=0.0
      DO 70 J=1,NQT2
      SUM=SUM+WSLD1((J-1)*NQT1+I)
   70 CONTINUE
      IF(ABS(SUM-PROB1(I)).GT.1.0E-4) THEN
         CALL XABORT('LIBCOR: BAD NORMALIZATION EXCEPTION(1).')
      ENDIF
   80 CONTINUE
      DO 100 I=1,NQT2
      SUM=0.0
      DO 90 J=1,NQT1
      SUM=SUM+WSLD1((I-1)*NQT1+J)
   90 CONTINUE
      IF(ABS(SUM-PROB2(I)).GT.1.0E-4) THEN
         CALL XABORT('LIBCOR: BAD NORMALIZATION EXCEPTION(2).')
      ENDIF
  100 CONTINUE
      DEALLOCATE(PROB2,PROB1)
*
      CALL LCMPUT(KP1,HNAMIS2,NQT1*NQT2,2,WSLD1)
      DO 115 I=1,NQT1
      DO 110 J=1,NQT2
      WSLD2((I-1)*NQT2+J)=WSLD1((J-1)*NQT1+I)
  110 CONTINUE
  115 CONTINUE
      CALL LCMPUT(KP2,HNAMIS1,NQT2*NQT1,2,WSLD2)
  120 LBIN=LBIN+NFS1(IGRP)
  130 CONTINUE
*
      CALL LCMSIX(IP2,' ',2)
      CALL LCMSIX(IP1,' ',2)
      DEALLOCATE(DBIN,EBIN,TBIN2,TBIN1)
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(NFS2,NFS1)
      RETURN
      END