Initial commit from Polytechnique Montreal

1 files changed, 504 insertions, 0 deletions

diff --git a/Donjon/src/MCRAGF.f b/Donjon/src/MCRAGF.f
new file mode 100644
index 0000000..3c780f6
--- /dev/null
+++ b/Donjon/src/MCRAGF.f
@@ -0,0 +1,504 @@
+*DECK MCRAGF
+      SUBROUTINE MCRAGF(IPMAC,IPMPO,IACCS,NMIL,NMIX,NGRP,NALBP,LALBG,
+     1 LADFM,IMPX,NCAL,TERP,MIXC,NSURFD,HEDIT,VOSAP,VOLMI2,IDF)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Build the macrolib by scanning the NCAL elementary calculations and
+* weighting them with TERP factors. ADF and physical albedos part.
+*
+*Copyright:
+* Copyright (C) 2022 Ecole Polytechnique de Montreal
+*
+*Author(s):
+* A. Hebert
+*
+*Parameters: input
+* IPMAC   address of the output macrolib LCM object.
+* IPMPO   address of the MPO file.
+* IACCS   =0 macrolib is created; =1 ... is updated.
+* NMIL    number of material mixtures in the MPO file.
+* NMIX    maximum number of material mixtures in the macrolib.
+* NGRP    number of energy groups.
+* NALBP   number of physical albedos per energy group.
+* LALBG   type of physical albedos (.true.: diagonal; .false.: GxG).
+* LADFM   type of discontinuity factors (.true.: diagonal; .false.: GxG).
+* IMPX    print parameter (equal to zero for no print).
+* NCAL    number of elementary calculations in the MPO file.
+* TERP    interpolation factors.
+* MIXC    mixture index in the MPO file corresponding to each macrolib
+*         mixture. Equal to zero if a macrolib mixture is not updated.
+* NSURFD  number of discontinuity factors.
+* HEDIT   name of output group for a (multigroup mesh, output geometry)
+*         couple (generally equal to 'output_0').
+* VOSAP   zone volumes in the MPO file.
+* VOLMI2  mixture volumes in the macrolib.
+*
+*Parameters: output
+* IDF     type of discontinuity factors (DF) in the macrolib (=0: not
+*         used; =3: DF used; =4 matrix DF used). 
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE hdf5_wrap
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPMAC,IPMPO
+      INTEGER IACCS,NMIL,NMIX,NGRP,NALBP,IMPX,NCAL,MIXC(NMIX),NSURFD,IDF
+      REAL TERP(NCAL,NMIX),VOSAP(NMIL),VOLMI2(NMIX)
+      LOGICAL LALBG,LADFM
+      CHARACTER(LEN=12) HEDIT
+*----
+*  LOCAL VARIABLES
+*----
+      REAL WEIGHT,FACTOR
+      CHARACTER RECNAM*80,HSMG*131
+      INTEGER IKEFF,IKINF,I,IBM,IBMOLD,ICAL,IGR,JGR,ILONG,ITYLCM,IAL,
+     1 RANK,NBYTE,TYPE,TYPE2,TYPE4,DIMSR(5),NSURFD_OLD,NITMA
+      LOGICAL LNEW
+      DOUBLE PRECISION DGAR1,DGAR2
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: GAR4,ZKINF,ZKEFF,VREAL
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: GAR6,ALBP,AVGFL2,DISFAC,VFLUX
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ADF2,ALBP2,ALBP_ERM,SFLUX
+      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: ADF2M,ALBP2_E
+      CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: HADF
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(GAR4(NGRP*NGRP),GAR6(NGRP,2),ALBP2(NMIX,NALBP,NGRP),
+     1 ALBP2_E(NMIX,NALBP,NGRP,NGRP),ZKINF(NMIX),ZKEFF(NMIX),
+     2 HADF(NSURFD),ADF2(NMIX,NGRP,NSURFD),ADF2M(NMIX,NGRP,NGRP,NSURFD),
+     3 AVGFL2(NMIX,NGRP))
+*----
+*  MICROLIB INITIALIZATION
+*----
+      IKINF=0
+      IKEFF=0
+      IDF=0
+      LNEW=.TRUE.
+      IF(NSURFD.GT.0) THEN
+        WRITE(RECNAM,'(8H/output/,A,32H/statept_0/zone_0/discontinuity/)
+     &  ') TRIM(HEDIT)
+        LNEW=hdf5_group_exists(IPMPO,TRIM(RECNAM))
+        IF(LNEW) THEN
+*         new specification
+          CALL hdf5_info(IPMPO,TRIM(RECNAM)//"DFACTOR",RANK,TYPE2,
+     &    NBYTE,DIMSR)
+          CALL hdf5_info(IPMPO,TRIM(RECNAM)//"DFACTORGxG",RANK,TYPE4,
+     &    NBYTE,DIMSR)
+          IF(TYPE2.NE.99) THEN
+            IDF=3 ! discontinuity factor information
+          ELSE IF(TYPE4.NE.99) THEN
+            IDF=4 ! matrix discontinuity factor information
+          ELSE
+            CALL hdf5_list(IPMPO,TRIM(RECNAM))
+            CALL XABORT('MCRAGF: UNABLE TO SET TYPE OF DF.')
+          ENDIF
+        ELSE
+*         old specification
+          IDF=3 ! discontinuity factor information
+        ENDIF
+        ADF2(:NMIX,:NGRP,:NSURFD)=0.0
+        ADF2M(:NMIX,:NGRP,:NGRP,:NSURFD)=0.0
+      ENDIF
+      AVGFL2(:NMIX,:NGRP)=0.0
+      IF(NALBP.NE.0) ALBP2(:NMIX,:NALBP,:NGRP)=0.0
+      ZKINF(:NMIX)=0.0
+      ZKEFF(:NMIX)=0.0
+      NSURFD_OLD=NSURFD
+      IF(IACCS.NE.0) THEN
+        ! Recover ADF, GFF and physical albedos previously computed
+        CALL LCMLEN(IPMAC,'VOLUME',ILONG,ITYLCM)
+        IF(ILONG.EQ.0) CALL XABORT('MCRAGF: NO VOLUMES IN MACROLIB.')
+        CALL LCMGET(IPMAC,'VOLUME',VOLMI2)
+        IF(NALBP.GT.0) THEN
+          CALL LCMLEN(IPMAC,'ALBEDO',ILONG,ITYLCM)
+          IF(ILONG.EQ.NALBP*NGRP) THEN
+*           diagonal albedo matrix
+            ALLOCATE(ALBP(NALBP,NGRP))
+            CALL LCMGET(IPMAC,'ALBEDO',ALBP)
+            DO IBM=1,NMIX ! mixtures in Macrolib
+              ALBP2(IBM,:NALBP,:NGRP)=ALBP(:NALBP,:NGRP)
+            ENDDO
+            DEALLOCATE(ALBP)
+          ELSE IF(ILONG.EQ.NALBP*NGRP*NGRP) THEN
+*           GxG albedo matrix
+            ALLOCATE(ALBP_ERM(NALBP,NGRP,NGRP))
+            CALL LCMGET(IPMAC,'ALBEDO',ALBP_ERM)
+            DO IBM=1,NMIX ! mixtures in Macrolib
+              ALBP2_E(IBM,:NALBP,:NGRP,:NGRP)=
+     &                  ALBP_ERM(:NALBP,:NGRP,:NGRP)
+            ENDDO
+            DEALLOCATE(ALBP_ERM)
+          ENDIF
+        ENDIF
+        IF(NSURFD_OLD.GT.0) THEN
+          CALL LCMLEN(IPMAC,'ADF',ILONG,ITYLCM)
+          IF(ILONG.EQ.0) THEN
+            CALL LCMLIB(IPMAC)
+            CALL XABORT('MCRAGF: UNABLE TO FIND DIRECTORY ADF.')
+          ENDIF
+          CALL LCMSIX(IPMAC,'ADF',1)
+          CALL LCMGET(IPMAC,'NTYPE',NSURFD)
+          IF(NSURFD.GT.NSURFD_OLD) THEN
+            DEALLOCATE(ADF2M,ADF2,HADF)
+            ALLOCATE(HADF(NSURFD),ADF2(NMIX,NGRP,NSURFD),
+     &      ADF2M(NMIX,NGRP,NGRP,NSURFD))
+            ADF2(:NMIX,:NGRP,:NSURFD)=0.0
+            ADF2M(:NMIX,:NGRP,:NGRP,:NSURFD)=0.0
+          ENDIF
+          CALL LCMGTC(IPMAC,'HADF',8,NSURFD,HADF)
+          DO I=1,NSURFD
+            CALL LCMLEN(IPMAC,HADF(I),ILONG,ITYLCM)
+            IF(ILONG.EQ.NMIX*NGRP) THEN
+              CALL LCMGET(IPMAC,HADF(I),ADF2(1,1,I))
+            ELSE IF(ILONG.EQ.NMIX*NGRP*NGRP) THEN
+              CALL LCMGET(IPMAC,HADF(I),ADF2M(1,1,1,I))
+            ENDIF
+          ENDDO
+          CALL LCMGET(IPMAC,'AVG_FLUX',AVGFL2)
+          CALL LCMSIX(IPMAC,' ',2)
+        ENDIF
+        DO IBM=1,NMIX
+          IF(MIXC(IBM).NE.0) THEN
+            IF(NALBP.NE.0) THEN
+              IF(LALBG) THEN
+                ALBP2(IBM,:NALBP,:NGRP)=0.0
+              ELSE
+                ALBP2_E(IBM,:NALBP,:NGRP,:NGRP)=0.0
+              ENDIF
+            ENDIF
+            IF((NSURFD.GT.0).AND.(IDF.EQ.3)) THEN
+              ADF2(IBM,:NGRP,:NSURFD)=0.0
+            ELSE IF((NSURFD.GT.0).AND.(IDF.EQ.4)) THEN
+              ADF2M(IBM,:NGRP,:NGRP,:NSURFD)=0.0
+            ENDIF
+            AVGFL2(IBM,:NGRP)=0.0
+            ZKINF(IBM)=0.0
+            ZKEFF(IBM)=0.0
+          ENDIF
+        ENDDO
+      ENDIF
+*----
+*  OVERALL ELEMENTARY CALCULATION LOOP
+*----
+      DO 40 ICAL=1,NCAL
+      IF(NSURFD_OLD.GT.0) THEN
+        IF(LNEW) THEN
+          ALLOCATE(SFLUX(NMIL,NGRP**2,NSURFD_OLD),VFLUX(NMIL,NGRP))
+          DO IBMOLD=1,NMIL
+            WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)
+     &      ') TRIM(HEDIT),ICAL-1,IBMOLD-1
+            CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ZONEFLUX",VREAL)
+            VFLUX(IBMOLD,:NGRP)=VREAL(:NGRP)/VOSAP(IBMOLD)
+            DEALLOCATE(VREAL)
+            WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,
+     &      15H/discontinuity/)') TRIM(HEDIT),ICAL-1,IBMOLD-1
+            CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NSURF",NITMA)
+            IF(NITMA.NE.NSURFD_OLD) THEN
+              WRITE(HSMG,'(32HMCRAGF: THE NUMBER OF SURFACES (,I5,
+     &        12H) IN MIXTURE,I5,31H IS DIFFERENT FROM THE NUMBER (,I5,
+     &        15H) IN MIXTURE 1.)') NITMA,IBMOLD,NSURFD_OLD
+              CALL XABORT(HSMG)
+            ENDIF
+            IF(IDF.EQ.3) THEN
+              CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"DFACTOR",DISFAC)
+              DO I=1,NSURFD_OLD
+                SFLUX(IBMOLD,:NGRP,I)=DISFAC(I,:NGRP)
+              ENDDO
+              DEALLOCATE(DISFAC)
+            ELSE IF(IDF.EQ.4) THEN
+              CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"DFACTORGxG",
+     &        DISFAC)
+              DO I=1,NSURFD_OLD
+                SFLUX(IBMOLD,:NGRP**2,I)=DISFAC(I,:NGRP**2)
+              ENDDO
+              DEALLOCATE(DISFAC)
+            ENDIF
+          ENDDO
+        ELSE
+          ALLOCATE(SFLUX(NMIL,NGRP,NSURFD_OLD),VFLUX(NMIL,NGRP))
+          CALL SPHMOL(IPMPO,ICAL,NMIL,NGRP,NSURFD_OLD,HEDIT,VOSAP,
+     &    SFLUX,VFLUX)
+        ENDIF
+        DO IBM=1,NMIX ! mixtures in Macrolib
+          IBMOLD=MIXC(IBM)
+          IF(IBMOLD.GT.NMIL) CALL XABORT('MCRAGF: NMIL OVERFLOW.')
+          IF(IBMOLD.EQ.0) CYCLE
+          WEIGHT=TERP(ICAL,IBM)
+          IF(WEIGHT.EQ.0.0) CYCLE
+          IF(IDF.EQ.3) THEN
+            DO I=1,NSURFD_OLD
+              WRITE(HADF(I),'(3HFD_,I5.5)') I
+              ADF2(IBM,:NGRP,I)=ADF2(IBM,:NGRP,I)+WEIGHT*
+     &        SFLUX(IBMOLD,:NGRP,I)*VFLUX(IBMOLD,:NGRP)
+            ENDDO
+          ELSE IF(IDF.EQ.4) THEN
+            DO I=1,NSURFD_OLD
+              WRITE(HADF(I),'(3HFD_,I5.5)') I
+              DO JGR=1,NGRP
+                DO IGR=1,NGRP
+                  ADF2M(IBM,IGR,JGR,I)=ADF2M(IBM,IGR,JGR,I)+WEIGHT*
+     &            SFLUX(IBMOLD,(JGR-1)*NGRP+IGR,I)*VFLUX(IBMOLD,IGR)
+                ENDDO
+              ENDDO
+            ENDDO
+          ENDIF
+          AVGFL2(IBM,:NGRP)=AVGFL2(IBM,:NGRP)+WEIGHT*VFLUX(IBMOLD,:NGRP)
+        ENDDO
+        DEALLOCATE(VFLUX,SFLUX)
+      ENDIF
+*----
+*  PROCESS PHYSICAL ALBEDO INFORMATION
+*----
+      WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/flux/)')
+     & TRIM(HEDIT),ICAL-1
+      IF((NALBP.GT.0).AND.LALBG) THEN
+*       diagonal albedo matrix
+        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ALBEDO",ALBP)
+        DO 20 IBM=1,NMIX ! mixtures in Macrolib
+        IBMOLD=MIXC(IBM)
+        IF(IBMOLD.EQ.0) GO TO 20
+        WEIGHT=TERP(ICAL,IBM)
+        IF(WEIGHT.EQ.0.0) GO TO 20
+        DO IGR=1,NGRP
+          DO IAL=1,NALBP
+            FACTOR=(1.0-ALBP(IAL,IGR))/(1.0+ALBP(IAL,IGR))
+            ALBP2(IBM,IAL,IGR)=ALBP2(IBM,IAL,IGR)+WEIGHT*FACTOR
+          ENDDO
+        ENDDO
+   20     CONTINUE
+        DEALLOCATE(ALBP)        
+      ELSE IF(NALBP.GT.0) THEN
+*       GxG albedo matrix
+        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ALBEDOGxG",ALBP_ERM)
+        DO 25 IBM=1,NMIX ! mixtures in Macrolib
+        IBMOLD=MIXC(IBM)
+        IF(IBMOLD.EQ.0) GO TO 25
+        WEIGHT=TERP(ICAL,IBM)
+        IF(WEIGHT.EQ.0.0) GO TO 25
+        DO IGR=1,NGRP
+          DO JGR=1,NGRP
+            DO IAL=1,NALBP
+              FACTOR=(1.0-ALBP_ERM(IAL,IGR,JGR))/(1.0+
+     1        ALBP_ERM(IAL,IGR,JGR))
+              ALBP2_E(IBM,IAL,IGR,JGR)=ALBP2_E(IBM,IAL,IGR,JGR)+WEIGHT*
+     1        FACTOR
+            ENDDO
+          ENDDO
+        ENDDO
+   25   CONTINUE
+        DEALLOCATE(ALBP_ERM)        
+      ENDIF
+*----
+*  PROCESS KINF AND KEFF
+*----
+      DO 30 IBM=1,NMIX ! mixtures in Macrolib
+      IBMOLD=MIXC(IBM)
+      IF(IBMOLD.EQ.0) GO TO 30
+      WEIGHT=TERP(ICAL,IBM)
+      IF(WEIGHT.EQ.0.0) GO TO 30
+      WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,8H/addons/)')
+     & TRIM(HEDIT),ICAL-1
+      CALL hdf5_info(IPMPO,TRIM(RECNAM)//"KINF",RANK,TYPE,NBYTE,DIMSR)
+      IF(TYPE.NE.99) THEN
+        IKINF=1
+        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"KINF",FACTOR)
+        ZKINF(IBM)=ZKINF(IBM)+WEIGHT*FACTOR
+      ENDIF
+      CALL hdf5_info(IPMPO,TRIM(RECNAM)//"KEFF",RANK,TYPE,NBYTE,DIMSR)
+      IF(TYPE.NE.99) THEN
+        IKEFF=1
+        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"KEFF",FACTOR)
+        ZKEFF(IBM)=ZKEFF(IBM)+WEIGHT*FACTOR
+      ENDIF
+   30 CONTINUE
+   40 CONTINUE
+*----
+*  SAVE ADF INFORMATION
+*----
+      IF((NSURFD.GT.0).AND.(IDF.EQ.3)) THEN
+        CALL LCMSIX(IPMAC,'ADF',1)
+        CALL LCMPUT(IPMAC,'NTYPE',1,1,NSURFD)
+        CALL LCMPTC(IPMAC,'HADF',8,NSURFD,HADF)
+        DO IBM=1,NMIX
+          IF(MIXC(IBM).EQ.0) CYCLE
+          DO IGR=1,NGRP
+            ADF2(IBM,IGR,:NSURFD)=ADF2(IBM,IGR,:NSURFD)/AVGFL2(IBM,IGR)
+          ENDDO
+          IF(NSURFD.GT.NSURFD_OLD) THEN
+            IF(NSURFD_OLD.GT.2) CALL XABORT('MCRAGF: INVALID NSURFD.')
+*           assign the same ADF to all sides.
+            DO I=2,NSURFD
+              ADF2(IBM,:NGRP,I)=ADF2(IBM,:NGRP,1)
+            ENDDO
+          ENDIF
+        ENDDO
+        IF(LADFM) THEN
+          DO I=1,NSURFD
+            CALL LCMPUT(IPMAC,HADF(I),NMIX*NGRP,2,ADF2(1,1,I))
+          ENDDO
+        ELSE
+*         write non-matrix DF into a matrix DF
+          DO I=1,NSURFD
+            DO IBM=1,NMIX
+              IF(MIXC(IBM).EQ.0) CYCLE
+              ADF2M(IBM,:NGRP,:NGRP,I)=0.0
+              IF(IDF.EQ.3) THEN
+                DO IGR=1,NGRP
+                  ADF2M(IBM,IGR,IGR,I)=ADF2(IBM,IGR,I)
+                ENDDO
+              ELSE IF(IDF.EQ.4) THEN
+                DO IGR=1,NGRP
+                  ADF2M(IBM,IGR,IGR,I)=ADF2(IBM,IGR,I)
+                ENDDO
+              ENDIF
+            ENDDO
+            CALL LCMPUT(IPMAC,HADF(I),NMIX*NGRP*NGRP,2,ADF2M(1,1,1,I))
+          ENDDO
+          IDF=4
+        ENDIF
+        CALL LCMPUT(IPMAC,'AVG_FLUX',NMIX*NGRP,2,AVGFL2)
+        CALL LCMSIX(IPMAC,' ',2)
+        IF(IMPX.GT.1) THEN
+          DO IBM=1,NMIX
+            IF(MIXC(IBM).EQ.0) CYCLE
+            WRITE(6,'(/40H MCRAGF: DISCONTINUITY FACTORS - MIXTURE,I5)')
+     1      IBM
+            DO I=1,NSURFD
+              WRITE(6,'(1X,A,1H:,1P,(5X,10E12.4))') TRIM(HADF(I)),
+     1        (ADF2(IBM,IGR,I)/AVGFL2(IBM,IGR),IGR=1,NGRP)
+            ENDDO
+          ENDDO
+        ENDIF
+      ELSE IF((NSURFD.GT.0).AND.(IDF.EQ.4)) THEN
+        CALL LCMSIX(IPMAC,'ADF',1)
+        CALL LCMPUT(IPMAC,'NTYPE',1,1,NSURFD)
+        CALL LCMPTC(IPMAC,'HADF',8,NSURFD,HADF)
+        DO IBM=1,NMIX
+          IF(MIXC(IBM).EQ.0) CYCLE
+          DO JGR=1,NGRP
+            DO IGR=1,NGRP
+              ADF2M(IBM,IGR,JGR,:NSURFD)=ADF2M(IBM,IGR,JGR,:NSURFD)/
+     1        AVGFL2(IBM,IGR)
+            ENDDO
+          ENDDO
+          IF(NSURFD.GT.NSURFD_OLD) THEN
+            IF(NSURFD_OLD.GT.2) CALL XABORT('MCRAGF: INVALID NSURFD.')
+*           assign the same matrix ADF to all sides.
+            DO I=2,NSURFD
+              ADF2M(IBM,:NGRP,:NGRP,I)=ADF2M(IBM,:NGRP,:NGRP,1)
+            ENDDO
+          ENDIF
+        ENDDO
+        DO I=1,NSURFD
+          CALL LCMPUT(IPMAC,HADF(I),NMIX*NGRP*NGRP,2,ADF2M(1,1,1,I))
+        ENDDO
+        CALL LCMPUT(IPMAC,'AVG_FLUX',NMIX*NGRP,2,AVGFL2)
+        CALL LCMSIX(IPMAC,' ',2)
+        IF(IMPX.GT.1) THEN
+          DO IBM=1,NMIX
+            IF(MIXC(IBM).EQ.0) CYCLE
+            WRITE(6,'(/44H MCRAGF: MATRIX DISCONTINUITY FACTORS - MIXT,
+     1      3HURE,I5)') IBM
+            DO I=1,NSURFD
+              WRITE(6,'(1X,A,1H:,1P,(5X,10E12.4))') TRIM(HADF(I)),
+     1        ((ADF2M(IBM,IGR,JGR,I),IGR=1,NGRP),JGR=1,NGRP)
+            ENDDO
+          ENDDO
+        ENDIF
+      ENDIF
+*----
+*  AVERAGE PHYSICAL ALBEDO INFORMATION
+*----
+      IF((NALBP.GT.0).AND.LALBG) THEN
+*       diagonal albedo matrix
+        ALLOCATE(ALBP(NALBP,NGRP))
+        DO IGR=1,NGRP
+          DO IAL=1,NALBP
+            DGAR1=0.0D0
+            DGAR2=0.0D0
+            DO IBM=1,NMIX
+              IF(VOLMI2(IBM).EQ.0.0) CYCLE
+              DGAR1=DGAR1+ALBP2(IBM,IAL,IGR)*VOLMI2(IBM)
+              DGAR2=DGAR2+VOLMI2(IBM)
+            ENDDO
+            ALBP(IAL,IGR)=REAL((1.D0-DGAR1/DGAR2)/(1.D0+DGAR1/DGAR2))
+          ENDDO
+        ENDDO
+        IF(LADFM) THEN
+          CALL LCMPUT(IPMAC,'ALBEDO',NALBP*NGRP,2,ALBP)
+        ELSE
+*         write non-matrix albedo into a matrix albedo
+          ALLOCATE(ALBP_ERM(NALBP,NGRP,NGRP))
+          ALBP_ERM(:NALBP,:NGRP,:NGRP)=0.0
+          DO IGR=1,NGRP
+            DO IAL=1,NALBP
+              ALBP_ERM(IAL,IGR,IGR)=ALBP(IAL,IGR)
+            ENDDO
+          ENDDO
+          CALL LCMPUT(IPMAC,'ALBEDO',NALBP*NGRP*NGRP,2,ALBP_ERM)
+          DEALLOCATE(ALBP_ERM)
+        ENDIF
+        DEALLOCATE(ALBP)
+      ELSE IF(NALBP.GT.0) THEN
+*       GxG albedo matrix
+        ALLOCATE(ALBP_ERM(NALBP,NGRP,NGRP))
+        DO IGR=1,NGRP
+          DO JGR=1,NGRP
+            DO IAL=1,NALBP
+              DGAR1=0.0D0
+              DGAR2=0.0D0
+              DO IBM=1,NMIX
+                IF(VOLMI2(IBM).EQ.0.0) CYCLE
+                DGAR1=DGAR1+ALBP2_E(IBM,IAL,IGR,JGR)*VOLMI2(IBM)
+                DGAR2=DGAR2+VOLMI2(IBM)
+              ENDDO
+              ALBP_ERM(IAL,IGR,JGR)=REAL((1.D0-DGAR1/DGAR2)/(1.D0+
+     1        DGAR1/DGAR2))
+            ENDDO
+          ENDDO
+        ENDDO
+        CALL LCMPUT(IPMAC,'ALBEDO',NALBP*NGRP*NGRP,2,ALBP_ERM)
+        DEALLOCATE(ALBP_ERM)
+      ENDIF
+*----
+*  AVERAGE KINF
+*----
+      IF(IKINF.EQ.1) THEN
+        DGAR1=0.0D0
+        DGAR2=0.0D0
+        DO IBM=1,NMIX
+          DGAR1=DGAR1+ZKINF(IBM)*VOLMI2(IBM)
+          DGAR2=DGAR2+VOLMI2(IBM)
+        ENDDO
+        FACTOR=REAL(DGAR1/DGAR2)
+        CALL LCMPUT(IPMAC,'K-INFINITY',1,2,FACTOR)
+      ENDIF
+*----
+*  AVERAGE KEFF
+*----
+      IF(IKEFF.EQ.1) THEN
+        DGAR1=0.0D0
+        DGAR2=0.0D0
+        DO IBM=1,NMIX
+          DGAR1=DGAR1+ZKEFF(IBM)*VOLMI2(IBM)
+          DGAR2=DGAR2+VOLMI2(IBM)
+        ENDDO
+        FACTOR=REAL(DGAR1/DGAR2)
+        CALL LCMPUT(IPMAC,'K-EFFECTIVE',1,2,FACTOR)
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(AVGFL2,ADF2M,ADF2,HADF,ZKEFF,ZKINF,ALBP2_E,ALBP2,
+     1 GAR6,GAR4)
+      RETURN
+      END