Initial commit from Polytechnique Montreal

1 files changed, 288 insertions, 0 deletions

diff --git a/Trivac/src/ERRDRV.f b/Trivac/src/ERRDRV.f
new file mode 100755
index 0000000..b5ac203
--- /dev/null
+++ b/Trivac/src/ERRDRV.f
@@ -0,0 +1,288 @@
+*DECK ERRDRV
+      SUBROUTINE ERRDRV(IMPX,IPMAC1,IPMAC2,NREG,NREG2,NGRP,HREAC,ERAMAX,
+     1 ERASUM,ERQMAX,ERQSUM,ERGMARR,ERGSARR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform reaction rate statistics between two extended macrolibs.
+*
+*Copyright:
+* Copyright (C) 2002 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
+* IPMAC1  pointer to the reference extended macrolib.
+* IPMAC2  pointer to the approximate extended macrolib.
+* NREG    number of regions in the macrolib.
+* NREG2   number of regions used for statistics.
+* NGRP    number of energy groups in the macrolib.
+* HREAC   nuclear reaction used to compute power map
+*
+*Parameters: output
+* ERAMAX  maximum relative error on absorption rates.
+* ERASUM  average relative error on absorption rates.
+* ERQMAX  maximum relative error on QUANDRY powers.
+* ERQSUM  average relative error on QUANDRY powers.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPMAC1,IPMAC2
+      INTEGER IMPX,NREG,NGRP
+      REAL ERAMAX,ERASUM,ERQMAX,ERQSUM,ERGMARR(NGRP),ERGSARR(NGRP)
+      CHARACTER HREAC*8
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (NSTATE=40)
+      CHARACTER HSMG*131
+      INTEGER IDATA(NSTATE)
+      TYPE(C_PTR) JPMAC1,KPMAC1,JPMAC2,KPMAC2
+      REAL, DIMENSION(:), ALLOCATABLE :: VOL1,VOL2,TOTAL,GAR,FLUX,
+     1 QUAN1,QUAN2,TRABS1,TRABS2
+      REAL, DIMENSION(:,:), ALLOCATABLE :: TRA1,TRA2,XABS1,XABS2
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(TRA1(NREG2,NGRP),TRA2(NREG,NGRP),XABS1(NREG,NGRP),
+     1 XABS2(NREG,NGRP),VOL1(NREG),VOL2(NREG),TOTAL(NREG),GAR(NREG),
+     2 FLUX(NREG),QUAN1(NREG),QUAN2(NREG),TRABS1(NREG),TRABS2(NREG))
+*----
+*  RECOVER REFERENCE REACTION RATES:
+*----
+      CALL LCMGET(IPMAC1,'STATE-VECTOR',IDATA)
+      IF((NREG.NE.IDATA(2)).OR.(NGRP.NE.IDATA(1))) THEN
+         CALL XABORT('ERRDRV: INVALID VALUE OF NREG OR NGRP.')
+      ENDIF
+      CALL LCMGET(IPMAC1,'VOLUME',VOL1)
+      VOL1T=0.0
+      PWR1T=0.0
+      DO 10 I=1,NREG2
+      TRABS1(I)=0.0
+      QUAN1(I)=0.0
+      VOL1T=VOL1T+VOL1(I)
+   10 CONTINUE
+      TRA1(:NREG2,:NGRP)=0.0
+      CALL ERRABS(IPMAC1,NREG2,NREG,NGRP,XABS1)
+      JPMAC1=LCMGID(IPMAC1,'GROUP')
+      DO 35 IGR=1,NGRP
+      KPMAC1=LCMGIL(JPMAC1,IGR)
+      CALL LCMGET(KPMAC1,'NTOT0',TOTAL)
+      CALL LCMGET(KPMAC1,'SIGW00',GAR)
+      CALL LCMGET(KPMAC1,'FLUX-INTG',FLUX)
+      DO 20 I=1,NREG2
+      IF(VOL1(I).EQ.0.0) GO TO 20
+      TRA1(I,IGR)=(TOTAL(I)-GAR(I))*FLUX(I)/VOL1(I)
+      TRABS1(I)=TRABS1(I)+XABS1(I,IGR)*FLUX(I)/VOL1(I)
+   20 CONTINUE
+      CALL LCMLEN(KPMAC1,HREAC,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) THEN
+         WRITE(HSMG,'(32HERRDRV: UNABLE TO FIND REACTION ,A,1H.)') HREAC
+         CALL XABORT(HSMG)
+      ENDIF
+      CALL LCMGET(KPMAC1,HREAC,GAR)
+      DO 30 I=1,NREG2
+      QUAN1(I)=QUAN1(I)+GAR(I)*FLUX(I)
+      PWR1T=PWR1T+QUAN1(I)
+   30 CONTINUE
+   35 CONTINUE
+*----
+*  RECOVER APPROXIMATE REACTION RATES:
+*----
+      CALL LCMGET(IPMAC2,'STATE-VECTOR',IDATA)
+      IF((NREG.NE.IDATA(2)).OR.(NGRP.NE.IDATA(1))) THEN
+         CALL XABORT('ERRDRV: INVALID VALUE OF NREG OR NGRP.')
+      ENDIF
+      CALL LCMGET(IPMAC2,'VOLUME',VOL2)
+      VOL2T=0.0
+      PWR2T=0.0
+      DO 50 I=1,NREG2
+      TRABS2(I)=0.0
+      QUAN2(I)=0.0
+      VOL2T=VOL2T+VOL2(I)
+   50 CONTINUE
+      CALL ERRABS(IPMAC2,NREG2,NREG,NGRP,XABS2)
+      JPMAC2=LCMGID(IPMAC2,'GROUP')
+      DO 80 IGR=1,NGRP
+      KPMAC2=LCMGIL(JPMAC2,IGR)
+      CALL LCMGET(KPMAC2,'NTOT0',TOTAL)
+      CALL LCMGET(KPMAC2,'SIGW00',GAR)
+      CALL LCMGET(KPMAC2,'FLUX-INTG',FLUX)
+      DO 60 I=1,NREG2
+      IF(VOL2(I).EQ.0.0) GO TO 60
+      TRA2(I,IGR)=(TOTAL(I)-GAR(I))*FLUX(I)/VOL2(I)
+      TRABS2(I)=TRABS2(I)+XABS2(I,IGR)*FLUX(I)/VOL2(I)
+   60 CONTINUE
+      IF(ILONG.NE.0) THEN
+         CALL LCMGET(KPMAC2,HREAC,GAR)
+         DO 70 I=1,NREG2
+         QUAN2(I)=QUAN2(I)+GAR(I)*FLUX(I)
+         PWR2T=PWR2T+QUAN2(I)
+   70    CONTINUE
+      ENDIF
+   80 CONTINUE
+*----
+*  COMPUTE QUANDRY TYPE NORMALIZED POWER DENSITIES.
+*----
+      IF(ILONG.GT.0) THEN
+         DO 90 I=1,NREG2
+         IF(VOL1(I).NE.0.0) QUAN1(I)=QUAN1(I)/VOL1(I)
+         IF(VOL2(I).NE.0.0) QUAN2(I)=QUAN2(I)/VOL2(I)
+         IF(PWR1T.NE.0.0) QUAN1(I)=QUAN1(I)*VOL1T/PWR1T
+         IF(PWR2T.NE.0.0) QUAN2(I)=QUAN2(I)*VOL2T/PWR2T
+   90    CONTINUE
+      ENDIF
+*----
+*  PRINT STATISTICS ON GROUPWISE REMOVAL RATES.
+*----
+      WRITE(6,'(/47H ERRDRV: STATISTICS ON GROUPWISE REMOVAL RATES:)')
+      SUMREF=0.0
+      SUM=0.0
+      DO 125 IGR=1,NGRP
+      DO 120 I=1,NREG2
+      SUMREF=SUMREF+TRA1(I,IGR)*VOL1(I)
+      SUM=SUM+TRA2(I,IGR)*VOL2(I)
+  120 CONTINUE
+  125 CONTINUE
+      DO 150 IGR=1,NGRP
+      WRITE (6,'(/17H PROCESSING GROUP,I3)') IGR
+      ERGMAX=0.0
+      ERGSUM=0.0
+      VOLTOT=0.0
+      DO 130 I=1,NREG2
+      TRA2(I,IGR)=TRA2(I,IGR)*(SUMREF/SUM)*(VOL2T/VOL1T)
+      IF(TRA1(I,IGR).NE.0.0) THEN
+         VOLTOT=VOLTOT+VOL1(I)
+         GAR(I)=100.0*(TRA2(I,IGR)-TRA1(I,IGR))/TRA1(I,IGR)
+      ELSE
+         GAR(I)=0.0
+      ENDIF
+      ERGSUM=ERGSUM+VOL1(I)*ABS(GAR(I))
+      ERGMAX=MAX(ERGMAX,ABS(GAR(I)))
+  130 CONTINUE
+      ERGSUM=ERGSUM/VOLTOT
+      ERGMARR(IGR)=ERGMAX
+      ERGSARR(IGR)=ERGSUM
+      IF(IMPX.GT.1) WRITE (6,'(/8X,9HREFERENCE,7X,6HAPPROX,7X,5HERROR)')
+      DO 140 I=1,NREG2
+      IF(IMPX.GT.1) WRITE (6,'(4X,I4,1X,1P,2E13.5,0P,F9.3,2H %)')
+     1 I,TRA1(I,IGR),TRA2(I,IGR),GAR(I)
+  140 CONTINUE
+      WRITE(6,300) IGR, ERGMAX,ERGMAX,ERGMAX
+      WRITE(6,310) IGR, ERGSUM,ERGSUM,ERGSUM
+  150 CONTINUE
+      WRITE(6,400) MAXVAL(ERGMARR), MAXVAL(ERGMARR),
+     1   MAXVAL(ERGMARR)
+      WRITE(6,410) MAXVAL(ERGSARR), MAXVAL(ERGSARR),
+     1   MAXVAL(ERGSARR)
+*----
+*  PRINT STATISTICS ON CONDENSED ABSORPTION RATES.
+*----
+      WRITE(6,'(/40H ERRDRV: STATISTICS ON ABSORPTION RATES:)')
+      SUMREF=0.0
+      SUM=0.0
+      DO 160 I=1,NREG2
+      SUMREF=SUMREF+TRABS1(I)*VOL1(I)
+      SUM=SUM+TRABS2(I)*VOL2(I)
+  160 CONTINUE
+      ERAMAX=0.0
+      ERASUM=0.0
+      VOLTOT=0.0
+      DO 165 I=1,NREG2
+      TRABS2(I)=TRABS2(I)*(SUMREF/SUM)*(VOL2T/VOL1T)
+      IF(TRABS1(I).NE.0.0) THEN
+         VOLTOT=VOLTOT+VOL1(I)
+         GAR(I)=100.0*(TRABS2(I)-TRABS1(I))/TRABS1(I)
+      ELSE
+         GAR(I)=0.0
+      ENDIF
+      ERASUM=ERASUM+VOL1(I)*ABS(GAR(I))
+      ERAMAX=MAX(ERAMAX,ABS(GAR(I)))
+  165 CONTINUE
+      ERASUM=ERASUM/VOLTOT
+      IF(IMPX.GT.1) WRITE (6,'(/8X,9HREFERENCE,7X,6HAPPROX,7X,5HERROR)')
+      DO 170 I=1,NREG2
+      IF(IMPX.GT.1) WRITE (6,'(4X,I4,1X,1P,2E13.5,0P,F9.3,2H %)')
+     1 I,TRABS1(I),TRABS2(I),GAR(I)
+  170 CONTINUE
+      WRITE(6,420) ERAMAX,ERAMAX,ERAMAX
+      WRITE(6,430) ERASUM,ERASUM,ERASUM
+*----
+*  PRINT STATISTICS ON QUANDRY TYPE NORMALIZED POWER DENSITIES.
+*----
+      IF(ILONG.NE.0) THEN
+         WRITE(6,'(/48H ERRDRV: STATISTICS ON QUANDRY TYPE NORMALIZED P,
+     1   15HOWER DENSITIES:)')
+         ERQMAX=0.0
+         ERQSUM=0.0
+         VOLTOT=0.0
+         DO 180 I=1,NREG2
+         ERR=ABS(VOL1(I)/VOL1T-VOL2(I)/VOL2T)
+         IF(ERR.GT.1.0E-4*ABS(VOL1(I)/VOL1T)) THEN
+            WRITE(HSMG,'(37HERRDRV: INCONSISTENT VOLUME IN REGION,I5,
+     1      3H BY,F7.2,2H %)') I,ERR*100.0
+            CALL XABORT(HSMG)
+         ENDIF
+         GAR(I)=0.0
+         IF(QUAN1(I).EQ.0.0) GO TO 180
+         VOLTOT=VOLTOT+VOL1(I)
+         GAR(I)=100.0*(QUAN2(I)-QUAN1(I))/QUAN1(I)
+         ERQSUM=ERQSUM+VOL1(I)*ABS(QUAN1(I)-QUAN2(I))/QUAN1(I)
+         ERQMAX=MAX(ERQMAX,ABS(GAR(I)))
+  180    CONTINUE
+         IF(VOLTOT.NE.0.0) ERQSUM=100.0*ERQSUM/VOLTOT
+         IF(IMPX.GT.1)
+     1      WRITE(6,'(/8X,9HREFERENCE,7X,6HAPPROX,7X,5HERROR)')
+         DO 190 I=1,NREG2
+         IF((QUAN1(I).NE.0.0).OR.(QUAN2(I).NE.0.0)) THEN
+            IF(IMPX.GT.1) WRITE(6,'(4X,I4,1X,1P,2E13.5,0P,F9.3,2H %)')
+     1      I,QUAN1(I),QUAN2(I),GAR(I)
+         ENDIF
+  190    CONTINUE
+         WRITE(6,440) ERQMAX,ERQMAX,ERQMAX
+         WRITE(6,450) ERQSUM,ERQSUM,ERQSUM
+      ENDIF
+*----
+*  PRINT STATISTICS ON K-EFFECTIVE.
+*----
+      CALL LCMLEN(IPMAC1,'K-EFFECTIVE',LENGT,ITYLCM)
+      IF(LENGT.EQ.1) THEN
+         CALL LCMGET(IPMAC1,'K-EFFECTIVE',FKEFF1)
+         CALL LCMGET(IPMAC2,'K-EFFECTIVE',FKEFF2)
+         WRITE(6,'(/5X,22HREFERENCE K-EFFECTIVE=,F9.6/8X,11HAPPROX K-EF,
+     1   8HFECTIVE=,F9.6,8H  ERROR=,F9.1,4H PCM)') FKEFF1,FKEFF2,
+     2   (FKEFF2-FKEFF1)*1.0E5
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(XABS2,XABS1,TRA1,TRA2,VOL1,VOL2,TOTAL,GAR,FLUX,QUAN2,
+     1 QUAN1,TRABS2,TRABS1)
+      RETURN
+*
+  300 FORMAT(/1X,37HGROUPWISE REM. RATE MAX ERR FOR GROUP,I4,2H =,
+     1   F9.3,2H %,F9.2,2H %,F9.1,2H %)
+  310 FORMAT( 1X,37HGROUPWISE REM. RATE AV. ERR FOR GROUP,I4,2H =,
+     1   F9.3,2H %,F9.2,2H %,F9.1,2H %/)
+  400 FORMAT(/1X,30HMAXIMUM ERROR OVER ALL GROUPS=,F9.3,2H %,F9.2,2H %,
+     1   F9.1,2H %)
+  410 FORMAT( 1X,30HAVERAGE ERROR OVER ALL GROUPS=,F9.3,2H %,F9.2,2H %,
+     1   F9.1,2H %/)
+  420 FORMAT(/1X,30HABSORPTION RATE MAXIMUM ERROR=,F9.3,2H %,F9.2,2H %,
+     1   F9.1,2H %)
+  430 FORMAT( 1X,30HABSORPTION RATE AVERAGE ERROR=,F9.3,2H %,F9.2,2H %,
+     1   F9.1,2H %/)
+  440 FORMAT(/1X,28HPOWER DENSITY MAXIMUM ERROR=,F9.3,2H %,F9.2,2H %,
+     1   F9.1,2H %)
+  450 FORMAT( 1X,28HPOWER DENSITY AVERAGE ERROR=,F9.3,2H %,F9.2,2H %,
+     1   F9.1,2H %/)
+      END