path: root/Dragon/src/FLU2DR.f

*DECK FLU2DR
      SUBROUTINE FLU2DR(IPRT,IPMACR,IPFLUX,IPSYS,IPTRK,IPFLUP,IPSOU,
     1 IGPT,IFTRAK,CXDOOR,TITLE,NUNKNO,NREG,NSOUT,NANIS,NLF,NLIN,NFUNL,
     2 NGRP,NMAT,NIFIS,LFORW,LEAKSW,MAXINR,EPSINR,MAXOUT,EPSUNK,EPSOUT,
     3 NCPTL,NCPTA,ITYPEC,IPHASE,ITPIJ,ILEAK,OPTION,REFKEF,MATCOD,
     4 KEYFLX,VOL,XSTRC,XSDIA,XSNUF,XSCHI,LREBAL,INITFL,NMERG,IMERG)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Fixed source problem or inverse power method for K-effective or
* buckling iteration. Perform thermal iterations.
*
*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): R. Roy
*
*Parameters: input
* IPRT    print flag.
* IPMACR  pointer to the macrolib LCM object.
* IPFLUX  pointer to the flux LCM object.
* IPSYS   pointer to the system LCM object.
* IPTRK   pointer to the tracking LCM object.
* IPFLUP  pointer to the unperturbed flux LCM object (if ITYPEC=1).
* IPSOU   pointer to the fixed source LCM object (if ITYPEC=0 or 1).
* IGPT    index of the fixed source eigenvalue problem to solve.
* IFTRAK  tracking file unit number.
* CXDOOR  character name of the flux solution door.
* TITLE   title.
* NUNKNO  number of unknowns per energy group including spherical
*         harmonic terms, interface currents and fundamental
*         currents.
* NREG    number of regions.
* NSOUT   number of outer surfaces.
* NANIS   maximum cross section Legendre order in object IPMACR.
* NLF     number of Legendre orders for the flux.
* NLIN    number of polynomial components in flux spatial expansion.
* NFUNL   number of spherical harmonics components.
* NGRP    number of energy groups.
* NMAT    number of mixtures in the macrolib.
* NIFIS   number of fissile isotopes.
* LFORW   flag set to .false. to solve an adjoint problem.
* LEAKSW  leakage flag (=.true. if leakage is present on the outer
*         surface).
* MAXINR  maximum number of thermal iterations.
* EPSINR  thermal iterations epsilon.
* MAXOUT  maximum number of outer iterations.
* EPSUNK  outer iterations eigenvector epsilon.
* EPSOUT  outer iterations eigenvalue epsilon.
* NCPTL   number of free iterations in an acceleration cycle.
* NCPTA   number of accelerated iterations in an acceleration cycle.
* ITYPEC  type of flux evaluation:
*         =-2 Fourier analysis;
*         =-1 skip the flux calculation;
*         =0 fixed sources;
*         =1 fixed source eigenvalue problem (GPT type);
*         =2 fission sources/k effective convergence;
*         =3 fission sources/k effective convergence/
*             db2 buckling evaluation;
*         =4 fission sources/db2 buckling convergence;
*         =5 b2 sources/db2 buckling convergence.
* IPHASE  type of flux solution door (1 for asm 2 for pij).
* ITPIJ   type of cp available:
*         =1 scatt mod pij (wij);
*         =2 stand. pij;
*         =3 scatt mod pij+pijk (wij,wijk);
*         =4 standard pij+pijk.
* ILEAK   method used to include DB2 effect:
*         =1 the scattering modified cp matrix is multiplied by PNLR;
*         =2 the reduced cp matrix is multiplied by PNL;
*         =3 sigs0-db2 approximation;
*         =4 albedo approximation;
*         =5 Todorova-type isotropic streaming model;
*         =6 Ecco-type isotropic streaming model;
*         >6 Tibere type anisotropic streaming model.
* OPTION  type of leakage coefficients:
*         'LKRD' (recover leakage coefficients in Macrolib);
*         'RHS' (recover leakage coefficients in RHS flux object);
*         'B0' (B-0), 'P0' (P-0), 'B1' (B-1),
*         'P1' (P-1), 'B0TR' (B-0 with transport correction) or 'P0TR'
*         (P-0 with transport correction).
* REFKEF  target effective multiplication factor.
* MATCOD  mixture indices.
* KEYFLX  index of L-th order flux components in unknown vector.
* VOL     volumes.
* XSTRC   transport-corrected macroscopic total cross sections.
* XSDIA   transport-corrected macroscopic within-group scattering cross
*         sections.
* XSNUF   nu*macroscopic fission cross sections.
* XSCHI   fission spectrum.
* LREBAL  thermal iteration rebalancing flag (=.true. if thermal
*         rebalancing required).
* INITFL  flux initialization flag (=0/1/2: uniform flux/LCM/DSA).
* NMERG   number of leakage zones.
* IMERG   leakage zone index in each material mixture zone.
*
*-----------------------------------------------------------------------
*
*----
*  INTERNAL PARAMETERS:
*   SYBILF : SYBIL FLUX SOLUTION DOOR                 EXT ROUTINE
*   TRFICF : DEFAULT CP FLUX SOLUTION DOOR            EXT ROUTINE
*   BIVAF  : DEFAULT 2D DIFFUSION FLUX SOLUTION DOOR  EXT ROUTINE
*   TRIVAF : DEFAULT 3D DIFFUSION FLUX SOLUTION DOOR  EXT ROUTINE
*   PNF    : DEFAULT PN/SPN FLUX SOLUTION DOOR        EXT ROUTINE
*   SNF    : DEFAULT SN FLUX SOLUTION DOOR            EXT ROUTINE
*----
*
      USE GANLIB
      USE DOORS_MOD
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPMACR,IPFLUX,IPSYS,IPTRK,IPFLUP,IPSOU
      INTEGER IPRT,IGPT,IFTRAK,NUNKNO,NREG,NSOUT,NANIS,NLF,NLIN,NFUNL,
     1 NGRP,NMAT,NIFIS,MAXINR,MAXOUT,NCPTL,NCPTA,ITYPEC,IPHASE,ITPIJ,
     2 ILEAK,MATCOD(NREG),KEYFLX(NREG,NLIN,NFUNL),INITFL,NMERG,
     3 IMERG(NMAT)
      REAL EPSINR,EPSUNK,EPSOUT,VOL(NREG),XSTRC(0:NMAT,NGRP),
     1 XSDIA(0:NMAT,0:NANIS,NGRP),XSNUF(0:NMAT,NIFIS,NGRP),
     2 XSCHI(0:NMAT,NIFIS,NGRP)
      CHARACTER CXDOOR*12,TITLE*72,OPTION*4,HLEAK*6
      LOGICAL LFORW,LEAKSW,LREBAL,CFLI,CEXE
      DOUBLE PRECISION REFKEF
*----
*  LOCAL VARIABLES
*----
      PARAMETER(NSTATE=40,PI=3.141592654)
      TYPE(C_PTR) IPREB,J1,JPSOU,JPFLUX,JPMACR,KPMACR,JPSYS,KPSYS,IPSTR,
     1 JPSTR,KPSTR,JPFLUP1,JPFLUP2,JPSOUR
      INTEGER JPAR(NSTATE),KEYSPN(NREG)
      CHARACTER CAN(0:19)*2,MESSIN*8,MESSOU*5,HTYPE(0:5)*4
      INTEGER INDD(3)
      DOUBLE PRECISION AKEEP(8),FISOUR,OLDBIL,AKEFF,AKEFFO,AFLNOR,
     1 BFLNOR,DDELN1,DDELD1,PROD,FLXIN
      LOGICAL LSCAL,LEXAC,REBFLG
      REAL ALBEDO(6),FLUXC(NREG),B2(4)
*
************************************************************************
*                                                                      *
*   ICHAR       : COUNTER FOR NUM. OF OUTER ITERATIONS                 *
*   ICTOT       : TOTAL NUMBER OF FLUX CALCULATIONS                    *
*                                                                      *
************************************************************************
*
*----
*  ALLOCATABLE ARRAYS
*----
      INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS,NPSYS,KEYCUR,
     1 MATALB
      REAL, ALLOCATABLE, DIMENSION(:) :: DHOM,FXSOR,XSCAT,GAMMA,V,FL,DFL
      REAL, ALLOCATABLE, DIMENSION(:,:) :: DIFHET,SFNU
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: FLUX
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: XCSOU
*----
*  FOR NUMERICAL FOURIER ANALYSIS
*----
      REAL OMEGA,XLEN,EKKK,EVALRHO,SPECR
      REAL, ALLOCATABLE, DIMENSION(:) :: ARRAYRHO1
      REAL, ALLOCATABLE, DIMENSION(:) :: XXX
