Initial commit from Polytechnique Montreal

1 files changed, 435 insertions, 0 deletions

diff --git a/Dragon/src/SHISN2.f b/Dragon/src/SHISN2.f
new file mode 100644
index 0000000..edd4fee
--- /dev/null
+++ b/Dragon/src/SHISN2.f
@@ -0,0 +1,435 @@
+*DECK SHISN2
+      SUBROUTINE SHISN2 (IPLIB,IPTRK,IFTRAK,NGRO,NBISO,NBM,NREG,NUN,
+     1 CDOOR,NRES,NBM2,IMPX,ISONAM,MIX,DEN,SN,SB,LSHI,IPHASE,MAT,VOL,
+     2 KEYFLX,LEAKSW,TITR,START,SIGT,SIGT3,NOCONV,BIEFF,LGC,SIGE)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform one multidimensional self-shielding iteration using the
+* generalized Stamm'ler algorithm without Nordheim (PIC) approximation.
+*
+*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
+* IPLIB   pointer to the internal microscopic cross section library
+*         (L_LIBRARY signature).
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IFTRAK  unit number of the sequential binary tracking file.
+* NGRO    number of energy groups.
+* NBISO   number of isotopes present in the calculation domain.
+* NBM     number of mixtures in the macrolib.
+* NREG    number of volumes.
+* NUN     number of unknowns in the flux or source vector in one
+*         energy group.
+* CDOOR   name of the geometry/solution module.
+* NRES    number of resonant mixtures.
+* NBM2    number of resonant isotopes.
+* IMPX    print flag.
+* ISONAM  alias name of isotopes.
+* MIX     mix number of each isotope (can be zero).
+* DEN     density of each isotope.
+* LSHI    resonant region number associated with each isotope.
+*         Infinite dilution will be assumed if LSHI(i)=0.
+* IPHASE  type of flux solution (=1 use a native flux solution door;
+*         =2 use collision probabilities).
+* MAT     index-number of the mixture type assigned to each volume.
+* VOL     volumes.
+* KEYFLX  pointers of fluxes in unknown vector.
+* LEAKSW  leakage flag (.TRUE. only if leakage is present on the outer
+*         surface).
+* TITR    title.
+* START   beginning-of-iteration flag (.TRUE. if SHISN2 is called
+*         for the first time).
+* SIGT3   transport correction.
+* NOCONV  mixture convergence flag (.TRUE. if mixture IBM
+*         is not converged in group L).
+* BIEFF   Livolant-Jeanpierre normalization flag (.TRUE. to
+*         activate).
+* LGC     Goldstein-Cohen approximation flag (.TRUE. to activate).
+*
+*Parameters: output
+* SN      on input, estimate of the dilution cross section in each 
+*         energy group of each isotope. A value of 1.0e10 is used 
+*         for infinite dilution.
+*         On output, computed dilution cross section in each energy 
+*         group of each isotope.
+* SIGT    total macroscopic cross sections as modified by Shiba.
+*
+*Parameters: output
+* SB      dilution cross section as used in Livolant-Jeanpierre
+*         normalization.
+* SIGE    computed macroscopic dilution cross section in each resonant
+*         mixture and each energy group.
+*
+*Reference:
+* A. Hebert and G. Marleau, Generalization of the Stamm'ler Method
+* for the Self-Shielding of Resonant isotopes in Arbitrary Geometries,
+* Nucl. Sci. Eng. 108, 230 (1991).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      PARAMETER (NALPHA=5)
+      TYPE(C_PTR) IPLIB,IPTRK
+      INTEGER IFTRAK,NGRO,NBISO,NBM,NREG,NUN,NRES,NBM2,IMPX,
+     1 ISONAM(3,NBISO),MIX(NBISO),LSHI(NBISO),IPHASE,MAT(NREG),
+     2 KEYFLX(NREG)
+      REAL DEN(NBISO),SN(NGRO,NBISO),SB(NGRO,NBISO),VOL(NREG),
+     1 SIGT(NBM,NGRO),SIGT3(NBM,NGRO),SIGE(NRES,NGRO)
+      CHARACTER CDOOR*12,TITR*72,CGRPNM*12
+      LOGICAL LEAKSW,START,NOCONV(NBM,NGRO),BIEFF,LGC
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER HSMG*131
+      LOGICAL LOGDO
+      COMPLEX COEF(3),DENOM(3),EAV
+      PARAMETER (NRAT=(NALPHA+1)/2)
+      TYPE(C_PTR) KPLIB
+      REAL FACT(NALPHA),SIGX(NALPHA)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IRES,MIX2,IRNBM,NPSYS
+      REAL, ALLOCATABLE, DIMENSION(:) :: GAR,SIGRES,DILAV,FUN
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: SIG0,SIG1,SIG3,TOTAL,SIGOLD,
+     1 DILUT
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKI
+      TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO
+*----
+*  DATA STATEMENTS
+*----
+      DATA FACT/0.01,0.1,1.0,10.0,100.0/
+*----
+*  SCRATCH STORAGE ALLOCATION
+*   SIG0    macroscopic xs of the resonant isotopes as interpolated.
+*   SIG1    macroscopic xs of the resonant isotopes at various SIGX.
+*   SIG3    macroscopic transport correction.
+*----
+      ALLOCATE(IRES(NBM),MIX2(NBISO),IRNBM(NBM),NPSYS(NGRO))
+      ALLOCATE(SIG0(NBM,NGRO),SIG1(NBM,NGRO),SIG3(NBM,NGRO),
+     1 TOTAL(NGRO,NBM2),SIGOLD(NGRO,NBM2),GAR(NGRO),SIGRES(NBM),
+     2 DILAV(NGRO),DILUT(NALPHA,NGRO))
+      ALLOCATE(MASKI(NBISO))
+      ALLOCATE(IPISO(NBISO))
+*----
+*  SET THE LCM MICROLIB ISOTOPEWISE DIRECTORIES.
+*----
+      CALL LIBIPS(IPLIB,NBISO,IPISO)
+*----
+*  UNLOAD MICROSCOPIC X-S FROM LCM TO SCRATCH STORAGE
+*----
+      IBM=0
+      DO 20 ISO=1,NBISO
+      MIX2(ISO)=0
+      IF(LSHI(ISO).GT.0) THEN
+         IBM=IBM+1
+         MIX2(ISO)=IBM
+         KPLIB=IPISO(ISO) ! set ISO-th isotope
+         CALL LCMGET(KPLIB,'NTOT0',TOTAL(1,IBM))
+         CALL LCMLEN(KPLIB,'NGOLD',LENGT,ITYLCM)
+         IF((LENGT.EQ.NGRO).AND.(.NOT.START).AND.LGC) THEN
+            IF(IMPX.GE.5) WRITE (6,390) (ISONAM(I0,ISO),I0=1,3)
+            CALL LCMGET(KPLIB,'SIGS00',SIGOLD(1,IBM))
+            CALL LCMGET(KPLIB,'NGOLD',GAR)
+            DO 10 LLL=1,NGRO
+            SIGOLD(LLL,IBM)=(1.0-GAR(LLL))*SIGOLD(LLL,IBM)
+   10       CONTINUE
+         ELSE
+            SIGOLD(:NGRO,IBM)=0.0
+         ENDIF
+      ENDIF
+   20 CONTINUE
+*----
+*  LOOP OVER RESONANT REGIONS. THE CP ARE STORED ON DIRECTORY SHIBA
+*----
+      CALL LCMSIX(IPLIB,'SHIBA',1)
+      DO 260 INRS=1,NRES
+*----
+*  FIND THE RESONANT MIXTURE NUMBERS (IRNBM) ASSOCIATED WITH REGION INRS
+*----
+      NBNRS=0
+      DO 50 IBM=1,NBM
+      IRES(IBM)=0
+      DO 40 ISO=1,NBISO
+      IF((MIX(ISO).EQ.IBM).AND.(LSHI(ISO).EQ.INRS)) THEN
+         NBNRS=NBNRS+1
+         IRNBM(NBNRS)=IBM
+         IRES(IBM)=1
+         GO TO 50
+      ENDIF
+   40 CONTINUE
+   50 CONTINUE
+      IF(NBNRS.EQ.0) THEN
+         IF(START.AND.(IMPX.GE.1)) WRITE(6,385) 'SHISN2',INRS
+         GO TO 260
+      ELSE IF(START.AND.(NBNRS.GT.1).AND.(IMPX.GE.5)) THEN
+         WRITE (6,380) NBNRS,INRS
+      ENDIF
+*
+      NPSYS(:NGRO)=0
+      DO 120 LLL=1,NGRO
+      LOGDO=.FALSE.
+      DO 60 I=1,NBNRS
+      LOGDO=LOGDO.OR.NOCONV(IRNBM(I),LLL)
+   60 CONTINUE
+      IF(LOGDO) THEN
+         NPSYS(LLL)=LLL
+*
+*        COMPUTE THE LIGHT AND RESONANT COMPONENTS OF THE MACROSCOPIC
+*        CROSS SECTIONS IN EACH RESONANT MIXTURE.
+         DO 80 I=1,NBNRS
+         SIGRES(I)=0.0
+         DO 70 ISO=1,NBISO
+         IF((MIX(ISO).EQ.IRNBM(I)).AND.(LSHI(ISO).EQ.INRS)) THEN
+            SIGRES(I)=SIGRES(I)+TOTAL(LLL,MIX2(ISO))*DEN(ISO)
+         ENDIF
+   70    CONTINUE
+         SIGT(IRNBM(I),LLL)=SIGT(IRNBM(I),LLL)-SIGRES(I)
+   80    CONTINUE
+         DO 90 IBM=1,NBM
+         SIG0(IBM,LLL)=0.0
+         SIG1(IBM,LLL)=0.0
+         SIG3(IBM,LLL)=SIGT3(IBM,LLL)
+   90    CONTINUE
+         DO 110 I=1,NBNRS
+         SIG0(IRNBM(I),LLL)=SIGRES(I)
+         SIG3(IRNBM(I),LLL)=0.0
+  110    CONTINUE
+         IF(IMPX.GE.10) THEN
+            WRITE (6,400) LLL,(SIG0(I,LLL),I=1,NBM)
+            WRITE (6,410) LLL,(SIGT(I,LLL),I=1,NBM)
+            WRITE (6,420) LLL,(SIGT3(I,LLL),I=1,NBM)
+         ENDIF
+      ENDIF
+  120 CONTINUE
+*----
+*  SET UP VECTORS DILUT AND SIGX.
+*----
+      DILAV(:NGRO)=0.0
+      IF(START) THEN
+*        USE A VERY CHEAP APPROXIMATION TO START ITERATIONS.
+         ALLOCATE(FUN(NUN*NGRO))
+         CALL LCMSIX(IPLIB,'--AVERAGE--',1)
+         CALL SHIDST (IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG,
+     1   NUN,NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG0,SIGT,
+     2   SIGT3(1,1),TITR,FUN,DILAV)
+         CALL LCMSIX(IPLIB,' ',2)
+         DEALLOCATE(FUN)
+         DO 135 LLL=1,NGRO
+         DO 130 IALP=1,NALPHA
+         DILUT(IALP,LLL)=DILAV(LLL)
+  130    CONTINUE
+  135    CONTINUE
+      ELSE
+         DO 165 IALP=1,NALPHA
+         DO 150 LLL=1,NGRO
+         IF(NPSYS(LLL).NE.0) THEN
+            DO 140 I=1,NBNRS
+            SIG1(IRNBM(I),LLL)=FACT(IALP)*SIGE(INRS,LLL)
+  140       CONTINUE
+         ENDIF
+  150    CONTINUE
+         ALLOCATE(FUN(NUN*NGRO))
+         WRITE(CGRPNM,'(8H--BAND--,I4.4)') IALP
+         CALL LCMSIX(IPLIB,CGRPNM,1)
+         CALL SHIDST (IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG,
+     1   NUN,NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG1,SIGT,
+     2   SIG3(1,1),TITR,FUN,DILAV)
+         CALL LCMSIX(IPLIB,' ',2)
+         DEALLOCATE(FUN)
+         DO 160 LLL=1,NGRO
+         DILUT(IALP,LLL)=DILAV(LLL)
+  160    CONTINUE
+  165    CONTINUE
+      ENDIF
+*----
+*  COMPUTE AVERAGE MACROSCOPIC DILUTION X-S (SIGE) USING A THREE-TERM
+*  RATIONAL APPROXIMATION.
+*----
+      DO 200 LLL=1,NGRO
+      IF(NPSYS(LLL).NE.0) THEN
+         DO 170 IALP=1,NALPHA
+         SIGX(IALP)=FACT(IALP)*SIGE(INRS,LLL)
+  170    CONTINUE
+         IMPX2=IMPX
+         IF(START) IMPX2=MAX(0,IMPX-10)
+*        **********************************************************
+         CALL SHIRAT(IMPX2,NRAT,SIGX,DILUT(1,LLL),LLL,A,COEF,DENOM)
+*        **********************************************************
+         EAV=(COEF(1)*SQRT(DENOM(1))+COEF(2)*SQRT(DENOM(2))+
+     1   COEF(3)*SQRT(DENOM(3)))**2
+         SIGE(INRS,LLL)=REAL(EAV)
+         IF((.NOT.START).AND.(BIEFF).AND.(NBNRS.EQ.1)) THEN
+*           COMPUTE DILAV FOR THE L-J NORMALIZATION.
+            SIGXX=SIG0(IRNBM(1),LLL)
+            PXX=REAL(COEF(1)/(SIGXX+DENOM(1))+COEF(2)/(SIGXX+DENOM(2))
+     1      +COEF(3)/(SIGXX+DENOM(3)))
+            DILAV(LLL)=1.0/PXX-SIGXX
+         ENDIF
+*----
+*  COMPUTE THE ISOTOPE DILUTION MICROSCOPIC CROSS SECTIONS (SN) USED
+*  FOR LIBRARY INTERPOLATION.
+*----
+         DO 190 ISO=1,NBISO
+         IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND.
+     1   (DEN(ISO).NE.0.)) THEN
+            SUM=0.0
+            DO 180 JSO=1,NBISO
+            IBM=MIX(JSO)
+            IF((LSHI(JSO).EQ.INRS).AND.(IBM.EQ.MIX(ISO)).AND.
+     1      (ISO.NE.JSO)) SUM=SUM+(TOTAL(LLL,MIX2(JSO))-
+     2      SIGOLD(LLL,MIX2(JSO)))*DEN(JSO)
+  180       CONTINUE
+            SN(LLL,ISO)=REAL(COEF(1)*SQRT(DENOM(1)+SUM)+COEF(2)*
+     1      SQRT(DENOM(2)+SUM)+COEF(3)*SQRT(DENOM(3)+SUM))**2/DEN(ISO)
+            IF(SN(LLL,ISO).LE.0.0) THEN
+               WRITE (HSMG,510) (ISONAM(I0,ISO),I0=1,3),SN(LLL,ISO),LLL
+               CALL XABORT(HSMG)
+            ENDIF
+         ELSE IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND.
+     1   (DEN(ISO).EQ.0.)) THEN
+            SN(LLL,ISO)=1.0E10
+         ENDIF
+  190    CONTINUE
+         IF((.NOT.START).AND.(IMPX.GE.10)) THEN
+            DO 195 I=1,NBNRS
+            PP=A-SIGT(IRNBM(I),LLL)
+            QQ=SIGE(INRS,LLL)-SIGT(IRNBM(I),LLL)
+            IF(ABS(PP).GT.1.0E-4*SIGT(IRNBM(I),LLL)) THEN
+               BEL=QQ/PP
+            ELSE
+               BEL=0.0
+            ENDIF
+            WRITE (6,610) I,SIGE(INRS,LLL),BEL
+  195       CONTINUE
+         ENDIF
+      ENDIF
+  200 CONTINUE
+*----
+*  COMPUTE THE ISOTOPE DILUTION MICROSCOPIC CROSS SECTIONS (SB) USED
+*  IN L-J NORMALIZATION.
+*----
+      IF((.NOT.START).AND.(BIEFF).AND.(NBNRS.GT.1)) THEN
+*        COMPUTE DILAV FOR THE L-J NORMALIZATION.
+         ALLOCATE(FUN(NUN*NGRO))
+         CALL LCMSIX(IPLIB,'--AVERAGE--',1)
+         CALL SHIDST (IPLIB,NPSYS,IPTRK,IFTRAK,CDOOR,IMPX,NBM,NREG,
+     1   NUN,NGRO,IPHASE,MAT,VOL,KEYFLX,LEAKSW,IRES,SIG0,SIGT,
+     2   SIGT3(1,1),TITR,FUN,DILAV)
+         CALL LCMSIX(IPLIB,' ',2)
+         DEALLOCATE(FUN)
+      ENDIF
+      DO 250 LLL=1,NGRO
+      IF(NPSYS(LLL).NE.0) THEN
+         DO 220 ISO=1,NBISO
+         IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND.
+     1   (DEN(ISO).NE.0.)) THEN
+            SUM=0.0
+            DO 210 JSO=1,NBISO
+            IBM=MIX(JSO)
+            IF((LSHI(JSO).EQ.INRS).AND.(IBM.EQ.MIX(ISO)).AND.
+     1      (ISO.NE.JSO)) SUM=SUM+TOTAL(LLL,MIX2(JSO))*DEN(JSO)
+  210       CONTINUE
+            IF(START.OR.(.NOT.BIEFF)) THEN
+               SB(LLL,ISO)=SN(LLL,ISO)
+            ELSE
+               SB(LLL,ISO)=(DILAV(LLL)+SUM)/DEN(ISO)
+               IF(SB(LLL,ISO).LT.0.0) THEN
+                  WRITE (HSMG,515) (ISONAM(I0,ISO),I0=1,3),SB(LLL,ISO),
+     1            LLL
+                  CALL XABORT(HSMG)
+               ELSE IF(SB(LLL,ISO).LT.SN(LLL,ISO)) THEN
+                  IF(SB(LLL,ISO).LT.0.99*SN(LLL,ISO)) WRITE (6,520)
+     1            (ISONAM(I0,ISO),I0=1,3),SB(LLL,ISO)/SN(LLL,ISO),LLL
+                  SB(LLL,ISO)=SN(LLL,ISO)
+               ENDIF
+            ENDIF
+         ELSE IF((LSHI(ISO).EQ.INRS).AND.(IRES(MIX(ISO)).EQ.1).AND.
+     1   (DEN(ISO).EQ.0.)) THEN
+            SB(LLL,ISO)=1.0E10
+         ENDIF
+  220    CONTINUE
+*
+*        RESTORE SIGT ARRAY.
+         DO 240 I=1,NBNRS
+         SIGRES(I)=0.0
+         DO 230 ISO=1,NBISO
+         IF((MIX(ISO).EQ.IRNBM(I)).AND.(LSHI(ISO).EQ.INRS)) THEN
+            SIGRES(I)=SIGRES(I)+TOTAL(LLL,MIX2(ISO))*DEN(ISO)
+         ENDIF
+  230    CONTINUE
+         SIGT(IRNBM(I),LLL)=SIGT(IRNBM(I),LLL)+SIGRES(I)
+  240    CONTINUE
+      ENDIF
+  250 CONTINUE
+  260 CONTINUE
+      CALL LCMSIX(IPLIB,' ',2)
+*----
+*  SAVE THE GROUP- AND ISOTOPE-DEPENDENT DILUTIONS
+*----
+      CALL LCMPUT(IPLIB,'ISOTOPESDSB',NBISO*NGRO,2,SB)
+      CALL LCMPUT(IPLIB,'ISOTOPESDSN',NBISO*NGRO,2,SN)
+*----
+*  COMPUTE THE SELF-SHIELDED MICROSCOPIC CROSS SECTIONS AND UPDATE
+*  VECTOR SIGT
+*----
+      DO 290 ISO=1,NBISO
+      LOGDO=START.OR.(DEN(ISO).NE.0.)
+      MASKI(ISO)=(LSHI(ISO).GT.0).AND.LOGDO
+  290 CONTINUE
+      IMPX2=MAX(0,IMPX-1)
+      CALL LIBLIB (IPLIB,NBISO,MASKI,IMPX2)
+      DO 320 ISO=1,NBISO
+      IBM=MIX(ISO)
+      IF((LSHI(ISO).GT.0).AND.(IBM.GT.0).AND.(DEN(ISO).NE.0.)) THEN
+         KPLIB=IPISO(ISO) ! set ISO-th isotope
+         CALL LCMGET(KPLIB,'NTOT0',GAR)
+         DO 300 LLL=1,NGRO
+         TOTAL(LLL,MIX2(ISO))=TOTAL(LLL,MIX2(ISO))-GAR(LLL)
+  300    CONTINUE
+         DO 310 LLL=1,NGRO
+         IF(NOCONV(IBM,LLL)) SIGT(IBM,LLL)=SIGT(IBM,LLL)-DEN(ISO)*
+     1   TOTAL(LLL,MIX2(ISO))
+  310    CONTINUE
+      ENDIF
+  320 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IPISO)
+      DEALLOCATE(MASKI)
+      DEALLOCATE(DILUT,DILAV,SIGRES,GAR,SIGOLD,TOTAL,SIG3,SIG1,SIG0)
+      DEALLOCATE(NPSYS,IRNBM,MIX2,IRES)
+      RETURN
+*
+  380 FORMAT(/16H SHISN2: MERGING,I3,30H RESONANT MIXTURES IN RESONANT,
+     1 14H REGION NUMBER,I3,1H.)
+  385 FORMAT(A6,1X,': RESONANT REGION =',I10,1X,'NOT USED.')
+  390 FORMAT(/53H SHISN2: GOLDSTEIN AND COHEN APPROXIMATION USED FOR I,
+     1 8HSOTOPE ',3A4,2H'.)
+  400 FORMAT(1X,'TOTAL MACROSCOPIC CROSS SECTIONS OF THE RESONANT ',
+     1'MATERIALS IN EACH MIXTURE (GROUP',I5,'):'/(1X,1P,11E11.3))
+  410 FORMAT(1X,'TOTAL MACROSCOPIC CROSS SECTIONS OF THE OTHER ',
+     1'MATERIALS IN EACH MIXTURE (GROUP',I5,'):'/(1X,1P,11E11.3))
+  420 FORMAT(//1X,'TRANSPORT CORRECTION CROSS SECTIONS OF THE OTHER ',
+     1'MATERIALS IN EACH MIXTURE (GROUP',I5,'):'/(1X,1P,11E11.3))
+  510 FORMAT(30HSHISN2: THE RESONANT ISOTOPE ',3A4,14H' HAS A NEGATI,
+     1 27HVE DILUTION CROSS-SECTION (,1P,E14.4,0P,10H) IN GROUP,I4,1H.)
+  515 FORMAT(30HSHISN2: THE RESONANT ISOTOPE ',3A4,14H' HAS A NEGATI,
+     1 22HVE L-J CROSS-SECTION (,1P,E14.4,0P,10H) IN GROUP,I4,1H.)
+  520 FORMAT(54H SHISN2: THE L-J EQUIVALENCE FACTOR OF RESONANT ISOTOP,
+     1 3HE ',3A4,18H' WAS CHANGED FROM,F6.3,16H TO 1.0 IN GROUP,I4,1H.)
+  610 FORMAT(8X,8HAVERAGE(,I2,1H),1P,E13.5/8X,11HBELL FACTOR,E13.5)
+      END