Initial commit from Polytechnique Montreal

1 files changed, 179 insertions, 0 deletions

diff --git a/Dragon/src/MPOIDF.f b/Dragon/src/MPOIDF.f
new file mode 100644
index 0000000..5a458e9
--- /dev/null
+++ b/Dragon/src/MPOIDF.f
@@ -0,0 +1,179 @@
+*DECK MPOIDF
+      SUBROUTINE MPOIDF(IPMPO,IPEDIT,HEDIT,NG,NMIL,ICAL,IDF,NALBP,
+     1 FNORM,VOLMIL,FLXMIL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* To store discontinuity factor and albedo information in the MPO file.
+*
+*Copyright:
+* Copyright (C) 2022 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
+* IPMPO   pointer to the MPO file.
+* IPEDIT  pointer to the edition object (L_EDIT signature).
+* HEDIT   name of output group for a (multigroup mesh, output geometry)
+*         couple (generally equal to 'output_0').
+* NG      number of condensed energy groups.
+* NMIL    number of mixtures.
+* ICAL    index of the current elementary calculation.
+* IDF     type of surfacic information (2/3: boundary flux/DF).
+* NALBP   number of physical albedos per energy group.
+* FNORM   flux normalization factor.
+* VOLMIL  mixture volumes.
+* FLXMIL  averaged flux of mixtures.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE hdf5_wrap
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPMPO,IPEDIT
+      INTEGER NG,NMIL,ICAL,IDF,NALBP
+      REAL FNORM,VOLMIL(NMIL),FLXMIL(NMIL,NG)
+      CHARACTER(LEN=12) HEDIT
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER HSMG*131,RECNAM*80
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: SURF
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: VREAL,ALBP
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: DISFAC,ALBP2
+      CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: HADF
+*----
+*  RECOVER DISCONTINUITY FACTOR INFORMATION FROM MACROLIB
+*----
+      CALL LCMSIX(IPEDIT,'MACROLIB',1)
+      CALL LCMLEN(IPEDIT,'ADF',ILONG,ITYLCM)
+      IF(ILONG.NE.0) THEN
+        CALL LCMSIX(IPEDIT,'ADF',1)
+        CALL LCMGET(IPEDIT,'NTYPE',NSURFD)
+        NGG=0
+        IF((IDF.EQ.2).OR.(IDF.EQ.3)) THEN
+          NGG=NG
+        ELSE IF(IDF.EQ.4) THEN
+          NGG=NG*NG
+        ELSE
+          CALL XABORT('MPOIDF: INVALID ADF OPTION.')
+        ENDIF
+        ALLOCATE(DISFAC(NSURFD,NGG,NMIL),SURF(NMIL*NGG),HADF(NSURFD))
+        CALL LCMGTC(IPEDIT,'HADF',8,NSURFD,HADF)
+        DO I=1,NSURFD
+          CALL LCMLEN(IPEDIT,HADF(I),ILONG,ITYLCM)
+          IF(IDF.EQ.2) THEN
+*           boundary flux information
+            IF(ILONG.NE.NMIL*NG) THEN
+              WRITE(HSMG,'(16HMPOIDF: INVALID ,A,8H LENGTH=,I5,
+     1        10H EXPECTED=,I5,4H.(1))') HADF(I),ILONG,NMIL*NG
+              CALL XABORT(HSMG)
+            ENDIF
+            CALL LCMGET(IPEDIT,HADF(I),SURF)
+            DO IMIL=1,NMIL
+              DO IGR=1,NG
+                IF(FNORM.NE.1.0) THEN
+                  DISFAC(I,IGR,IMIL)=SURF((IGR-1)*NMIL+IMIL)*
+     1            FNORM*1.0E13*VOLMIL(IMIL)/FLXMIL(IMIL,IGR)
+                ELSE
+                  DISFAC(I,IGR,IMIL)=SURF((IGR-1)*NMIL+IMIL)*
+     1            VOLMIL(IMIL)/FLXMIL(IMIL,IGR)
+                ENDIF
+              ENDDO
+            ENDDO
+          ELSE IF(IDF.EQ.3) THEN
+*           discontinuity factor information
+            IF(ILONG.NE.NMIL*NG) THEN
+              WRITE(HSMG,'(16HMPOIDF: INVALID ,A,8H LENGTH=,I5,
+     1        10H EXPECTED=,I5,4H.(2))') HADF(I),ILONG,NMIL*NG
+              CALL XABORT(HSMG)
+            ENDIF
+            CALL LCMGET(IPEDIT,HADF(I),SURF)
+            DO IMIL=1,NMIL
+              DO IGR=1,NG
+                IOF=(IGR-1)*NMIL+IMIL
+                DISFAC(I,IGR,IMIL)=SURF(IOF)
+              ENDDO
+            ENDDO
+          ELSE IF(IDF.EQ.4) THEN
+*           matrix discontinuity factor information
+            IF(ILONG.NE.NMIL*NG*NG) THEN
+              WRITE(HSMG,'(16HMPOIDF: INVALID ,A,8H LENGTH=,I5,
+     1        10H EXPECTED=,I5,4H.(3))') HADF(I),ILONG,NMIL*NG*NG
+              CALL XABORT(HSMG)
+            ENDIF
+            CALL LCMGET(IPEDIT,HADF(I),SURF)
+            DO IMIL=1,NMIL
+              DO IGR=1,NG
+                DO JGR=1,NG
+                  IOF=((JGR-1)*NG+IGR-1)*NMIL+IMIL
+                  DISFAC(I,(JGR-1)*NG+IGR,IMIL)=SURF(IOF)
+                ENDDO
+              ENDDO
+            ENDDO
+          ENDIF
+        ENDDO
+        DEALLOCATE(HADF,SURF)
+        CALL LCMSIX(IPEDIT,' ',2)
+*----
+*  MOVE TO THE /statept_id/zone_id/discontinuity GROUP.
+*----
+        DO IMIL=1,NMIL
+          WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,
+     1    15H/discontinuity/)') TRIM(HEDIT),ICAL-1,IMIL-1
+          CALL hdf5_create_group(IPMPO,TRIM(RECNAM))
+          CALL hdf5_write_data(IPMPO,TRIM(RECNAM)//"NSURF",NSURFD)
+          IF((IDF.EQ.2).OR.(IDF.EQ.3)) THEN
+            ALLOCATE(VREAL(NSURFD,NG))
+            VREAL(:NSURFD,:NG)=DISFAC(:NSURFD,:NG,IMIL)
+            CALL hdf5_write_data(IPMPO,TRIM(RECNAM)//"DFACTOR",VREAL)
+          ELSE IF(IDF.EQ.4) THEN
+            ALLOCATE(VREAL(NSURFD,NG*NG))
+            VREAL(:NSURFD,:NG*NG)=DISFAC(:NSURFD,:NG*NG,IMIL)
+            CALL hdf5_write_data(IPMPO,TRIM(RECNAM)//"DFACTORGxG",VREAL)
+          ENDIF
+          DEALLOCATE(VREAL)
+        ENDDO
+        DEALLOCATE(DISFAC)
+      ENDIF
+*----
+*  MOVE TO THE /statept_id/flux GROUP.
+*----
+      IF(NALBP.NE.0) THEN
+        WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/flux/)')
+     1  TRIM(HEDIT),ICAL-1
+        CALL hdf5_create_group(IPMPO,TRIM(RECNAM))
+*----
+*  RECOVER AND SAVE ALBEDO INFORMATION
+*----
+        CALL hdf5_write_data(IPMPO,TRIM(RECNAM)//"NALBP",NALBP)
+        CALL LCMLEN(IPEDIT,'ALBEDO',ILONG,ITYLCM)
+        IF(ILONG.EQ.NALBP*NG) THEN
+*         diagonal physical albedos
+          ALLOCATE(ALBP(NALBP,NG))
+          CALL LCMGET(IPEDIT,'ALBEDO',ALBP)
+          CALL hdf5_write_data(IPMPO,TRIM(RECNAM)//"ALBEDO",ALBP)
+          DEALLOCATE(ALBP)
+        ELSE IF(ILONG.EQ.NALBP*NG*NG) THEN
+*         matrix physical albedos
+          ALLOCATE(ALBP2(NALBP,NG,NG))
+          CALL LCMGET(IPEDIT,'ALBEDO',ALBP2)
+          CALL hdf5_write_data(IPMPO,TRIM(RECNAM)//"ALBEDOGxG",ALBP2)
+          DEALLOCATE(ALBP2)
+        ELSE
+          CALL XABORT('MPOIDF: INCONSISTENT ALBEDO INFORMATION.')
+        ENDIF
+      ENDIF
+      CALL LCMSIX(IPEDIT,' ',2)
+      RETURN
+      END