*----
*  DATA STATEMENTS
*----
      SAVE CAN,HTYPE
      DATA CAN /'00','01','02','03','04','05','06','07','08','09',
     >          '10','11','12','13','14','15','16','17','18','19'/
      DATA (HTYPE(JJ),JJ=0,5)/'S   ','P   ',2*'K   ','B   ','L   '/
*----
*  SCRATCH STORAGE ALLOCATION
*   DHOM    homogeneous leakage coefficients.
*   DIFHET  heterogeneous leakage coefficients.
*   FLUX    iteration flux:
*           FLUX(:,:,1) <=old     outer;
*           FLUX(:,:,2) <=present outer;
*           FLUX(:,:,3) <=new     outer;
*           FLUX(:,:,4) <=source  outer;
*           FLUX(:,:,5) <=old     inner;
*           FLUX(:,:,6) <=present inner;
*           FLUX(:,:,7) <=new     inner;
*           FLUX(:,:,8) <=source  inner.
*----
      ALLOCATE(IJJ(0:NMAT),NJJ(0:NMAT),IPOS(0:NMAT),NPSYS(NGRP))
      ALLOCATE(FLUX(NUNKNO,NGRP,8),XSCAT(0:NMAT*NGRP),GAMMA(NGRP),
     1 DHOM(NGRP),DIFHET(NMERG,NGRP),XCSOU(NGRP))
*
      REBFLG=.TRUE.
      IPREB=IPMACR
*
      AKEEP(:8)=0.0D0
      ICHAR=0
      ICTOT=0
*----
*  RECOVER INDEX FOR THE CURRENTS IN FLUX, NUMERICAL SURFACES,
*  ALBEDO IF NEEDED BY THE REBALANCING.
*----
      ICREB=0
      NNN=0
      INSB=0
      IBFP=0
      NDIM=0
      NFOU=0
      LX=0
      ITYPE=0
      IF(CXDOOR.EQ.'MCCG') THEN
         CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR)
         NANIS_TRK=JPAR(6)
         NDIM=JPAR(16)
         INSB=JPAR(22)
         CALL LCMLEN(IPTRK,'KEYCUR$MCCG',ICREB,ITYLCM)
         IF(ICREB.GT.0) THEN
            CALL LCMLEN(IPTRK,'NZON$MCCG',ILONG,ITYLCM)
            NNN=ILONG-ICREB
            ALLOCATE(KEYCUR(NSOUT),V(ILONG),MATALB(ILONG))
            CALL LCMGET(IPTRK,'KEYCUR$MCCG',KEYCUR)
            CALL LCMGET(IPTRK,'V$MCCG',V)
            CALL LCMGET(IPTRK,'ALBEDO',ALBEDO)
            CALL LCMGET(IPTRK,'NZON$MCCG',MATALB)
         ENDIF
      ELSE IF(CXDOOR.EQ.'SN') THEN
         CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR)
         ITYPE=JPAR(6)
         NDIM=JPAR(9)
         LX=JPAR(12)
         LY=JPAR(13)
         INSB=JPAR(27)
         IBFP=JPAR(31)
         NFOU=JPAR(34)
      ENDIF
*----
*  SELECT THE CALCULATION DOORS FOR WHICH A GROUP-BY-GROUP SCALAR
*  PROCEDURE WILL BE USED. A VECTORIAL APPROACH WILL BE USED WITH
*  OTHER DOORS. LSCAL is true for CXDOOR = 'TRAFIC'.
*----
      LSCAL=(INSB.EQ.0)
*
      CALL KDRCPU(CPU0)
      IF(ILEAK.LT.6) THEN
        INORM=1
      ELSE IF(ILEAK.EQ.6) THEN
        INORM=2 ! Ecco
      ELSE IF(ILEAK.GE.7) THEN
        INORM=3 ! Tibere
      ENDIF
      LEXAC=.FALSE.
      AKEEP(5)=1.0D0
      AKEEP(6)=1.0D0
      AKEEP(7)=1.0D0
      DIFHET(:NMERG,:NGRP)=0.0
      GAMMA(:NGRP)=1.0
*----
*  EXTERNAL FLUX(:,:,2) INITIALISATION AND FIXED-EXTERNAL SOURCE IN
*  FLUX(:,:,4)
*----
      IF(ITYPEC.GE.3) THEN
         CALL LCMGET(IPFLUX,'B2  B1HOM',B2(4))
         IF(ILEAK.GE.7) CALL LCMGET(IPFLUX,'B2  HETE',B2)
      ELSE
         B2(:4)=0.0
      ENDIF
      AKEFFO=0.0D0
      IF(LFORW) THEN
         IF(ITYPEC.EQ.1) THEN
           J1=LCMGID(IPSOU,'DSOUR')
           JPSOU=LCMGIL(J1,IGPT)
           J1=LCMGID(IPFLUX,'DFLUX')
           JPFLUX=LCMGIL(J1,IGPT)
         ELSE
           IF(C_ASSOCIATED(IPSOU)) THEN
             J1=LCMGID(IPSOU,'DSOUR')
             JPSOU=LCMGIL(J1,1)
           ENDIF
           JPFLUX=LCMGID(IPFLUX,'FLUX')
           JPSOUR=LCMLID(IPFLUX,'SOUR',NGRP)
         ENDIF
      ELSE
         IF(ITYPEC.EQ.1) THEN
           J1=LCMGID(IPSOU,'ASOUR')
           JPSOU=LCMGIL(J1,IGPT)
           J1=LCMGID(IPFLUX,'ADFLUX')
           JPFLUX=LCMGIL(J1,IGPT)
         ELSE
           IF(C_ASSOCIATED(IPSOU)) THEN
             J1=LCMGID(IPSOU,'ASOUR')
             JPSOU=LCMGIL(J1,1)
           ENDIF
           JPFLUX=LCMGID(IPFLUX,'AFLUX')
           JPSOUR=LCMLID(IPFLUX,'SOUR',NGRP)
         ENDIF
      ENDIF
      ALLOCATE(FXSOR(0:NMAT))
      JPMACR=LCMGID(IPMACR,'GROUP')
      DO 20 IG=1,NGRP
      FLUX(:NUNKNO,IG,2)=0.0
      FLUX(:NUNKNO,IG,4)=0.0
      CALL LCMLEL(JPFLUX,1,ILINIT,ITYLCM)
      IF(LFORW) THEN
         CALL LCMGDL(JPFLUX,IG,FLUX(1,IG,2))
      ELSE
         CALL LCMGDL(JPFLUX,NGRP-IG+1,FLUX(1,IG,2))
      ENDIF
*
      IF(((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)).AND.
     1 (.NOT.C_ASSOCIATED(IPSOU))) THEN
         KPMACR=LCMGIL(JPMACR,IG)
         FXSOR(0)=0.0
         CALL LCMGET(KPMACR,'FIXE',FXSOR(1))
         CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,FXSOR,
     >   FLUX(1,IG,4))
      ELSE IF(((ITYPEC.EQ.0).OR.(ITYPEC.EQ.1).OR.(ITYPEC.EQ.-2))
     1 .AND.C_ASSOCIATED(IPSOU))THEN
         IF(LFORW) THEN
            CALL LCMGDL(JPSOU,IG,FLUX(1,IG,4))
         ELSE
            CALL LCMGDL(JPSOU,NGRP-IG+1,FLUX(1,IG,4))
         ENDIF
      ENDIF
   20 CONTINUE
      DEALLOCATE(FXSOR)
