Initial commit from Polytechnique Montreal

1 files changed, 177 insertions, 0 deletions

diff --git a/Dragon/src/LIBCOR.f b/Dragon/src/LIBCOR.f
new file mode 100644
index 0000000..2fe872d
--- /dev/null
+++ b/Dragon/src/LIBCOR.f
@@ -0,0 +1,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