Initial commit from Polytechnique Montreal

1 files changed, 88 insertions, 0 deletions

diff --git a/Dragon/src/DUO005.f b/Dragon/src/DUO005.f
new file mode 100644
index 0000000..62dcac4
--- /dev/null
+++ b/Dragon/src/DUO005.f
@@ -0,0 +1,88 @@
+*DECK DUO005
+      SUBROUTINE DUO005(NGRP,DLK,ALK,FLUX,AFLUX,V,W)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute a consistent direct-adjoint leakage multigroup term using
+* the Villarino-Stamm'ler normalization method.
+*
+*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
+* NGRP    number of energy groups.
+* DLK     leakage rates for the direct balance equation.
+* ALK     leakage rates for the adjoint balance equation.
+* FLUX    integrated direct flux.
+* AFLUX   integrated adjoint flux.
+*
+*Parameters: output
+* V       first Villarino-Stamm'ler normalization vector.
+* W       second Villarino-Stamm'ler normalization vector.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NGRP
+      REAL DLK(NGRP),ALK(NGRP),FLUX(NGRP),AFLUX(NGRP),V(NGRP),W(NGRP)
+*----
+*  LOCAL VARIABLES
+*----
+      DOUBLE PRECISION ZNUMD,ZNUMA,GAR
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DFF,DFF2
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(DFF(2*NGRP+1,2*NGRP),DFF2(2*NGRP,2*NGRP+1))
+*
+      ZNUMD=0.0D0
+      ZNUMA=0.0D0
+      DO IGR=1,NGRP
+        ZNUMD=ZNUMD+FLUX(IGR)
+        ZNUMA=ZNUMA+AFLUX(IGR)
+      ENDDO
+      DFF(:2*NGRP+1,:2*NGRP)=0.0D0
+      DO IGR=1,NGRP
+        DFF(IGR,IGR)=ZNUMD
+        DFF(NGRP+IGR,NGRP+IGR)=ZNUMA
+        DO JGR=1,NGRP
+          DFF(IGR,NGRP+JGR)=FLUX(JGR)
+          DFF(NGRP+IGR,JGR)=AFLUX(JGR)
+        ENDDO
+        DFF(2*NGRP+1,IGR)=AFLUX(IGR)
+        DFF(2*NGRP+1,NGRP+IGR)=FLUX(IGR)
+      ENDDO
+      DFF2(:2*NGRP,:2*NGRP+1)=0.0D0
+      DO IGR=1,2*NGRP
+        GAR=0.0D0
+        DO JGR=1,NGRP
+          GAR=GAR+DFF(JGR,IGR)*DLK(JGR)+DFF(NGRP+JGR,IGR)*ALK(JGR)
+        ENDDO
+        DFF2(IGR,2*NGRP+1)=GAR
+        DO JGR=1,2*NGRP
+          DO KGR=1,2*NGRP+1
+            DFF2(IGR,JGR)=DFF2(IGR,JGR)+DFF(KGR,IGR)*DFF(KGR,JGR)
+          ENDDO
+        ENDDO
+      ENDDO
+      CALL ALSBD(2*NGRP,1,DFF2,IER,2*NGRP)
+      IF(IER.NE.0) CALL XABORT('DUO005: SINGULAR MATRIX.')
+      DO IGR=1,NGRP
+        V(IGR)=REAL(DFF2(IGR,2*NGRP+1))
+        W(IGR)=REAL(DFF2(NGRP+IGR,2*NGRP+1))
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(DFF2,DFF)
+      RETURN
+      END