*-------
*  IF IMPORTED FLUX PRESENT FOR SN, REORDER FLUX.
*-------
      IF((CXDOOR.EQ.'SN').AND.(INITFL.EQ.2)) THEN
         CALL LCMLEN(IPFLUX,'KEYFLX',ILINIT,ITYLCM)
         IF(ILINIT.NE.NREG) THEN
            WRITE(*,*) NREG, ILINIT
            CALL XABORT('FLU2DR: NUMBER OF REGIONS FROM SPN CALCULATION'
     1      //' (OBTAINED FROM LENGTH OF KEYFLX) DOES NOT MATCH NUMBER '
     2      //'OF REGIONS IN SN CALCULATION. CHECK INPUT FILE FOR '
     3      //'POTENTIAL ERRORS.')
         ENDIF
         KEYSPN(:) = 0
         CALL LCMGET(IPFLUX,'KEYFLX',KEYSPN)
         DO 25 IG=1,NGRP
            CALL SNEST(IPTRK,IPRT,NREG,NUNKNO,MATCOD,IG,KEYFLX,KEYSPN,
     1      FLUX(:,IG,2))
   25    CONTINUE
      ENDIF
*----
*  FOURIER ANALYSIS FLUX INITIALISATION
*----
      ALLOCATE(ARRAYRHO1(NFOU**NDIM))
      ARRAYRHO1(:)=0.0
      IFACOUNT=-1
   26 IFACOUNT=IFACOUNT+1
*
      IF(ITYPEC.EQ.-2) THEN
         IF(NFOU.EQ.0)
     >      CALL XABORT('FLU2DR: NEED TO SPECIFY FOURIER ANALYSIS '
     >      //'KEYWORD NFOU IN TRACKING, AS WELL AS NUMBER OF '
     >      //'FREQUENCIES TO INVESTIGATE.')
         IF(CXDOOR.NE.'SN')
     >      CALL XABORT('FLU2DR: FOURIER ANALYSIS ONLY MEANT FOR SN')
         IF(NGRP.NE.1)
     >      CALL XABORT('FLU2DR: FOURIER ANALYSIS NOT MEANT FOR MULTI-'
     >      //'GROUP PROBLEMS. CONSIDER ADDING THAT FUNCTIONALITY. ')
         IGR=1

         FLUX(:,:,:) = 0.0

         SUMXSTRC=0.0
         DO IR=1,NREG
            SUMXSTRC = SUMXSTRC + XSTRC(MATCOD(IR),1)
         ENDDO
         AVXSTRC = (SUMXSTRC/NREG)

         ALLOCATE(XXX(LX+1))
         CALL LCMGET(IPTRK,'XXX',XXX)
         CALL LCMGET(IPTRK,'XLEN',XLEN)

         OMEGA = (2*PI)/(XLEN*AVXSTRC)
         EKKK  =(REAL(IFACOUNT)/NFOU)

         IF(ITYPE.EQ.2)THEN
            PARTX=0.0
            DO IX=1,LX
               PARTX = XXX(IX) + ((XXX(IX+1)-XXX(IX))/2)
               IND=KEYFLX(IX,1,1)
               IF(IND.GT.0)
     >            FLUX(IND,IGR,2) = COS(EKKK*OMEGA*AVXSTRC*PARTX)
            ENDDO
         ELSE
            CALL XABORT('FLU2DR: FOURIER ANALYSIS FOR GEOMETRIES OTHER '
     >         //'THAN CARTESIAN 1D NOT AVAILABLE.')
         ENDIF
         DEALLOCATE(XXX)
      ENDIF
*----
*  COMPUTE FIRST K-EFFECTIVE
*----
      IF((ITYPEC.EQ.0).OR.(ITYPEC.EQ.5).OR.(ITYPEC.EQ.-2)) THEN
         AKEFFO=1.0D0
         AKEFF=1.0D0
         AFLNOR=1.0D0
      ELSE IF(ITYPEC.EQ.1) THEN
         CALL LCMGET(IPFLUP,'STATE-VECTOR',JPAR)
         IF(JPAR(6).GE.2) THEN
            CALL LCMGET(IPFLUP,'K-EFFECTIVE',RKEFF)
            CALL LCMGET(IPFLUP,'K-INFINITY',CUREIN)
         ENDIF
         AKEFF=RKEFF
         IF(JPAR(6).GE.3) THEN
            B2(:4)=0.0
            CALL LCMGET(IPFLUP,'B2  B1HOM',B2(4))
         ENDIF
         IF((JPAR(6).GT.2).AND.(JPAR(7).GE.6)) THEN
            CALL LCMGET(IPFLUP,'B2  HETE',B2)
         ENDIF
         IF((JPAR(6).GT.2).AND.(JPAR(7).GE.5)) THEN
            CALL LCMGET(IPFLUP,'GAMMA',GAMMA)
         ENDIF
         AKEFFO=AKEFF
         AKEEP(2)=AKEFF
         AFLNOR=1.0D0/RKEFF
      ELSE
         OLDBIL=0.0D0
         CALL FLUKEF(IPRT,IPMACR,NGRP,NREG,NUNKNO,NMAT,NIFIS,NANIS,
     1   MATCOD(1),VOL,KEYFLX(1,1,1),XSTRC,XSDIA,XSNUF,XSCHI,NMERG,
     2   IMERG,DIFHET,FLUX(1,1,2),B2,ILEAK,LEAKSW,OLDBIL,AKEFF,AFLNOR)
         AKEFFO=AKEFF
         AKEEP(2)=AKEFF
      ENDIF
      B2VALO=B2(4)
*
      NCTOT=NCPTA+NCPTL
      IF(NCPTA.EQ.0) THEN
         NCPTM=NCTOT+1
      ELSE
         NCPTM=NCPTL
      ENDIF
      MESSOU='     '
      IF(IPRT.GT.0) WRITE(6,1090) 0,1.0,EPSOUT,AKEFFO,B2(4)
*----
*  CALCULATION OF THE INITIAL LEAKAGE COEFFICIENTS
*----
      IF(ITYPEC.GT.2) THEN
         DIFHET(:NMERG,:NGRP)=0.0
         IF(OPTION.EQ.'LKRD') THEN
            CALL LCMGET(IPMACR,'STATE-VECTOR',JPAR)
            IF(JPAR(2).NE.NMAT) THEN
               CALL XABORT('FLU2DR: INVALID NMAT IN THE MACROLIB.')
            ELSE IF(JPAR(9).NE.1) THEN
               CALL XABORT('FLU2DR: INVALID LEAKAGE IN THE MACROLIB.')
            ENDIF
            ALLOCATE(FL(NMAT),DFL(NMAT))
            JPMACR=LCMGID(IPMACR,'GROUP')
            DO IG=1,NGRP
               KPMACR=LCMGIL(JPMACR,IG)
               CALL LCMLEN(KPMACR,'DIFF',ILONG,ITYLCM)
               IF(ILONG.EQ.0) CALL XABORT('FLU2DR: UNABLE TO RECOVER T'
     >         //'HE DIFF RECORD IN THE MACROLIB OBJECT.')
               IF(NMERG.EQ.NMAT) THEN
                 CALL LCMGET(KPMACR,'DIFF',DIFHET(1,IG))
               ELSE
                 CALL LCMGET(KPMACR,'DIFF',DFL)
                 CALL LCMGET(KPMACR,'FLUX-INTG',FL)
                 DO INM=1,NMERG
                   DDELN1=0.D0
                   DDELD1=0.D0
                   DO IBM=1,NMAT
                     IF(IMERG(IBM).EQ.INM) THEN
                       DDELN1=DDELN1+FL(IBM)/DFL(IBM)
                       DDELD1=DDELD1+FL(IBM)
                     ENDIF
                   ENDDO
                   IF(DDELN1.EQ.0.D0) CALL XABORT('FLU2DR: DDELN1=0.')
                   DIFHET(INM,IG)=REAL(DDELD1/DDELN1)
                 ENDDO
               ENDIF
            ENDDO
            DEALLOCATE(DFL,FL)
            GAMMA(:NGRP)=1.0
         ELSE IF(OPTION.EQ.'RHS') THEN
            CALL LCMLEN(IPFLUX,'DIFFHET',ILONG,ITYLCM)
            IF(ILONG.EQ.NMERG*NGRP) THEN
              IF(LFORW) THEN
                 CALL LCMGET(IPFLUX,'DIFFHET',DIFHET)
              ELSE
*                Permute the diffusion coefficients if the LKRD option
*                is set for an adjoint calculation
                 CALL LCMGET(IPFLUX,'DIFFHET',DIFHET)
                 DO INM=1,NMERG
                   GAMMA(:NGRP)=DIFHET(INM,:NGRP)
                   DO IG=1,NGRP
                     DIFHET(INM,IG)=GAMMA(NGRP-IG+1)
                   ENDDO
                 ENDDO
              ENDIF
            ELSE
               CALL XABORT('FLU2DR: UNABLE TO RECOVER THE DIFFHET RECO'
     >         //'RD IN THE FLUX OBJECT.(1)')
            ENDIF
            GAMMA(:NGRP)=1.0
         ELSE IF((LEAKSW).OR.(ILEAK.EQ.5)) THEN
*           Todorova heterogeneous leakage model.
            CALL FLULPN(IPMACR,NUNKNO,OPTION,'DIFF',NGRP,NREG,NMAT,
     1      VOL,MATCOD,NMERG,IMERG,KEYFLX(1,1,1),FLUX(1,1,2),B2(4),
     2      IPRT,DIFHET,DHOM)
            GAMMA(:NGRP)=1.0
         ELSE
*           FUNDAMENTAL MODE CONDITION.
            IF(NMERG.NE.1) CALL XABORT('FLU2DR: ONE LEAKAGE ZONE EXPEC'
     1      //'TED.(1)')
            CALL B1HOM(IPMACR,LEAKSW,NUNKNO,OPTION,'DIFF',NGRP,NREG,
     1      NMAT,NIFIS,VOL,MATCOD,KEYFLX(1,1,1),FLUX(1,1,2),REFKEF,
     2      IPRT,DIFHET(1,1),GAMMA,AKEFF,INORM,B2)
         ENDIF
         CALL LCMPUT(IPFLUX,'B2  B1HOM',1,2,B2(4))
         CALL LCMPUT(IPFLUX,'DIFFHET',NMERG*NGRP,2,DIFHET)
      ENDIF
*
****  OUTER LOOP  ******************************************
      IGDEB=1
      CFLI=.FALSE.
      CEXE=.FALSE.
      DO 400 IT=1,MAXOUT
      CALL KDRCPU(CPU1)
      MESSIN='     '
*----
*  FISSION SOURCE CALCULATION IN FLUX(:,:,4)
*----
      IF(((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2))
     1   .AND.(.NOT.C_ASSOCIATED(IPSOU))) THEN
        ALLOCATE(FXSOR(0:NMAT))
        JPMACR=LCMGID(IPMACR,'GROUP')
        FLUX(:NUNKNO,:NGRP,4)=0.0
        DO IG=1,NGRP
          KPMACR=LCMGIL(JPMACR,IG)
          FXSOR(0)=0.0
          CALL LCMGET(KPMACR,'FIXE',FXSOR(1))
          CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,FXSOR,
     >    FLUX(1,IG,4))
        ENDDO
        DEALLOCATE(FXSOR)
      ELSE IF((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)) THEN
        DO IG=1,NGRP
          IF(LFORW) THEN
            CALL LCMGDL(JPSOU,IG,FLUX(1,IG,4))
          ELSE
            CALL LCMGDL(JPSOU,NGRP-IG+1,FLUX(1,IG,4))
          ENDIF
        ENDDO
      ELSE IF(ITYPEC.EQ.1) THEN
        DO IG=1,NGRP
          IF(LFORW) THEN
            CALL LCMGDL(JPSOU,IG,FLUX(1,IG,4))
          ELSE
            CALL LCMGDL(JPSOU,NGRP-IG+1,FLUX(1,IG,4))
          ENDIF
          DO IUN=1,NUNKNO
            FLUX(IUN,IG,4)=-FLUX(IUN,IG,4)
          ENDDO
        ENDDO
      ELSE
        FLUX(:NUNKNO,:NGRP,4)=0.0
      ENDIF
      IF(NIFIS.GT.0) THEN
        IF(CXDOOR.EQ.'BIVAC') THEN
          CALL BIVFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL,
     >    XSCHI,XSNUF,FLUX(1,1,2),FLUX(1,1,4))
        ELSE IF(CXDOOR.EQ.'TRIVAC') THEN
          CALL TRIFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL,
     >    XSCHI,XSNUF,FLUX(1,1,2),FLUX(1,1,4))
        ELSE
          ALLOCATE(FXSOR(NUNKNO))
          DO IS=1,NIFIS
            FXSOR(:NUNKNO)=0.0
            DO IG=1,NGRP
              CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,XSNUF(0,IS,IG),
     >        FXSOR,FLUX(1,IG,2))
            ENDDO
            DO IR=1,NREG
              IBM=MATCOD(IR)
              IF(IBM.EQ.0) CYCLE
              DO IE=1,NLIN
                IND=KEYFLX(IR,IE,1)
                IF(IND.EQ.0) CYCLE
                DO IG=1,NGRP
                  FLUX(IND,IG,4)=FLUX(IND,IG,4)+XSCHI(IBM,IS,IG)*
     >            FXSOR(IND)
                ENDDO
              ENDDO
            ENDDO
          ENDDO ! IS
          DEALLOCATE(FXSOR)
        ENDIF
        FLUX(:NUNKNO,:NGRP,4)=FLUX(:NUNKNO,:NGRP,4)*REAL(AFLNOR)
      ENDIF
*----
*  VOLUME-INTEGRATED SOURCE CALCULATION FOR USE IN NEUTRON BALANCE
*----
      DO IG=1,NGRP
        XCSOU(IG)=0.0D0
        DO IR=1,NREG
          IND=KEYFLX(IR,1,1)
          IF((CXDOOR.EQ.'BIVAC').OR.(CXDOOR.EQ.'TRIVAC')) THEN
            ! volumes are already included in the sources
            IF(IND.GT.0) XCSOU(IG)=XCSOU(IG)+FLUX(IND,IG,4)
          ELSE
            IF(IND.GT.0) XCSOU(IG)=XCSOU(IG)+FLUX(IND,IG,4)*VOL(IR)
          ENDIF
        ENDDO
      ENDDO
      ISBS=-1
      IF(C_ASSOCIATED(IPSOU)) CALL LCMLEN(IPSOU,'NBS',ISBS,ITYLCM)
*----
*  SET THE STARTING ENERGY GROUP
*----
      DO 40 IG=1,NGRP
      IGDEB=IG
      IF(XCSOU(IG).NE.0.0.OR.ISBS.NE.0) GO TO 50
   40 CONTINUE
*----
*  DOWNLOAD FROM EXTERNAL FLUX(:,:,2) TO PRESENT INTERNAL FLUX(:,:,6)
*----
   50 FLUX(:NUNKNO,:NGRP,6)=FLUX(:NUNKNO,:NGRP,2)
*
****  INNER LOOP  ******************************************
      DO 270 JT=1,MAXINR
      FLUX(:NUNKNO,IGDEB:NGRP,7)=FLUX(:NUNKNO,IGDEB:NGRP,6)
      FLUX(:NUNKNO,IGDEB:NGRP,8)=FLUX(:NUNKNO,IGDEB:NGRP,4)
      JPMACR=LCMGID(IPMACR,'GROUP')
      DO 140 IG=IGDEB,NGRP
*----
*  PROCESS SELF-SCATTERING REDUCTION IN INNER SOURCES.
*----
      IF((ITPIJ.EQ.2).OR.(ITPIJ.EQ.4)) THEN
         IF((CXDOOR.EQ.'BIVAC').OR.(CXDOOR.EQ.'TRIVAC')) THEN
            CALL XABORT('FLU2DR: SCATTERING REDUCTION IS MANDATORY.')
         ENDIF
         CALL DOORS(CXDOOR,IPTRK,NMAT,NANIS,NUNKNO,XSDIA(0,0,IG),
     >   FLUX(1,IG,8),FLUX(1,IG,7))
      ENDIF
      IF(ILEAK.EQ.6) THEN
*        ECCO ISOTROPIC STREAMING MODEL.
         CCLBD=0.0
         IF((ITPIJ.EQ.1).OR.(ITPIJ.EQ.3).AND.(OPTION.EQ.'B1')) 
     >    CCLBD=1.0-GAMMA(IG)    
         DO 75 IE=1,NLIN
         DO 70 IR=1,NREG
         IBM=MATCOD(IR)
         IND=NUNKNO/2+KEYFLX(IR,IE,1)
         IF(IND.EQ.NUNKNO/2) GO TO 70
         IF(OPTION(2:2).EQ.'1') THEN
*           B1 OR P1 CASE.
            IF(ITPIJ.EQ.2) THEN
               FLUX(IND,IG,8)=FLUX(IND,IG,8)+XSDIA(IBM,1,IG)*
     >         FLUX(IND,IG,7)
            ENDIF
         ELSE IF(ITPIJ.EQ.1) THEN
*           B0, P0, B0TR OR P0TR CASE.
            FLUX(IND,IG,8)=FLUX(IND,IG,8)-XSDIA(IBM,1,IG)*
     >      FLUX(IND,IG,7)*GAMMA(IG)
         ENDIF
         FLUX(IND,IG,8)=FLUX(IND,IG,8)+CCLBD*XSDIA(IBM,1,IG)*
     >   FLUX(IND,IG,7)
  70     CONTINUE
  75     CONTINUE
      ELSE IF(ILEAK.GE.7) THEN
*        TIBERE ANISOTROPIC STREAMING MODEL.
         CCLBD=0.0
         IF((ITPIJ.EQ.3).AND.(OPTION.EQ.'B1')) CCLBD=1.0-GAMMA(IG)
         DO 86 IE=1,NLIN
         DO 85 IR=1,NREG
         IND0=KEYFLX(IR,IE,1)
         IF(IND0.EQ.0) GO TO 85
         IBM=MATCOD(IR)
         INDD(1)=NUNKNO/4+IND0
         INDD(2)=NUNKNO/2+IND0
         INDD(3)=3*NUNKNO/4+IND0
         DO 80 IDIR=1,3
         IND=INDD(IDIR)
         IF(OPTION(2:2).EQ.'1') THEN
*           B1 OR P1 CASE.
            IF(ITPIJ.EQ.4) THEN
               FLUX(IND,IG,8)=FLUX(IND,IG,8)+XSDIA(IBM,1,IG)*
     >         FLUX(IND,IG,7)
            ENDIF
         ELSE IF(ITPIJ.EQ.3) THEN
*           B0, P0, B0TR OR P0TR CASE.
            FLUX(IND,IG,8)=FLUX(IND,IG,8)-XSDIA(IBM,1,IG)*
     >      FLUX(IND,IG,7)*GAMMA(IG)
         ENDIF
         FLUX(IND,IG,8)=FLUX(IND,IG,8)+CCLBD*XSDIA(IBM,1,IG)*
     >   FLUX(IND,IG,7)
  80     CONTINUE
  85     CONTINUE
  86     CONTINUE
      ENDIF
*----
*  COMPUTE INNER SOURCES ASSUMING SELF-SCATTERING REDUCTION.
*----
      IF(.NOT.LSCAL) THEN 
         KPMACR=LCMGIL(JPMACR,IG)
         NUNK2=NUNKNO
         IF(ILEAK.EQ.6) NUNK2=NUNKNO/2
         IF(ILEAK.GE.7) NUNK2=NUNKNO/4
         IF((CXDOOR.EQ.'SN').AND.(IBFP.EQ.0)) THEN
            NUNK2=NUNKNO
            CALL SNSOUR(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2,
     1      NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8))
         ELSE IF(CXDOOR.EQ.'SN') THEN
            NUNK2=NUNKNO
            IPSTR=LCMGID(IPSYS,'STREAMING')
            JPSTR=LCMGID(IPSTR,'GROUP')
            KPSYS=LCMGIL(JPSTR,IG)
            CALL SNSBFP(IG,IPTRK,KPMACR,KPSYS,NANIS,NLF,NREG,NMAT,
     1      NUNK2,NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8))
         ELSE
            HLEAK='      '
            CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS,
     1      NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8))
         ENDIF
         IF((ILEAK.EQ.6).AND.(OPTION(2:2).EQ.'1')) THEN
