Initial commit from Polytechnique Montreal

1 files changed, 421 insertions, 0 deletions

diff --git a/Dragon/src/FMTSUD.f b/Dragon/src/FMTSUD.f
new file mode 100644
index 0000000..b41b72d
--- /dev/null
+++ b/Dragon/src/FMTSUD.f
@@ -0,0 +1,421 @@
+*DECK FMTSUD
+      SUBROUTINE FMTSUD(NENTRY,IENTRY,KENTRY,IPRINT,NOPT,IOPT,IKFLU,
+     >                  NREG  ,NGROUP,NDIM  ,POLOAQ,AZMOAQ,
+     >                  VOLUME,XPOL  ,WPOL  ,XAZI  ,WAZI  ,
+     >                  FLUX  ,AFLUX ,TFLUX ,WGHT  ,MU    ,ETA   )
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* To process the angular fluxes and generate the SUS3D file.
+*
+*Copyright:
+* Copyright (C) 2008 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):
+* G. Marleau
+*
+*Parameters: input
+* NENTRY  number of data structures transfered to this module.
+* IENTRY  data structure type where:
+*         =1 for LCM memory object;
+*         =2 for XSM file;
+*         =3 for sequential binary file;
+*         =4 for sequential ASCII file.
+* KENTRY  data structure pointer.
+* IPRINT  print level.
+* NOPT    number of options.
+* IOPT    processing option.
+* IKFLU   pointer to the FLUX data structure.
+* NREG    number of regions for problem.
+* NGROUP  number of groups for problem.
+* NDIM    number of dimensions of problem.
+* POLOAQ  polar quadrature order.
+* AZMOAQ  azimuthal quadrature order.
+* VOLUME  regional volumes.
+* XPOL    polar quadrature points.
+* WPOL    polar quadrature weights.
+* XAZI    azimuthal quadrature points.
+* WAZI    azimuthal quadrature weights.
+*
+*Parameters: output
+* FLUX    direct and adjoint flux.
+* AFLUX   angular components of the direct and adjoint flux.
+* TFLUX   temporary flux vector.
+* WGHT    temporary weight.
+* MU      temporary mu.
+* ETA     temporary eta.
+*
+*----------
+*
+      USE              GANLIB
+      IMPLICIT         NONE
+*----
+*  Subroutine arguments
+*----
+      TYPE(C_PTR)      KENTRY(NENTRY)
+      INTEGER          NENTRY,IENTRY(NENTRY)
+      INTEGER          IPRINT,NOPT,IOPT(NOPT),IKFLU
+      INTEGER          NREG,NGROUP,NDIM,POLOAQ,AZMOAQ
+      REAL             VOLUME(NREG),XPOL(POLOAQ,2),WPOL(POLOAQ),
+     >                 XAZI(NDIM,AZMOAQ),WAZI(AZMOAQ)
+      REAL             FLUX(NREG,2,NGROUP,2)
+      DOUBLE PRECISION AFLUX(NREG,POLOAQ,AZMOAQ*2,2,NGROUP)
+      DOUBLE PRECISION TFLUX(NREG+1,POLOAQ,AZMOAQ*2)
+      REAL             WGHT(POLOAQ*AZMOAQ*2),MU(POLOAQ*AZMOAQ*2),
+     >                 ETA(POLOAQ*AZMOAQ*2)
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='FMTSUD')
+      INTEGER          ILCMUP,ILCMDN
+      PARAMETER       (ILCMUP=1,ILCMDN=2)
+*----
+*  Local variables
+*----
+      TYPE(C_PTR)      IPU
+      INTEGER          IFPU,IAPU,IFORM,IGROUP,IKU,IAKU,NAZ,IA,IP,IR,IQUA
+      CHARACTER*12     NAMREC
+      REAL             RELERR,ERRMAX,ETATMP,FNORM
+*----
+*  Initialize FLUX vectors
+*----
+      NAZ=AZMOAQ
+      DO IP=1,POLOAQ
+        XPOL(IP,2)=SQRT(1.-XPOL(IP,1)*XPOL(IP,1))
+      ENDDO
+      IF(IPRINT .GE. 10) THEN
+        WRITE(IOUT,6000) NAMSBR,NGROUP,NDIM,2*AZMOAQ,POLOAQ,NREG
+        WRITE(IOUT,6001)
+        WRITE(IOUT,6011)
+     >  (XAZI(1,IA),XAZI(2,IA),WAZI(IA),IA=1,AZMOAQ)
+        WRITE(IOUT,6002)
+        WRITE(IOUT,6011)
+     >  (XPOL(IP,1),XPOL(IP,2),WPOL(IP),IP=1,POLOAQ)
+      ENDIF
+      FLUX(:NREG,:2,:NGROUP,:2)=0.0
+      AFLUX(:NREG,:POLOAQ,:AZMOAQ*2,:2,:NGROUP)=0.0D0
+*----
+*  Get information from FLUX data structure.
+*  1. Flux
+*  2. Angular flux
+*  3. Adjoint
+*  4. Angular adjoint
+*----
+      IPU=KENTRY(IKFLU)
+      CALL LCMSIX(IPU,'FLUXDIRECT  ',ILCMUP)
+      DO IGROUP=1,NGROUP
+        WRITE(NAMREC,'(A4,I3,5X)') 'FLUX',IGROUP
+        CALL LCMGET(IPU,NAMREC,FLUX(1,1,IGROUP,1))
+      ENDDO
+      CALL LCMSIX(IPU,'ANGULAR DIR ',ILCMUP)
+      DO IGROUP=1,NGROUP
+        WRITE(NAMREC,'(A5,I3,4X)') 'AFLUX',IGROUP
+        CALL LCMGET(IPU,NAMREC,AFLUX(1,1,1,1,IGROUP))
+      ENDDO
+      CALL LCMSIX(IPU,'ANGULAR DIR ',ILCMDN)
+      CALL LCMSIX(IPU,'FLUXDIRECT  ',ILCMDN)
+      CALL LCMSIX(IPU,'FLUXADJOINT ',ILCMUP)
+      DO IGROUP=1,NGROUP
+        WRITE(NAMREC,'(A4,I3,5X)') 'FLUX',IGROUP
+        CALL LCMGET(IPU,NAMREC,FLUX(1,2,IGROUP,1))
+      ENDDO
+      CALL LCMSIX(IPU,'ANGULAR ADJ ',ILCMUP)
+      DO IGROUP=1,NGROUP
+        WRITE(NAMREC,'(A5,I3,4X)') 'AFLUX',IGROUP
+        CALL LCMGET(IPU,NAMREC,AFLUX(1,1,1,2,IGROUP))
+      ENDDO
+      CALL LCMSIX(IPU,'ANGULAR ADJ ',ILCMDN)
+      CALL LCMSIX(IPU,'FLUXADJOINT ',ILCMDN)
+*----
+*  Create first SUS file
+*  Volume, and tracking directions
+*----
+      IFPU=FILUNIT(KENTRY(1))
+      WRITE(IFPU,1000) NREG
+      WRITE(IFPU,1001) 2*AZMOAQ*POLOAQ
+*----
+*  -\Omega (\varphi+\pi,\pi-\theta)
+*----
+      IQUA=0
+      DO IA=1,AZMOAQ
+        ETATMP=XAZI(2,NAZ+1-IA)
+        IF(ETATMP .EQ. 0.0) THEN
+          ETATMP=1.0E-10
+        ENDIF
+        DO IP=1,POLOAQ
+          IQUA=IQUA+1
+          MU(IQUA)=-XPOL(IP,2)*XAZI(1,NAZ+1-IA)
+          ETA(IQUA)=-XPOL(IP,2)*ETATMP
+          WGHT(IQUA)=WPOL(IP)/WAZI(NAZ+1-IA)
+        ENDDO
+      ENDDO
+*----
+*  +\Omega (\varphi,\theta)
+*----
+      DO IA=1,AZMOAQ
+        ETATMP=XAZI(2,IA)
+        IF(ETATMP .EQ. 0.0) THEN
+          ETATMP=1.0E-10
+        ENDIF
+        DO IP=1,POLOAQ
+          IQUA=IQUA+1
+          MU(IQUA)=XPOL(IP,2)*XAZI(1,IA)
+          ETA(IQUA)=XPOL(IP,2)*ETATMP
+          WGHT(IQUA)=WPOL(IP)/WAZI(IA)
+        ENDDO
+      ENDDO
+      WRITE(IFPU,1010) 'wght',
+     >                (WGHT(IQUA),IQUA=1,POLOAQ*AZMOAQ*2)
+      WRITE(IFPU,1010) 'mu  ',
+     >                (MU(IQUA),IQUA=1,POLOAQ*AZMOAQ*2)
+      WRITE(IFPU,1010) 'eta ',
+     >                (ETA(IQUA),IQUA=1,POLOAQ*AZMOAQ*2)
+      WRITE(IFPU,1010) 'vol ',
+     >                (VOLUME(IR),IR=1,NREG)
+*----
+*  Print angular flux
+*----
+      IFPU=FILUNIT(KENTRY(2))
+      IKU=IENTRY(2)
+      IFORM=IOPT(2)
+      ERRMAX=0.0
+      IF(IKU .EQ. 3) THEN
+        DO IGROUP=1,NGROUP
+          WRITE(IFPU) 'Group',IGROUP
+        ENDDO
+        DO IGROUP=NGROUP+1,NGROUP+6
+          WRITE(IFPU) 'Comments'
+        ENDDO
+      ELSE
+        DO IGROUP=1,NGROUP
+          WRITE(IFPU,'(A8,4X,I8)') 'Groups =',IGROUP
+        ENDDO
+        DO IGROUP=NGROUP+1,NGROUP+6
+          WRITE(IFPU,'(A8)') 'Comments'
+        ENDDO
+      ENDIF
+      DO IGROUP=1,NGROUP
+        DO IA=1,AZMOAQ*2
+          DO IP=1,POLOAQ
+            IF(IFORM .EQ. 1) THEN
+              DO IR=1,NREG
+                TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,1,IGROUP)
+                IF(IA .LE. AZMOAQ) THEN
+                  FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+                ELSE
+                  FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+                ENDIF
+              ENDDO
+              IF(IPRINT .GE. 100) THEN
+                WRITE(IOUT,6003) IP,IA,IGROUP
+                WRITE(IOUT,6010) (TFLUX(IR,IP,IA),IR=1,NREG)
+              ENDIF
+            ELSE
+              IR=1
+              TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,1,IGROUP)
+              DO IR=1,NREG
+                TFLUX(IR+1,IP,IA)=
+     >            2.0D0*AFLUX(IR,IP,IA,1,IGROUP)-TFLUX(IR,IP,IA)
+                IF(IA .LE. AZMOAQ) THEN
+                  FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+                ELSE
+                  FLUX(IR,1,IGROUP,2)=FLUX(IR,1,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,1,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+                ENDIF
+              ENDDO
+              IF(IPRINT .GE. 100) THEN
+                WRITE(IOUT,6003) IP,IA,IGROUP
+                WRITE(IOUT,6010)
+     >          ((TFLUX(IR,IP,IA)+TFLUX(IR+1,IP,IA))/2.0,IR=1,NREG)
+              ENDIF
+            ENDIF
+          ENDDO
+        ENDDO
+        IF(IFORM .EQ. 1) THEN
+          IF(IKU .EQ. 3) THEN
+            WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ELSE
+            WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ENDIF
+        ELSE
+          IF(IKU .EQ. 3) THEN
+            WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IFPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ELSE
+            WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IFPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ENDIF
+        ENDIF
+      ENDDO
+      IF(IPRINT .GE. 10) THEN
+        FNORM=FLUX(1,1,1,2)/FLUX(1,1,1,1)
+        DO IGROUP=1,NGROUP
+          WRITE(6,6030) IGROUP
+          DO IR=1,NREG
+            RELERR=100.0*(FLUX(IR,1,IGROUP,2)
+     >            -FNORM*FLUX(IR,1,IGROUP,1))
+     >            /FLUX(IR,1,IGROUP,2)
+            ERRMAX=MAX(ERRMAX,ABS(RELERR))
+            IF(IPRINT .GE. 20) THEN
+              WRITE(6,6031) IR,FNORM*FLUX(IR,1,IGROUP,1),
+     >                      FLUX(IR,1,IGROUP,2),RELERR
+            ENDIF
+          ENDDO
+        ENDDO
+        WRITE(6,6020) ERRMAX
+      ENDIF
+*----
+*  Print angular adjoint
+*----
+      IAPU=FILUNIT(KENTRY(3))
+      IAKU=IENTRY(3)
+      ERRMAX=0.0
+      IF(IAKU .EQ. 3) THEN
+        DO IGROUP=NGROUP,1,-1
+          WRITE(IAPU) 'Group',IGROUP
+        ENDDO
+        DO IGROUP=NGROUP+1,NGROUP+6
+          WRITE(IAPU) 'Comments'
+        ENDDO
+      ELSE
+        DO IGROUP=NGROUP,1,-1
+          WRITE(IAPU,'(A8,4X,I8)') 'Groups =',IGROUP
+        ENDDO
+        DO IGROUP=NGROUP+1,NGROUP+6
+          WRITE(IAPU,'(A8)') 'Comments'
+        ENDDO
+      ENDIF
+      DO IGROUP=NGROUP,1,-1
+        DO IA=1,AZMOAQ*2
+          DO IP=1,POLOAQ
+            IF(IFORM .EQ. 1) THEN
+              DO IR=1,NREG
+                TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,2,IGROUP)
+                IF(IA .LE. AZMOAQ) THEN
+                  FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+                ELSE
+                  FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+                ENDIF
+              ENDDO
+              IF(IPRINT .GE. 100) THEN
+                WRITE(IOUT,6004) IP,IA,IGROUP
+                WRITE(IOUT,6010) (TFLUX(IR,IP,IA),IR=1,NREG)
+              ENDIF
+            ELSE
+              IR=1
+              TFLUX(IR,IP,IA)=AFLUX(IR,IP,IA,2,IGROUP)
+              DO IR=1,NREG
+                TFLUX(IR+1,IP,IA)=
+     >            2.0D0*AFLUX(IR,IP,IA,2,IGROUP)-TFLUX(IR,IP,IA)
+                IF(IA .LE. AZMOAQ) THEN
+                  FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(NAZ+1-IA))
+                ELSE
+                  FLUX(IR,2,IGROUP,2)=FLUX(IR,2,IGROUP,2)
+     >           +REAL(AFLUX(IR,IP,IA,2,IGROUP)*WPOL(IP)/WAZI(IA-NAZ))
+                ENDIF
+              ENDDO
+              IF(IPRINT .GE. 100) THEN
+                WRITE(IOUT,6004) IP,IA,IGROUP
+                WRITE(IOUT,6010)
+     >          ((TFLUX(IR,IP,IA)+TFLUX(IR+1,IP,IA))/2.0,IR=1,NREG)
+              ENDIF
+            ENDIF
+          ENDDO
+        ENDDO
+        IF(IFORM .EQ. 1) THEN
+          IF(IAKU .EQ. 3) THEN
+            WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ELSE
+            WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ENDIF
+        ELSE
+          IF(IAKU .EQ. 3) THEN
+            WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IAPU) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ELSE
+            WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=1,AZMOAQ)
+            WRITE(IAPU,1002) (((TFLUX(IR,IP,IA),IR=1,NREG+1),
+     >                   IP=1,POLOAQ),IA=AZMOAQ*2,AZMOAQ+1,-1)
+          ENDIF
+        ENDIF
+      ENDDO
+      IF(IPRINT .GE. 10) THEN
+        FNORM=FLUX(1,2,1,2)/FLUX(1,2,1,1)
+        DO IGROUP=NGROUP,1,-1
+          WRITE(6,6032) IGROUP
+          DO IR=1,NREG
+            RELERR=100.0*(FLUX(IR,2,IGROUP,2)
+     >            -FNORM*FLUX(IR,2,IGROUP,1))
+     >            /FLUX(IR,1,IGROUP,2)
+            ERRMAX=MAX(ERRMAX,ABS(RELERR))
+            IF(IPRINT .GE. 20) THEN
+              WRITE(6,6031) IR,FNORM*FLUX(IR,2,IGROUP,1),
+     >                      FLUX(IR,2,IGROUP,2),RELERR
+            ENDIF
+          ENDDO
+        ENDDO
+        WRITE(6,6020) ERRMAX
+      ENDIF
+*----
+*  Processing finished, return
+*----
+      RETURN
+*----
+*  Formats
+*----
+ 1000 FORMAT('Number of regions =',I5)
+ 1001 FORMAT('Nomber of angles  =',I5)
+ 1002 FORMAT(1P,5E20.10)
+ 1010 FORMAT(A4,1P/(3E15.7))
+ 6000 FORMAT('Output from routine ',A6/
+     >       'Number of groups  =',I5/
+     >       'Number of dimens  =',I5/
+     >       'Number of azimuth =',I5/
+     >       'Number of polar   =',I5/
+     >       'Number of regions =',I5)
+ 6001 FORMAT('Azimuthal quadrature')
+ 6002 FORMAT('Polar quadrature')
+ 6003 FORMAT('Polar =',I5,2X,'Azim  = ',I5,2X,'Group = ',I5/
+     >       'Direct angular flux per region ')
+ 6004 FORMAT('Polar =',I5,2X,'Azim  = ',I5,2X,'Group = ',I5/
+     >       'Adjoint angular flux per region')
+ 6010 FORMAT(1P,5E20.10)
+ 6011 FORMAT(1P,3E20.10)
+ 6020 FORMAT('Maximum relative on flux    (%) = ',F15.7)
+ 6030 FORMAT(' Flux for group    ',I5)
+ 6031 FORMAT(' Region ',I5,1P,2E20.10,5X,0P,F20.10)
+ 6032 FORMAT(' Adjoint for group ',I5)
+      END