Initial commit from Polytechnique Montreal

1 files changed, 206 insertions, 0 deletions

diff --git a/Dragon/src/DUO006.f b/Dragon/src/DUO006.f
new file mode 100644
index 0000000..c9ace05
--- /dev/null
+++ b/Dragon/src/DUO006.f
@@ -0,0 +1,206 @@
+*DECK DUO006
+      SUBROUTINE DUO006(IPLIB,IPRINT,NISOT,NGRP,HREAC,IDIV,RHS,
+     > FLUX,AFLUX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Processing one of the two microlibs and return the RHS matrix for
+* the single reaction HREAC.
+*
+*Copyright:
+* Copyright (C) 2013 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   microlib.
+* IPRINT  print parameter.
+* NISOT   number of isotopes.
+* NGRP    number of energy groups.
+* HREAC   character*8 reaction name of the reaction to process.
+* IDIV    type of divergence term processing (=0: no processing;
+*         =1: direct processing; =2: adjoint processing;
+*         =3: direct-adjoint processing).
+*
+*Parameters: output
+* RHS     macroscopic cross-section matrix.
+* FLUX    integrated direct flux.
+* AFLUX   integrated adjoint flux.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPLIB
+      INTEGER IPRINT,NISOT,NGRP,IDIV
+      CHARACTER HREAC*8
+      REAL RHS(NGRP,NGRP,NISOT),FLUX(NGRP,NISOT),AFLUX(NGRP,NISOT)
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) KPLIB
+      CHARACTER HSMG*131
+      DOUBLE PRECISION SUM
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IHUSED
+      REAL, ALLOCATABLE, DIMENSION(:) :: DENS,VOL,VECTOR,CHI,SIGS,DLK,
+     > ALK,V,W
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: SCAT,RATE
+      TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IHUSED(3,NISOT))
+      ALLOCATE(DENS(NISOT),VOL(NISOT),VECTOR(NGRP),SCAT(NGRP,NGRP),
+     > CHI(NGRP))
+      ALLOCATE(IPISO(NISOT))
+*----
+*  FIND ISOTOPE POINTERS IN INPUT MICROLIB
+*----
+      CALL LIBIPS(IPLIB,NISOT,IPISO)
+*----
+*  COMPUTE THE RHS AND LHS MATRICES
+*----
+      CALL LCMGET(IPLIB,'ISOTOPESUSED',IHUSED)
+      CALL LCMGET(IPLIB,'ISOTOPESDENS',DENS)
+      CALL LCMGET(IPLIB,'ISOTOPESVOL',VOL)
+      CALL LCMGET(IPLIB,'K-EFFECTIVE',ZKEFF)
+      IF(IPRINT.GT.4) WRITE(6,'(35H DUO006: EFFECTIVE MULTIPLICATION F,
+     > 6HACTOR=,1P,E12.5)') ZKEFF
+      RHS(:NGRP,:NGRP,:NISOT)=0.0
+      DO ISOT=1,NISOT
+        IF(IPRINT.GT.4) WRITE(6,'(29H DUO006: PROCESSING ISOTOPE '',
+     >  3A4,2H''.)') (IHUSED(I0,ISOT),I0=1,3)
+        KPLIB=IPISO(ISOT) ! set ISOT-th isotope
+        IF(.NOT.C_ASSOCIATED(KPLIB)) THEN
+          WRITE(HSMG,'(18H DUO006: ISOTOPE '',3A4,7H'' (ISO=,I8,4H) IS,
+     >    31H NOT AVAILABLE IN THE MICROLIB.)') (IHUSED(I0,ISOT),
+     >    I0=1,3),ISOT
+          CALL XABORT(HSMG)
+        ENDIF
+        CALL LCMGET(KPLIB,'NWT0',FLUX(1,ISOT))
+        CALL LCMLEN(KPLIB,'NWAT0',ILON,ITYLCM)
+        IF(ILON.NE.0) THEN
+          CALL LCMGET(KPLIB,'NWAT0',AFLUX(1,ISOT))
+        ELSE
+          AFLUX(:NGRP,ISOT)=1.0
+        ENDIF
+        DO IGR=1,NGRP
+          FLUX(IGR,ISOT)=FLUX(IGR,ISOT)*VOL(ISOT)
+          AFLUX(IGR,ISOT)=AFLUX(IGR,ISOT)*VOL(ISOT)
+        ENDDO
+        CALL LCMLEN(KPLIB,HREAC,ILONG,ITYLCM)
+        IF((ILONG.EQ.0).AND.(HREAC.NE.'LEAK')) CYCLE
+        IF(HREAC.EQ.'SCAT00') THEN
+          ALLOCATE(SIGS(NGRP))
+          CALL XDRLGS(KPLIB,-1,IPRINT,0,0,1,NGRP,SIGS,SCAT,ITYPRO)
+          DEALLOCATE(SIGS)
+          DO IGR=1,NGRP
+            DO JGR=1,NGRP
+              RHS(JGR,IGR,ISOT)=RHS(JGR,IGR,ISOT)+DENS(ISOT)*
+     >        SCAT(JGR,IGR)
+            ENDDO
+          ENDDO
+        ELSE IF((HREAC.EQ.'NUSIGF').OR.(HREAC.EQ.'CHI')) THEN
+          CALL LCMGET(KPLIB,'NUSIGF',VECTOR)
+          CALL LCMGET(KPLIB,'CHI',CHI)
+          DO IGR=1,NGRP
+            DO JGR=1,NGRP
+              RHS(JGR,IGR,ISOT)=RHS(JGR,IGR,ISOT)+DENS(ISOT)*CHI(JGR)*
+     >        VECTOR(IGR)
+            ENDDO
+          ENDDO
+        ELSE IF(HREAC(:3).EQ.'NWT') THEN
+          WRITE(HSMG,'(8HDUO006: ,A8,25H IS A FORBIDDEN REACTION.)')
+     >    HREAC
+          CALL XABORT(HSMG)
+        ELSE IF(HREAC.EQ.'LEAK') THEN
+          ALLOCATE(RATE(NGRP,NGRP))
+          RATE(:NGRP,:NGRP)=0.0
+          CALL LCMLEN(KPLIB,'NUSIGF',ILON,ITYLCM)
+          IF(ILON.GT.0) THEN
+            CALL LCMGET(KPLIB,'NUSIGF',VECTOR)
+            CALL LCMGET(KPLIB,'CHI',CHI)
+            DO IGR=1,NGRP
+              DO JGR=1,NGRP
+                RATE(JGR,IGR)=RATE(JGR,IGR)-DENS(ISOT)*CHI(JGR)*
+     >          VECTOR(IGR)/ZKEFF
+              ENDDO
+            ENDDO
+          ENDIF
+          ALLOCATE(SIGS(NGRP))
+          CALL XDRLGS(KPLIB,-1,IPRINT,0,0,1,NGRP,SIGS,SCAT,ITYPRO)
+          DEALLOCATE(SIGS)
+          CALL LCMGET(KPLIB,'NTOT0',VECTOR)
+          DO IGR=1,NGRP
+            DO JGR=1,NGRP
+              RATE(JGR,IGR)=RATE(JGR,IGR)-DENS(ISOT)*SCAT(JGR,IGR)
+            ENDDO
+            RATE(IGR,IGR)=RATE(IGR,IGR)+DENS(ISOT)*VECTOR(IGR)
+          ENDDO
+          IF(IDIV.EQ.1) THEN
+            DO JGR=1,NGRP
+              SUM=0.0D0
+              DO IGR=1,NGRP
+                SUM=SUM+RATE(JGR,IGR)*FLUX(IGR,ISOT)
+              ENDDO
+              RHS(JGR,JGR,ISOT)=-REAL(SUM)/FLUX(JGR,ISOT)
+            ENDDO
+          ELSE IF(IDIV.EQ.2) THEN
+            DO IGR=1,NGRP
+              SUM=0.0D0
+              DO JGR=1,NGRP
+                SUM=SUM+RATE(JGR,IGR)*AFLUX(JGR,ISOT)
+              ENDDO
+              RHS(IGR,IGR,ISOT)=-REAL(SUM)/AFLUX(IGR,ISOT)
+            ENDDO
+          ELSE IF(IDIV.EQ.3) THEN
+            ALLOCATE(DLK(NGRP),ALK(NGRP))
+            DO JGR=1,NGRP
+              SUM=0.0D0
+              DO IGR=1,NGRP
+                SUM=SUM+RATE(JGR,IGR)*FLUX(IGR,ISOT)
+              ENDDO
+              DLK(JGR)=REAL(SUM)
+            ENDDO
+            DO IGR=1,NGRP
+              SUM=0.0D0
+              DO JGR=1,NGRP
+                SUM=SUM+RATE(JGR,IGR)*AFLUX(JGR,ISOT)
+              ENDDO
+              ALK(IGR)=REAL(SUM)
+            ENDDO
+            ALLOCATE(V(NGRP),W(NGRP))
+            CALL DUO005(NGRP,DLK,ALK,FLUX(1,ISOT),AFLUX(1,ISOT),V,W)
+            DO IGR=1,NGRP
+              DO JGR=1,NGRP
+                RHS(IGR,JGR,ISOT)=RHS(IGR,JGR,ISOT)-V(IGR)-W(JGR)
+              ENDDO
+            ENDDO
+            DEALLOCATE(W,V,ALK,DLK)
+          ENDIF
+          DEALLOCATE(RATE)
+        ELSE
+          CALL LCMGET(KPLIB,HREAC,VECTOR)
+          DO IGR=1,NGRP
+            RHS(IGR,IGR,ISOT)=RHS(IGR,IGR,ISOT)+DENS(ISOT)*VECTOR(IGR)
+          ENDDO
+        ENDIF
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IPISO)
+      DEALLOCATE(CHI,SCAT,VECTOR,VOL,DENS)
+      DEALLOCATE(IHUSED)
+      RETURN
+      END