*           ECCO ISOTROPIC STREAMING MODEL.
            HLEAK='ECCO  '
            CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS,
     1      NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8))
         ELSE IF(ILEAK.GE.7) THEN
*           TIBERE ANISOTROPIC STREAMING MODEL.
            HLEAK='TIBERE'
            CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS,
     1      NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8))
         ENDIF
      ENDIF
  140 CONTINUE
*----
*  FLUX COMPUTATION
*----
      NPSYS(:NGRP)=0
      DO 150 IG=IGDEB,NGRP
      NPSYS(IG)=IG
  150 CONTINUE
      JPSTR=C_NULL_PTR
      IF(C_ASSOCIATED(IPSYS)) THEN
         JPSYS=LCMGID(IPSYS,'GROUP')
         IF(ILEAK.EQ.6.OR.((MOD(ILEAK,10).EQ.7).AND.(IPHASE.EQ.1))) THEN
            IPSTR=LCMGID(IPSYS,'STREAMING')
            JPSTR=LCMGID(IPSTR,'GROUP')
         ENDIF
      ENDIF
      IF((.NOT.LSCAL).AND.(ILEAK.EQ.0)) THEN
         IDIR=0
         CALL DOORFV(CXDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRT,NGRP,
     1   NMAT,IDIR,NREG,NUNKNO,IPHASE,LEXAC,MATCOD,VOL,KEYFLX,TITLE,
     2   FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,FLUXC,EVALRHO)
      ELSE IF(.NOT.LSCAL) THEN
         CALL FLUDBV(CXDOOR,IPHASE,JPSYS,JPSTR,NPSYS,IPTRK,IFTRAK,
     1   IPRT,NREG,NUNKNO,NFUNL,NGRP,NMAT,NANIS,LEXAC,MATCOD,VOL,KEYFLX,
     2   TITLE,ILEAK,LEAKSW,XSTRC,XSDIA,B2,NMERG,IMERG,DIFHET,GAMMA,
     3   FLUX(1,1,2),FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,FLUXC)
      ELSE
*        A GROUP-BY-GROUP SCALAR PROCEDURE IS BEEN USED.
         IF(.NOT.C_ASSOCIATED(IPSYS)) THEN
            CALL XABORT('FLU2DR: MISSING L_PIJ OBJECT.')
         ENDIF
         KPSTR=C_NULL_PTR
         DO 230 IG=IGDEB,NGRP
         IF(IPRT.GT.10) WRITE(6,'(/25H FLU2DR: PROCESSING GROUP,I5,
     >   6H WITH ,A,1H.)') IG,CXDOOR
         KPMACR=LCMGIL(JPMACR,IG)
         NUNK2=NUNKNO
         IF(ILEAK.EQ.6) NUNK2=NUNKNO/2
         IF(ILEAK.GE.7) NUNK2=NUNKNO/4
         IF(CXDOOR.EQ.'BIVAC') THEN
            CALL BIVSOU(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2,
     1      NGRP,MATCOD,VOL,FLUX(1,1,7),FLUX(1,1,8))
         ELSE IF(CXDOOR.EQ.'TRIVAC') THEN
            CALL TRIVSO(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2,
     1      NGRP,MATCOD,VOL,FLUX(1,1,7),FLUX(1,1,8))
         ELSE IF((CXDOOR.EQ.'SN').AND.(IBFP.EQ.0)) THEN
            CALL SNSOUR(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2,
     1      NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8))
         ELSE IF(CXDOOR.EQ.'SN') THEN
            JPSYS=LCMGID(IPSYS,'GROUP')
            KPSYS=LCMGIL(JPSYS,IG)
            CALL SNSBFP(IG,IPTRK,KPMACR,KPSYS,NANIS,NLF,NREG,NMAT,
     1      NUNK2,NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8))
         ELSE
            HLEAK='      '
            CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS,
     1      NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8))
         ENDIF
         IF((ILEAK.EQ.6).AND.(OPTION(2:2).EQ.'1')) THEN
*           ECCO ISOTROPIC STREAMING MODEL.
            HLEAK='ECCO  '
            CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS,
     1      NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8))
         ELSE IF(ILEAK.GE.7) THEN
*           TIBERE ANISOTROPIC STREAMING MODEL.
            HLEAK='TIBERE'
            CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS,
     1      NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8))
         ENDIF
*
         NPSYS(:NGRP)=0
         NPSYS(IG)=IG
         IF(ILEAK.EQ.0) THEN
           IDIR=0
           CALL DOORFV(CXDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRT,NGRP,
     1     NMAT,IDIR,NREG,NUNKNO,IPHASE,LEXAC,MATCOD,VOL,KEYFLX,TITLE,
     2     FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,FLUXC,EVALRHO)
         ELSE
           CALL FLUDBV(CXDOOR,IPHASE,JPSYS,JPSTR,NPSYS,IPTRK,IFTRAK,
     1     IPRT,NREG,NUNKNO,NFUNL,NGRP,NMAT,NANIS,LEXAC,MATCOD,VOL,
     2     KEYFLX,TITLE,ILEAK,LEAKSW,XSTRC,XSDIA,B2,NMERG,IMERG,DIFHET,
     3     GAMMA,FLUX(1,1,2),FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,
     4     FLUXC)
         ENDIF
  230    CONTINUE
      ENDIF
      IF(LREBAL.AND.(ITYPEC.NE.5)) THEN
         CALL FLUBAL(IPMACR,NGRP,ILEAK,NMAT,NREG,ICREB,NUNKNO,NANIS,
     1   MATCOD,VOL,KEYFLX(1,1,1),XSTRC,XSDIA,XCSOU,IGDEB,B2,NMERG,
     2   IMERG,DIFHET,KEYCUR,MATALB(NNN+1),ALBEDO,V(NNN+1),FLUX(1,1,7))
      ENDIF
*----
*  ACCELERATING INNER ITERATIONS CYCLICALLY DEPENDING ON PARAM.
*----
      IF(MOD(JT-1,NCTOT).GE.NCPTM) THEN
         CALL FLU2AC(NGRP,NUNKNO,IGDEB,FLUX(1,1,5),AKEEP(5),ZMU)
      ELSE
         ZMU=1.0
      ENDIF
*----
*  CALCULATING ERROR AND PREC BETWEEN PRESENT AND NEW FLUX FOR
*  EACH GROUP. RETAIN LARGEST ERROR BETWEEN ANY GROUP.
*----
      EINN=0.0
      ICHAR=ICHAR+1
      ICTOT=ICTOT+NGRP-IGDEB+1
      IGDEBO=IGDEB
      DO 260 IG=IGDEBO,NGRP
      GINN=0.0
      FINN=0.0
      DO 240 IR=1,NREG
      IND=KEYFLX(IR,1,1)
      IF(IND.EQ.0) GO TO 240
      GINN=MAX(GINN,ABS(FLUX(IND,IG,6)-FLUX(IND,IG,7)))
      FINN=MAX(FINN,ABS(FLUX(IND,IG,7)))
  240 CONTINUE
      FLUX(:NUNKNO,IG,5)=FLUX(:NUNKNO,IG,6)
      FLUX(:NUNKNO,IG,6)=FLUX(:NUNKNO,IG,7)
      GINN=GINN/FINN
      IF((GINN.LT.EPSINR).AND.(IGDEB.EQ.IG)) THEN
         IGDEB=IGDEB+1
      ELSEIF((IGDEB.EQ.IG).AND.(IG.LT.NGRP)) THEN
         ERRDEB1=GINN
      ENDIF
      EINN=MAX(EINN,GINN)
  260 CONTINUE
*
      ITERF=JT
      IF(IPRT.GT.0) WRITE(6,1080) JT,EINN,EPSINR,IGDEB,ZMU
      IF((IPRT.GT.0).AND.(IGDEB.GT.1).AND.(IGDEB.LE.NGRP)) THEN
         WRITE(6,1082) ERRDEB1 
      ENDIF
      IF(EINN.LT.EPSINR) THEN
*        thermal convergence is reached
         CFLI=CEXE
         GOTO 280
      ENDIF
*     near convergence (eps < 10.0 criterion) a new outer iteration
*     is started
      IF((IGDEB.GT.1).AND.(EINN.LT.10.*EPSINR)) GOTO 281
  270 CONTINUE
      MESSIN='*NOT*'
****  END OF INNER LOOP  ******************************************
*
  281 MESSIN='*NEARLY*'
  280 IF(LREBAL) THEN
         IF(LEAKSW) THEN
            IF(ICREB.EQ.0) THEN
               WRITE(6,*) ' *** INCOMPATIBILITY ON LEAKAGE SWITCH ***'
               CALL XABORT('FLU2DR: ERROR ON LEAKAGE SWITCH')
            ELSE
               IF(IPRT.GT.0) 
     &            WRITE(6,*) 'FLU2DR: LEAKAGE & ICREB -> REBALANCING ON'
            ENDIF
         ELSE
            IF(IPRT.GT.0) 
     &         WRITE(6,*) 'FLU2DR: NO LEAKAGE-> REBALANCING ON'
         ENDIF
      ELSE
         IF(IPRT.GT.0) WRITE(6,*) 'FLU2DR:    LEAKAGE-> REBALANCING OFF'
      ENDIF
      CALL KDRCPU(CPU2)
      IF(IPRT.GT.0) WRITE(6,1040) CPU2-CPU1,'INTERNAL',MESSIN,ITERF
*----
*  PROMOTE FROM NEW INTERNAL FLUX(,,,7) TO NEW EXTERNAL FLUX(,,,3)
*----
      FLUX(:NUNKNO,:NGRP,3)=FLUX(:NUNKNO,:NGRP,7)
      FLUX(:NUNKNO,:NGRP,4)=FLUX(:NUNKNO,:NGRP,8)
*----
*  HOTELLING DEFLATION IN GPT CASES.
*----
      IF(ITYPEC.EQ.1) THEN
        JPFLUP1=LCMGID(IPFLUP,'FLUX')
        JPFLUP2=LCMGID(IPFLUP,'AFLUX')
        DDELN1=0.0D0
        DDELD1=0.0D0
        DO 300 IG=1,NGRP
        IF(LFORW) THEN
          CALL LCMGDL(JPFLUP1,IG,FLUX(1,IG,5)) ! EVECT
          CALL LCMGDL(JPFLUP2,IG,FLUX(1,IG,6)) ! ADECT
        ELSE
          CALL LCMGDL(JPFLUP2,NGRP-IG+1,FLUX(1,IG,5)) ! ADECT
          CALL LCMGDL(JPFLUP1,NGRP-IG+1,FLUX(1,IG,6)) ! EVECT
        ENDIF
  300   CONTINUE
        FLUX(:NUNKNO,:NGRP,7)=0.0
        IF(CXDOOR.EQ.'BIVAC') THEN
          CALL BIVFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL,
     >    XSCHI,XSNUF,FLUX(1,1,6),FLUX(1,1,7))
        ELSE IF(CXDOOR.EQ.'TRIVAC') THEN
          CALL TRIFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL,
     >    XSCHI,XSNUF,FLUX(1,1,6),FLUX(1,1,7))
        ELSE
          ALLOCATE(FXSOR(NUNKNO))
          DO IS=1,NIFIS
            FXSOR(:NUNKNO)=0.0
            DO IG=1,NGRP
              CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,XSNUF(0,IS,IG),
     >        FXSOR,FLUX(1,IG,6))
            ENDDO
            DO IR=1,NREG
              IBM=MATCOD(IR)
              IF(IBM.EQ.0) CYCLE
              DO IE=1,NLIN
                IND=KEYFLX(IR,IE,1)
                IF(IND.EQ.0) CYCLE
                DO IG=1,NGRP
                  FLUX(IND,IG,7)=FLUX(IND,IG,7)+XSCHI(IBM,IS,IG)*
     >            FXSOR(IND)
                ENDDO
              ENDDO
            ENDDO
          ENDDO ! IS
          DEALLOCATE(FXSOR)
        ENDIF
*
        DO 335 IG=1,NGRP
        DO 330 IND=1,NUNKNO
        DDELN1=DDELN1+FLUX(IND,IG,7)*FLUX(IND,IG,3)
        DDELD1=DDELD1+FLUX(IND,IG,7)*FLUX(IND,IG,5)
  330   CONTINUE
  335   CONTINUE
        DO 345 IG=1,NGRP
        DO 340 IND=1,NUNKNO
        FLUX(IND,IG,3)=FLUX(IND,IG,3)-REAL(DDELN1/DDELD1)*FLUX(IND,IG,5)
  340   CONTINUE
  345   CONTINUE
      ENDIF
*
      IF(ITYPEC.EQ.2) THEN
*         NO B-N LEAKAGE CALCULATION REQUIRED
          IF(B2(4).NE.0.0) CALL XABORT('FLU2DR: NON ZERO BUCKLING.')
          CALL FLUKEF(IPRT,IPMACR,NGRP,NREG,NUNKNO,NMAT,NIFIS,NANIS,
     1    MATCOD(1),VOL,KEYFLX(1,1,1),XSTRC,XSDIA,XSNUF,XSCHI,NMERG,
     2    IMERG,DIFHET,FLUX(1,1,3),B2,ILEAK,LEAKSW,OLDBIL,AKEFF,AFLNOR)
      ELSE IF(ITYPEC.GT.2) THEN
*        PERFORM LEAKAGE CALCULATION.
         CALL LCMLEN(IPFLUX,'DIFFHET',ILONG,ITYLCM)
         IF(ILONG.EQ.0) THEN
            CALL XABORT('FLU2DR: UNABLE TO RECOVER THE DIFFHET RECORD '
     1      //'IN THE FLUX OBJECT.(2)')
         ENDIF
         CALL LCMGET(IPFLUX,'DIFFHET',DIFHET)
         GAMMA(:NGRP)=1.0
         IF(ILEAK.EQ.5) THEN
*           Todorova heterogeneous leakage model.
            CALL FLULPN(IPMACR,NUNKNO,OPTION,HTYPE(ITYPEC),NGRP,NREG,
     1      NMAT,VOL,MATCOD,NMERG,IMERG,KEYFLX(1,1,1),FLUX(1,1,3),B2(4),
     2      IPRT,DIFHET,DHOM)
            IF(.NOT.LEAKSW) THEN
              CALL B1HOM(IPMACR,LEAKSW,NUNKNO,'LKRD',HTYPE(ITYPEC),NGRP,
     1        NREG,NMAT,NIFIS,VOL,MATCOD,KEYFLX(1,1,1),FLUX(1,1,3),
     2        REFKEF,IPRT,DHOM(1),GAMMA,AKEFF,INORM,B2)
              GO TO 350
            ENDIF
         ENDIF
         IF(LEAKSW) THEN
            IF(HTYPE(ITYPEC).NE.'K') THEN
              CALL XABORT('FLU2DR: TYPE K EXPECTED WITH OPEN GEOMETRY.')
            ENDIF
            JPMACR=LCMGID(IPMACR,'GROUP')
            ALLOCATE(SFNU(NMAT,NIFIS))
            PROD=0.0D0
            DO IGR=1,NGRP
              KPMACR=LCMGIL(JPMACR,IGR)
              SFNU(:NMAT,:NIFIS)=0.0
              IF(NIFIS.GT.0) CALL LCMGET(KPMACR,'NUSIGF',SFNU)
              DO IBM=1,NMAT
                FLXIN=0.0D0
                DO I=1,NREG
                  IND=KEYFLX(I,1,1)
                  IF((MATCOD(I).EQ.IBM).AND.(IND.GT.0)) THEN
                    FLXIN=FLXIN+FLUX(IND,IGR,3)*VOL(I)
                  ENDIF
                ENDDO
                DO NF=1,NIFIS
                  PROD=PROD+SFNU(IBM,NF)*FLXIN
                ENDDO
              ENDDO
            ENDDO
            DEALLOCATE(SFNU)
            AKEFF=AKEFF*PROD/OLDBIL
            OLDBIL=PROD
            IF(IPRT.GT.0) WRITE (6,1150) B2(4),AKEFF
         ELSE
*           FUNDAMENTAL MODE CONDITION.
            IF(NMERG.NE.1) CALL XABORT('FLU2DR: ONE LEAKAGE ZONE EXPEC'
     1      //'TED.(2)')
            CALL B1HOM(IPMACR,LEAKSW,NUNKNO,OPTION,HTYPE(ITYPEC),NGRP,
     1      NREG,NMAT,NIFIS,VOL,MATCOD,KEYFLX(1,1,1),FLUX(1,1,3),
     2      REFKEF,IPRT,DIFHET(1,1),GAMMA,AKEFF,INORM,B2)
            IF(ILEAK.GE.7) THEN
*              COMPUTE THE DIRECTIONNAL BUCKLING COMPONENTS FOR TIBERE.
               IHETL=ILEAK/10-1
               IF(IHETL.GT.0) THEN
                  CALL FLUBLN(IPMACR,IPRT,NGRP,NMAT,NREG,NUNKNO,NIFIS,
     1            MATCOD,VOL,KEYFLX(1,1,1),FLUX(1,1,3),IHETL,REFKEF,B2)
               ENDIF
             ENDIF
          ENDIF
  350     CALL LCMPUT(IPFLUX,'B2  B1HOM',1,2,B2(4))
          CALL LCMPUT(IPFLUX,'DIFFHET',NMERG*NGRP,2,DIFHET)
      ENDIF
      IF(ITYPEC.GE.3) THEN
         IF(B2(4).EQ.0.0) THEN
           BFLNOR=1.0D0
         ELSE
           BFLNOR=1.0D0/ABS(B2(4))
         ENDIF
         EEXT=REAL(ABS(B2(4)-B2VALO)*BFLNOR)
         B2VALO=B2(4)
      ENDIF
      IEXTF=IT
      IF((ITYPEC.GT.1).AND.(ITYPEC.LT.5)) THEN
         IF(AKEFF.NE.0.0) AFLNOR=1.0D0/AKEFF
         EEXT=REAL(ABS(AKEFF-AKEFFO)/AKEFF)
      ELSE
         EEXT=0.0
      ENDIF
      AKEEP(3)=AKEFF
*----
*  ACCELERATING INNER ITERATIONS CYCLICALLY DEPENDING ON PARAM.
*----
      IF(MOD(IT-1,NCTOT).GE.NCPTM) THEN
         CALL FLU2AC(NGRP,NUNKNO,1,FLUX(1,1,1),AKEEP(1),ZMU)
      ELSE
         ZMU=1.0
      ENDIF
*
      EINN=0.0
      IF(IPRT.GT.0) WRITE(6,1090) IT,EEXT,EPSOUT,AKEFF,B2(4)
      IF(EEXT.LT.EPSOUT) THEN
*        COMPARE FLUX FOR OUTER ITERATIONS
         DO 370 IG=1,NGRP
         GINN=0.0
         FINN=0.0
         DO 360 IR=1,NREG
         IND=KEYFLX(IR,1,1)
         IF(IND.EQ.0) GO TO 360
         GINN=MAX(GINN,ABS(FLUX(IND,IG,2)-FLUX(IND,IG,3)))
         FINN=MAX(FINN,ABS(FLUX(IND,IG,3)))
  360    CONTINUE
         FLUX(:NUNKNO,IG,1)=FLUX(:NUNKNO,IG,2)
         FLUX(:NUNKNO,IG,2)=FLUX(:NUNKNO,IG,3)
         GINN=GINN/FINN
         EINN=MAX(EINN,GINN)
  370    CONTINUE
         IF(IPRT.GT.0) WRITE(6,1100) IT,EINN,EPSUNK,AFLNOR,ZMU
         CEXE=.TRUE.
      ELSE
         FLUX(:NUNKNO,:NGRP,1)=FLUX(:NUNKNO,:NGRP,2)
         FLUX(:NUNKNO,:NGRP,2)=FLUX(:NUNKNO,:NGRP,3)
         IF(IPRT.GT.0) WRITE(6,1110) IT,AFLNOR,ZMU
      ENDIF
      IF((ITYPEC.GE.2).AND.(AKEFF.NE.0.0)) THEN
         AFLNOR=(AKEFF/AKEEP(3))*AFLNOR
      ENDIF
      AKEEP(1)=AKEEP(2)
      AKEEP(2)=AKEEP(3)
*----
*  UPDATE KEFF
*----
      AKEFFO=AKEFF
      IF((EEXT.LT.EPSOUT).AND.(EINN.LT.EPSUNK).AND.(IT.GE.2)) GO TO 410
  400 CONTINUE
      WRITE(6,*) '*** FLU2DR: CONVERGENCE NOT REACHED ***'
      WRITE(6,*) '*** FLU2DR: CONVERGENCE NOT REACHED ***'
      WRITE(6,*) '*** FLU2DR: CONVERGENCE NOT REACHED ***'
      MESSOU='*NOT*'
*
****  CONVERGENCE REACHED  ******************************************
  410 RKEFF=REAL(AKEFF)
      IF(IPRT.GE.3) THEN
         WRITE(6,1010) (IR,IR=1,NREG)
         ALLOCATE(FL(NREG))
         DO 425 IG=1,NGRP
         WRITE(6,1070) IG
         FL(:NREG)=0.0
         DO 420 IR=1,NREG
         IND=KEYFLX(IR,1,1)
         IF(IND.GT.0) FL(IR)=FLUX(IND,IG,3)
  420    CONTINUE
         WRITE(6,1020) (FL(IR),IR=1,NREG)
  425    CONTINUE
         DEALLOCATE(FL)
      ENDIF
      IF(IPRT.GE.4) THEN
         ALLOCATE(FL(NREG))
         DO 445 IG=1,NGRP
         WRITE(6,1070) IG
         DO 440 IA=2,NFUNL
         FL(:NREG)=0.0
         DO 430 IR=1,NREG
         IND=KEYFLX(IR,1,IA)
         IF(IND.GT.0) FL(IR)=FLUX(IND,IG,3)
  430    CONTINUE
         WRITE(6,1030) IA,(FL(IR),IR=1,NREG)
  440    CONTINUE
  445    CONTINUE
         DEALLOCATE(FL)
      ENDIF
*----
*  COMPUTE K-INF
*----
      IF(ITYPEC.GE.2) THEN
        FISOUR=0.0D0
        OLDBIL=0.0D0
        DO 490 IG=1,NGRP
        DO 460 IR=1,NREG
        IND=KEYFLX(IR,1,1)
        IF(IND.EQ.0) GO TO 460
        DO 450 IS=1,NIFIS
        FISOUR=FISOUR+XSNUF(MATCOD(IR),IS,IG)*FLUX(IND,IG,3)*VOL(IR)
  450   CONTINUE
        OLDBIL=OLDBIL+XSTRC(MATCOD(IR),IG)*FLUX(IND,IG,3)*VOL(IR)
  460   CONTINUE
        KPMACR=LCMGIL(JPMACR,IG)
        CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
        CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
        CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
        CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
        DO 480 IR=1,NREG
        IBM=MATCOD(IR)
        IF(IBM.GT.0) THEN
          IND=KEYFLX(IR,1,1)
          JG=IJJ(IBM)
          DO 470 JND=1,NJJ(IBM)
          IF(JG.EQ.IG) THEN
            OLDBIL=OLDBIL-XSDIA(IBM,0,IG)*FLUX(IND,JG,3)*VOL(IR)
          ELSE
            OLDBIL=OLDBIL-XSCAT(IPOS(IBM)+JND-1)*FLUX(IND,JG,3)*VOL(IR)
          ENDIF
          JG=JG-1
  470     CONTINUE
        ENDIF
  480   CONTINUE
  490   CONTINUE
        CUREIN=0.0
        IF(FISOUR.NE.0.0) CUREIN=REAL(FISOUR/OLDBIL)
*
*       FLUX NORMALIZATION TO KEFF.
        IF(ITYPEC.LT.5) THEN
          FLUX(:NUNKNO,:NGRP,3)=FLUX(:NUNKNO,:NGRP,3)*REAL(AKEFF/FISOUR)
          FLUX(:NUNKNO,:NGRP,4)=FLUX(:NUNKNO,:NGRP,4)*REAL(AKEFF/FISOUR)
        ENDIF
      ENDIF
*----
*  PRINT TIME TAKEN
*----
      CALL KDRCPU(CPU1)
      IF(IPRT.GE.1) WRITE(6,1040) CPU1-CPU0,'EXTERNAL',MESSOU,IEXTF
*----
*  FOURIER ANALYSIS: STORE E-VALUE AND FIND SPECTRAL RADIUS
*----
      IF(ITYPEC.EQ.-2) THEN
         WRITE(6,1130)
         ARRAYRHO1(IFACOUNT+1)=EVALRHO
         IF(IFACOUNT.LT.(NFOU-1)) GO TO 26
         SPECR=MAXVAL(ARRAYRHO1)
         WRITE(6,1140) SPECR
         CALL LCMPUT(IPFLUX,'SPEC-RADIUS',1,2,SPECR)
      ENDIF
      DEALLOCATE(ARRAYRHO1)
*----
*  PRINT TRACKING INFORMATION
*----
      IF(IPRT.GE.1) THEN
         IF((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)) THEN
            WRITE(6,1050) ICHAR,EEXT
         ELSE
            WRITE(6,1060) ICHAR,CUREIN,AKEFF,B2(4),EEXT
         ENDIF
         WRITE(6,1120) ICTOT
      ENDIF
*----
*  RELEASE ARRAYS
*----
      IF(CXDOOR.EQ.'MCCG') THEN
         IF(ICREB.GT.0) DEALLOCATE(MATALB,V,KEYCUR)
      ENDIF
*----
*  SAVE THE SOLUTION
*----
      DO 510 IG=1,NGRP
      IF(LFORW) THEN
         CALL LCMPDL(JPFLUX,IG,NUNKNO,2,FLUX(1,IG,3))
         CALL LCMPDL(JPSOUR,IG,NUNKNO,2,FLUX(1,IG,4))
      ELSE
         CALL LCMPDL(JPFLUX,NGRP-IG+1,NUNKNO,2,FLUX(1,IG,3))
         CALL LCMPDL(JPSOUR,NGRP-IG+1,NUNKNO,2,FLUX(1,IG,4))
      ENDIF
  510 CONTINUE
      IF(C_ASSOCIATED(IPSOU)) THEN
        CALL LCMLEN(IPSOU,'NORM-FS',ILEN,ITYLCM)
        IF(ILEN.GT.0) THEN
          CALL LCMGET(IPSOU,'NORM-FS',ZNORM)
          CALL LCMPUT(IPFLUX,'NORM-FS',1,2,ZNORM)
          CALL LCMPUT(IPFLUX,'MATCOD',NREG,1,MATCOD)
        ENDIF
      ENDIF
      IF(IBFP.NE.0) THEN
        CALL LCMGET(IPSYS,'ECUTOFF',ECUTOFF)
        CALL LCMPUT(IPFLUX,'ECUTOFF',1,2,ECUTOFF)
        CALL LCMPUT(IPFLUX,'FLUXC',NREG,2,FLUXC)
      ENDIF
      IF(ITYPEC.GE.2) THEN
         CALL LCMPUT(IPFLUX,'K-EFFECTIVE',1,2,RKEFF)
         CALL LCMPUT(IPFLUX,'K-INFINITY',1,2,CUREIN)
      ENDIF
      IF(ITYPEC.GE.3) THEN
         CALL LCMPUT(IPFLUX,'B2  B1HOM',1,2,B2(4))
      ENDIF
      IF((ITYPEC.GT.2).AND.(ILEAK.GE.7)) THEN
         CALL LCMPUT(IPFLUX,'B2  HETE',3,2,B2)
      ENDIF
      IF((ITYPEC.GT.2).AND.(ILEAK.GE.6)) THEN
         CALL LCMPUT(IPFLUX,'GAMMA',NGRP,2,GAMMA)
      ENDIF
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(XCSOU,DIFHET,DHOM,GAMMA,XSCAT,FLUX)
      DEALLOCATE(NPSYS,IPOS,NJJ,IJJ)
      RETURN
*
 1010 FORMAT (/28H FLUXES AVERAGED           :/
     1 (9X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,
     2 I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,
     3 2H/=,I4))
 1020 FORMAT (7H FLUX  ,2H: ,1P,10E12.5/(9X,10E12.5))
 1030 FORMAT (5H CUR ,I2,2H: ,1P,10E12.5/(9X,10E12.5))
 1040 FORMAT (18H FLU2DR: CPU TIME=, F10.0,1X,A8,13H CONVERGENCE ,
     1        A8,14H REACHED AFTER ,I6,12H ITERATIONS.  )
 1050 FORMAT (/20H ++ TRACKING CALLED=,I4,6H TIMES ,
     1         11H PRECISION=,E9.2)
 1060 FORMAT (/20H ++ TRACKING CALLED=,I4,6H TIMES ,
     1         12H FINAL KINF=,1P,E13.6,
     2         12H FINAL KEFF=,E13.6,4H B2=,E12.5,
     3         11H PRECISION=,E9.2)
 1070 FORMAT (/14H ENERGY GROUP ,I6)
 1080 FORMAT (10X,3HIN(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2,
     1 7H IGDEB=, I13,6H ACCE=,0P,F12.5)
 1082 FORMAT (18X,28HFIRST UNCONVERGED GROUP PRC=,1P,E9.2)
 1090 FORMAT (5H OUT(,I3,6H) EIG:,5H PRC=,1P,E9.2,5H TAR=,E9.2,
     1 6H KEFF=,E13.6,6H BUCK=,E12.5)
 1100 FORMAT (5H OUT(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2,
     1 6H FNOR=,E13.6,6H ACCE=,0P,F12.5)
 1110 FORMAT (5H OUT(,I3,6H) FLX:,28X,6H FNOR=,1P,E13.6,6H ACCE=,
     1 0P,F12.5)
 1120 FORMAT (38H ++ TOTAL NUMBER OF FLUX CALCULATIONS=,I10)
 1130 FORMAT (24H CONVERGENCE NOT SOUGHT.)
 1140 FORMAT (49H FLU2DR: SPECTRAL RADIUS FOR FOURIER ANALYSIS IS ,
     1 E13.6)
 1150 FORMAT(/18H FLU2DR: BUCKLING=,1P,E13.5,15H K-EFFECTIVE  =,E13.5)
      END