*DECK XHXTRK SUBROUTINE XHXTRK(IPTRK ,IPGEOM,GEONAM,IDISP ,IFTEMP, > IPRT ,NDIM ,ITOPT ,NV ,NS ,NANGL , > ISYMM ,DENS ,PCORN ,MXSEG ,ICODE ,TITREC) * *----------------------------------------------------------------------- * *Purpose: * Call all routines for the production of tracks for the * hexagonal geometry. * *Copyright: * Copyright (C) 1991 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): M. Ouisloumen * *Parameters: input * IPTRK pointer to the excell tracking. * IPGEOM pointer to the geometry. * GEONAM geometry name. * IFTEMP temporary tracking file. * IPRT print option. * TITREC title for execution. * *Parameters: input/output * IDISP tracking file disposition: * = -2 no traking - only analyse geometry * then abort (option halt); * = -1 modify tracking file; * = 0 old tracking file; * = 1 new tracking file. * *Parameters: output * NDIM number of physical dimensions. * ITOPT tracking option: * = 0 finite; = 1 cyclic. * NV number of physical regions. * NS number of outer surface. * NANGL number of angles. * ISYMM symmetry factor. * DENS track density. * PCORN corner proximity. * MXSEG maximum segment length. * ICODE albedo associated with face. * *----------------------------------------------------------------------- * USE GANLIB IMPLICIT NONE INTEGER IOUT,NSTATE CHARACTER*6 NAMSBR DOUBLE PRECISION PI PARAMETER (IOUT=6,NSTATE=40,PI=3.141592653589793D0, > NAMSBR='XHXTRK') *---- * ROUTINE PARAMETERS *---- TYPE(C_PTR) IPTRK,IPGEOM INTEGER IDISP,IFTEMP,IPRT,NDIM,MXSEG,ITOPT,NV,NS,NANGL, > ISYMM,ICODE(6) CHARACTER GEONAM*12,TITREC*72 REAL DENS,PCORN *---- * LOCAL VARIABLES *---- INTEGER KDROPN,KDRCLS INTEGER IFILE,IER INTEGER ISTATE(NSTATE),NCODE(6),IPARAM(NSTATE) INTEGER NEIGHB LOGICAL SECTOR REAL ALBEDO(6),DENUSR,DENUSZ,EXTKOP(NSTATE) DOUBLE PRECISION DBLINP,AUX,SINT,SINF,COST,COSF, > SINT1,SINF1,COST1,COSF1 DOUBLE PRECISION COS1,COS2,COS3,PASY,PASZ,PPASY,PPASZ, > RAYON,POIDS,DZMIN,DZMAX,DSIDE INTEGER NSOUT,ILNLCM,ITPLCM,NSMIN,NSMAX,NCELA,MVOLUM, > MMESH,LINMAX,LTRK,NANGLE,INDLEC,NGLE,NPHI,NTETA, > IQUAD,INDATA,IHEX,LEVEL,NCL,NCOUR,IPLANZ,NCPHY, > NREGIO,LXI,MCODE,IPLANI,NCELAP,IDIM1,IDIM2, > IAUX,IA,KGAUS,IG,JG,JAUX,MAXCYL,I,J,MAXSEC,MMAXS, > NALBG,NSURF,MACP,MTT,MT0,NBLINE,KANG,ICC,ICC1, > ICORN,IFCC1,IXX,IYY,IYY0,NCEL2,IVV,II,IAUXS, > IAUXX,IANG,ISS,ISUR REAL POIDSH,RCUTOF,CUTOFX,DEN,REDATA,SIDE,ZMIN,SQRT3 CHARACTER TEDATA*8,COMENT*80 *---- * ALLOCATABLE ARRAYS *---- INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYMRG,MATALB,LATV,ICEL, + IPLAN,IFACB,IFFV,ISURB,IVSYM,ISSYM,IMAT,KCORN,IVOIS,MATRT,SURL, + IV0 REAL, ALLOCATABLE, DIMENSION(:) :: VOLSUR,MESH,XGS,WGS,XGS1, + WGS1,WW,A1,A2,A3,T0,T1,RAUX DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ANGLES,DENSTY,POP *---- * DATA AND COMMONS *---- CHARACTER CSTOP*8,CTISO*8, CTSPC*8,CNOTR*8, CBLAN*8, > CGAUS*8, CEQW*8 SAVE CSTOP,CTISO,CTSPC,CNOTR,CBLAN,CGAUS,CEQW DATA CSTOP,CTISO,CTSPC,CNOTR,CBLAN,CGAUS,CEQW > / 'STOP','TISO','TSPC', 'NOTR',' ' , > 'GAUS','EQW' / PCORN=0.0 DENS=0.0 ISYMM=1 * ** 2.1) TAKE THE ORIGINAL GEOMETRY. ISTATE(:NSTATE)=0 CALL LCMLEN(IPGEOM,'STATE-VECTOR',ILNLCM,ITPLCM) IF (ILNLCM .GT. 0 .AND. ILNLCM .LE. NSTATE) THEN CALL LCMGET(IPGEOM,'STATE-VECTOR',ISTATE) ELSE CALL XABORT('XHXTRK: INVALID STATE VECTOR IN GEOMETRY') ENDIF CALL LCMGET(IPGEOM,'NCODE',NCODE) NDIM=2 IF(ISTATE(5).GT.0) NDIM=3 * NSMIN = 0 NSMAX = 0 NCELA = 0 MVOLUM= 0 MMESH = 0 LINMAX= 0 LTRK = 0 NANGLE= 0 DENUSR=10.0 DENUSZ=10.0 RCUTOF= 0.0 CUTOFX= 0.0 TEDATA= CBLAN IQUAD=0 * ** 1) READ ALL USER INPUT. * * READ TRACKING PARAMETERS LTRK= 0 : NO TRACKING * LTRK= 1 : ISOTROPIC TRACKING (TISO) * LTRK= 2 : SPECULAR TRACKING (TSPC) IF( IDISP.GT.0 )THEN CALL REDGET( INDLEC, NGLE, DEN, TEDATA, DBLINP) IF( INDLEC.NE.3 ) THEN CALL XABORT('XHXTRK: TYPE OF INTEGRATION NOT SPECIFIED ') ENDIF LTRK=1 IF( TEDATA(1:4).EQ.CSTOP(1:4) ) THEN * * INITIALISATION PAR DEFAUT * NPHI=5 NTETA=5 IQUAD=1 GOTO 55 ELSEIF( TEDATA(1:4).EQ.CTISO(1:4) ) THEN LTRK=1 IQUAD=2 NTETA=1 CALL REDGET( INDLEC, NPHI , DEN, TEDATA, DBLINP) IF( INDLEC .EQ. 3 ) THEN IF(TEDATA(1:4) .EQ. CGAUS(1:4)) THEN IQUAD=1 ELSE IF(TEDATA(1:3) .NE. CEQW(1:3)) THEN CALL XABORT('XHXTRK: INVALID INTEGRATION TYPE') ENDIF CALL REDGET( INDLEC, NPHI , DEN, TEDATA, DBLINP) ELSE IF( INDLEC.EQ.2 ) THEN CALL XABORT('XHXTRK: INVALID DATA FOLLOWING TISO ') ENDIF IF(NDIM .EQ. 3) THEN IF(INDLEC .NE. 1) > CALL XABORT('XHXTRK: NTETA MUST BE INTEGER') IF(NPHI .LE. 1) CALL XABORT('XHXTRK: NTETA < 2 ') NTETA=NPHI CALL REDGET( INDLEC, NPHI , DENUSR, TEDATA, DBLINP) ENDIF IF(INDLEC .NE. 1) > CALL XABORT('XHXTRK: NPHI MUST BE INTEGER') IF(NPHI .LE. 1) CALL XABORT('XHXTRK: NPHI < 2 ') CALL REDGET( INDLEC, NPHI , DENUSR, TEDATA, DBLINP) IF (INDLEC .NE. 2) > CALL XABORT('XHXTRK: DENSITY MUST BE REAL') IF(DENUSR .LE. 0.0) CALL XABORT('XHXTRK: DENSITY < 0.0 ') IF(NDIM .EQ. 3) THEN CALL REDGET( INDLEC, NGLE , DENUSZ, TEDATA, DBLINP) IF(INDLEC .NE. 2) > CALL XABORT('XHXTRK: DENSITY MUST BE REAL') IF(DENUSZ .LE. 0.0) CALL XABORT('XHXTRK: DENSITY < 0 ') ENDIF IF(IQUAD .EQ. 1) THEN IF(NTETA .GT. 64 .OR. NPHI .GT. 64) THEN CALL XABORT('XHXTRK: NANGLE IS GREATER THAN 64.') ENDIF ENDIF IF(NDIM .EQ. 2)THEN NANGLE=NTETA*NPHI*2 ELSE NANGLE=NTETA*NPHI*4 ENDIF ELSEIF( TEDATA(1:4).EQ.CTSPC(1:4) ) THEN LTRK = 2 CALL XABORT('XHXTRK: *TSPC* NOT AVAILABLE FOR HEXAGONE' ) ELSEIF( TEDATA(1:4).EQ.CNOTR(1:4) )THEN LTRK=0 GOTO 54 ELSEIF( TEDATA(1:4).NE.CGAUS(1:4) + .OR.TEDATA(1:4).NE.CEQW(1:4) ) THEN CALL XABORT('XHXTRK: CHARACTER EXPECTED ') ENDIF 54 CALL REDGET( INDLEC, INDATA, REDATA, TEDATA, DBLINP) IF( INDLEC.NE.3.OR.TEDATA(1:1).NE.';' )THEN CALL XABORT( 'XHXTRK: ; IS SUPPOSED TO BE HERE' ) ENDIF 55 CONTINUE ENDIF ** EVERYTHING IS NOW READ. *---- * SAVE EXCELL SPECIFIC TRACKING INFORMATION. *---- ISTATE(:NSTATE)=0 EXTKOP(:NSTATE)=0.0 ISTATE(9)=LTRK-1 ISTATE(11)=NANGLE ISTATE(12)=ISYMM CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,ISTATE) EXTKOP(2)=DENUSR CALL LCMPUT(IPTRK,'EXCELTRACKOP',NSTATE,2,EXTKOP) * ** 2) PROCEED TO THE EXACT GEOMETRY TREATMENT. * * * STUDY TRACKING PARAMETERS. ARE THEY BASICALLY POSSIBLE ? IF( IPRT.GT.0 )THEN WRITE(IOUT,6000) GEONAM ENDIF * -------------*************************************-------------------- * / / * / TRAITEMENT DE LA GEOMETRIE HEXAGONALE / * / / *--------------*************************************-------------------- * * DUPLICATE AND UNFOLD THE GEOMETRY. * CALL LCMGET(IPGEOM,'IHEX',IHEX) LEVEL= 1 IF(IHEX.EQ.9) THEN IF(NCODE(5) .EQ. 5 .OR. NCODE(6) .EQ. 5 .OR. > NCODE(5) .EQ. 10 .OR. NCODE(6) .EQ. 10) THEN NCL=ISTATE(3) NCOUR=NINT( (4.+SQRT(1.+4.*FLOAT(NCL-1)/3.) + +SQRT(1.+4.*FLOAT(NCL-2)/3.))*.25) ALLOCATE(LATV(2*NCOUR)) ENDIF ENDIF CALL LHXUNH(IPTRK,IPGEOM,GEONAM,MMESH,NCELA,IPLANZ,NCPHY,ICODE, + ALBEDO,NV,NREGIO,NS,SIDE,ISTATE,NSMIN,NSMAX,MVOLUM, + IHEX,LXI,MCODE,IPLANI,LATV) CALL LCMSIX(IPTRK,'DATA_DUP',1) CALL LCMLEN(IPTRK,'SURL_HEX',ISUR,ITPLCM) ALLOCATE(SURL(MAX(1,ISUR))) IF(ISUR.GT.0) CALL LCMGET(IPTRK,'SURL_HEX',SURL) CALL LCMSIX(IPTRK,' ',2) * NCELAP=NCELA/IPLANZ IDIM1=2*NCELA+MVOLUM IDIM2=2*IPLANZ IAUX=NV+NS+1 ALLOCATE(MESH(MMESH),ICEL(IDIM1),IPLAN(IDIM2)) ALLOCATE(VOLSUR(IAUX),MATALB(IAUX),KEYMRG(IAUX)) IF(IQUAD.EQ.1.OR.IQUAD.EQ.2) THEN ALLOCATE(XGS(NPHI),WGS(NPHI)) IF(NDIM.EQ.3) THEN ALLOCATE(XGS1(NTETA),WGS1(NTETA)) ENDIF ENDIF ALLOCATE(WW(NANGLE),A1(NANGLE),A2(NANGLE)) IF(NDIM.EQ.3)ALLOCATE(A3(NANGLE)) IF(IQUAD.EQ.1) THEN * * GAUSS-LEGENDRE INTEGRATION POINTS * CALL ALGPT(NPHI,-1.,1.,XGS,WGS) IF(NDIM.EQ.3) > CALL ALGPT(NTETA,-1.,1.,XGS1,WGS1) ELSEIF(IQUAD.EQ.2) THEN * * EQUAL WEIGHT INTEGRATION POINTS. * DO 51 IA=1,NPHI XGS(IA)=(2.0*REAL(IA)-1.0)/REAL(NPHI)-1.0 WGS(IA)=2.0/REAL(NPHI) 51 CONTINUE IF(NDIM.EQ.3) THEN DO 61 IA=1,NTETA XGS1(IA)=(2.0*REAL(IA)-1.0)/REAL(NTETA)-1.0 WGS1(IA)=2.0/REAL(NTETA) 61 CONTINUE ENDIF * ENDIF IF(IQUAD.EQ.1.OR.IQUAD.EQ.2) THEN KGAUS=-1 IF(NDIM.EQ.3) THEN AUX=PI*0.125D0*0.125D0 DO 445 IG=0,NTETA-1 * * INTEGRATION OVER (0,PI/2) (COLATITUDE) * SINT=SIN(0.25D0*PI*(DBLE(XGS1(IG+1))+1.0D0)) COST=COS(0.25D0*PI*(DBLE(XGS1(IG+1))+1.0D0)) * INTEGRATION OVER (PI/2,PI) (COLATITUDE) SINT1=SIN(0.25D0*PI*(DBLE(XGS1(IG+1))+3.0D0)) COST1=COS(0.25D0*PI*(DBLE(XGS1(IG+1))+3.0D0)) DO 444 JG=0,NPHI-1 * INTEGRATION OVER (0,PI/2) (LATITUDE) SINF=SIN(0.25D0*PI*(DBLE(XGS(JG+1))+1.0D0)) COSF=COS(0.25D0*PI*(DBLE(XGS(JG+1))+1.0D0)) KGAUS=KGAUS+1 WW(KGAUS+1)=REAL(AUX*WGS1(IG+1)*WGS(JG+1)*SINT) A1(KGAUS+1)=REAL(SINT*COSF) A2(KGAUS+1)=REAL(SINT*SINF) A3(KGAUS+1)=REAL(COST) KGAUS=KGAUS+1 WW(KGAUS+1)=REAL(AUX*WGS1(IG+1)*WGS(JG+1)*SINT1) A1(KGAUS+1)=REAL(SINT1*COSF) A2(KGAUS+1)=REAL(SINT1*SINF) A3(KGAUS+1)=REAL(COST1) * INTEGRATION OVER (PI/2,PI) SINF1=SIN(0.25D0*PI*(DBLE(XGS(JG+1))+3.0D0)) COSF1=COS(0.25D0*PI*(DBLE(XGS(JG+1))+3.0D0)) KGAUS=KGAUS+1 WW(KGAUS+1)=REAL(AUX*WGS1(IG+1)*WGS(JG+1)*SINT) A1(KGAUS+1)=REAL(SINT*COSF1) A2(KGAUS+1)=REAL(SINT*SINF1) A3(KGAUS+1)=REAL(COST) KGAUS=KGAUS+1 WW(KGAUS+1)=REAL(AUX*WGS1(IG+1)*WGS(JG+1)*SINT1) A1(KGAUS+1)=REAL(SINT1*COSF1) A2(KGAUS+1)=REAL(SINT1*SINF1) A3(KGAUS+1)=REAL(COST1) 444 CONTINUE 445 CONTINUE DEALLOCATE(XGS1,WGS1) ELSE AUX=0.125D0 DO 544 JG=0,NPHI-1 * INTEGRATION OVER (0,PI/2) SINF=SIN(0.25D0*PI*(DBLE(XGS(JG+1))+1.0D0)) COSF=COS(0.25D0*PI*(DBLE(XGS(JG+1))+1.0D0)) KGAUS=KGAUS+1 WW(KGAUS+1)=REAL(AUX*WGS(JG+1)) A1(KGAUS+1)=REAL(COSF) A2(KGAUS+1)=REAL(SINF) * INTEGRATION OVER (PI/2,PI) SINF=SIN(0.25*PI*(DBLE(XGS(JG+1))+3.0D0)) COSF=COS(0.25*PI*(DBLE(XGS(JG+1))+3.0D0)) KGAUS=KGAUS+1 WW(KGAUS+1)=REAL(AUX*WGS(JG+1)) A1(KGAUS+1)=REAL(COSF) A2(KGAUS+1)=REAL(SINF) 544 CONTINUE ENDIF DEALLOCATE(WGS,XGS) ENDIF NANGL= NANGLE ALLOCATE(ANGLES(3*NANGL),DENSTY(NANGL)) DO 21 IANG=1,NANGL ANGLES((IANG-1)*NDIM+1)=A1(IANG) ANGLES((IANG-1)*NDIM+2)=A2(IANG) IF(NDIM.EQ.3) ANGLES((IANG-1)*NDIM+3)=0.0 DENSTY(IANG)=1.0/WW(IANG) 21 CONTINUE IF(NDIM.EQ.3)ALLOCATE(IFACB(2*NCELAP)) ALLOCATE(IFFV(NCELA),ISURB(NS-NSMIN-NSMAX),IVSYM(NV),ISSYM(NS)) CALL MESHST(IPTRK,IPGEOM,MESH,ICEL,IPLAN, + IPLAN(IPLANZ+1),NCELA,IPLANZ,ISTATE, + ICEL(NCELA+1),ICEL(2*NCELA+1),NCPHY, + VOLSUR,MATALB,SIDE,NCOUR, + NSMIN,NSMAX,NS,IFACB,IFFV, + ISURB,IVSYM,ISSYM,IHEX,LXI, + NV,MCODE,SURL,IPLANI,LATV,ZMIN) CALL LCMSIX(IPGEOM,' ',0) *d - IF(IHEX.EQ.9) THEN IF(MCODE.GT.0) DEALLOCATE(LATV) ENDIF DEALLOCATE(SURL) *---- * SET ARRAY MERGE CONTAINING MERGED VOLUMES AND SURFACES *---- NSOUT=0 JAUX=0 DO 20 J=NS-1, 0,-1 KEYMRG(JAUX+1)=-ISSYM(J+1) NSOUT=MAX(NSOUT,ISSYM(J+1)) JAUX= JAUX+1 20 CONTINUE KEYMRG(JAUX+1)=0 JAUX= JAUX+1 DO 29 J=0,NV-1 KEYMRG(JAUX+1)=IVSYM(J+1) JAUX= JAUX+1 29 CONTINUE DEALLOCATE(IVSYM,ISSYM) * * OPENING THE TRACK FILE * IFILE=KDROPN('DUMMYSQ2',0,2,0) IF( IFILE.LE.0 ) GO TO 998 IF( IFILE.EQ.IFTEMP ) CALL XABORT('XHXTRK: BAD TRACKING UNIT') IF( LTRK.NE.0 )THEN MAXCYL=1 DO 140 I=0,NCELA-1 MAXCYL=MAX(MAXCYL,ICEL(NCELA+I+1)) 140 CONTINUE MAXSEC=1 SECTOR=.FALSE. DO 141 I=0,MVOLUM-1 MMAXS=ICEL(2*NCELA+I+1) IF(MMAXS.GT.1) SECTOR=.TRUE. MAXSEC=MAX(MAXSEC,MMAXS) 141 CONTINUE ITOPT= LTRK-1 NALBG= 6 NSURF=2 IF(NDIM.EQ.3) NSURF=4 MXSEG= 2*NV+2*NSURF WRITE(IFILE) '$TRK',5,0,0 COMENT='CREATOR : DRAGON' WRITE(IFILE) COMENT COMENT='MODULE : XHXTRK' WRITE(IFILE) COMENT COMENT='TYPE : HEXAGONALE' WRITE(IFILE) COMENT COMENT='GEOMETRY : '//GEONAM WRITE(IFILE) COMENT COMENT=TITREC WRITE(IFILE) COMENT WRITE(IFILE) NDIM,ITOPT,NV,NS,NALBG,NSURF,NANGL,1,MXSEG WRITE(IFILE) (VOLSUR(J),J=1,IAUX) WRITE(IFILE) (MATALB(J),J=1,IAUX) WRITE(IFILE) (ICODE(J),J=1,NALBG) WRITE(IFILE) (ALBEDO(J),J=1,NALBG) WRITE(IFILE) (ANGLES(J),J=1,NDIM*NANGL) WRITE(IFILE) (DENSTY(J),J=1,NANGL) DEALLOCATE(DENSTY,ANGLES) MACP=(MAXCYL+1)*2*NCELA*(1+6*(MAXSEC-1)) ALLOCATE(POP(MACP),IMAT(MACP*2)) MTT=2*(MAXCYL+1) MT0=2*MTT*(1+6*(MAXSEC-1)) ALLOCATE(T0(2*MT0),T1(2*MT0),IV0(MT0)) * * RAYON DU CYLINDRE ENGLOBANT LA GEOMETRIE ET VALEURES EXTREMES DE Z * DSIDE=DBLE(SIDE) RAYON=0.5D0*DSIDE*SQRT(1.0D0+3.0D0*(DBLE(2*NCOUR-1))**2) * * FIND ODD NUMBER OF TRACKS THAT ENSURES MINIMUM REQUIRED * DENSITY IS RESPECTED * NBLINE=INT(RAYON*DENUSR+0.5D0) IF(MOD(NBLINE,2) .EQ. 0) NBLINE=NBLINE+1 PASY=1.D0/DBLE(NBLINE) PASZ=PASY IF(NDIM.EQ.3) THEN NBLINE=INT(RAYON*DENUSZ+0.5D0) IF(MOD(NBLINE,2) .EQ. 0) NBLINE=NBLINE+1 PASZ=1.D0/DBLE(NBLINE) ENDIF KANG=0 ICC=1 ICC1=1 IF(NCELA.GT.1) THEN ICC=2+(NCOUR-1)*(1+3*(NCOUR-2)) ICC1=2+(NCOUR-1)*(3*(NCOUR-2)) ENDIF DZMIN=DBLE(ZMIN) IF(NDIM.EQ.3) THEN DZMAX=DBLE(MESH(3*NCELA)) ELSE DZMAX=DZMIN ENDIF * SQRT3=SQRT(3.) ALLOCATE(RAUX(2*MAXCYL)) * *--CREATION DU VECTEUR CORN NECESSAIRE POUR LE STOCKAGE DES CELLULES *--PERIPHERIQUES UTILISEES DANS TRKHEX * IF(NDIM.EQ.3) THEN ICORN=IPLAN(1) IF(IPLANZ.GT.1) ICORN=ICORN+IPLAN(1) IF(IPLANZ.GT.2) ICORN=ICORN+(IPLANZ-2)*6*(NCOUR-1) ELSE ICORN=1 IF(NCOUR.GT.1)ICORN=6*(NCOUR-1) ENDIF ALLOCATE(KCORN(ICORN)) IFCC1=2+(NCOUR-1)*(3*(NCOUR-2)) IF(NCELA.EQ.1) IFCC1=1 IXX=0 IF(NDIM.EQ.3) THEN DO 147 I=1,IPLAN(1) KCORN(IXX+1)=I IXX=IXX+1 147 CONTINUE IF(IPLANZ.GT.1) THEN IAUX=IPLAN(1)+IFCC1 DO 149 I=2,IPLANZ-1 DO 148 J=IFCC1,IPLAN(1) KCORN(IXX+1)=IAUX IAUX=IAUX+1 IXX=IXX+1 148 CONTINUE IAUX=IPLAN(I)+IFCC1 149 CONTINUE IYY=IPLANZ*IPLAN(1) IYY0=(IPLANZ-1)*IPLAN(1)+1 DO 146 I=IYY0,IYY KCORN(IXX+1)=I IXX=IXX+1 146 CONTINUE ENDIF ELSE DO 145 J=IFCC1,NCELA KCORN(IXX+1)=J IXX=IXX+1 145 CONTINUE ENDIF * * LES CELLULES VOISINES A CHAQUE VOLUME SONT STOCKEES. * NCEL2=NCELA/IPLANZ ALLOCATE(IVOIS(6*NCEL2)) IVV=-1 DO 165 I=1,NCEL2 DO 160 J=1,6 IVV=IVV+1 IVOIS(IVV+1)=NEIGHB(I,J,9,NCEL2,POIDSH) 160 CONTINUE 165 CONTINUE IF(IPRT .GE. 100) THEN IF(NDIM .EQ. 3) THEN WRITE(IOUT,6300) NAMSBR,NCELA WRITE(IOUT,6301) (MESH(II+1),MESH(II+NCELA), > MESH(II+2*NCELA),II=1,NCELA) ELSE WRITE(IOUT,6200) NAMSBR,NCELA WRITE(IOUT,6201) (MESH(II),MESH(II+NCELA),II=1,NCELA) ENDIF ENDIF *---- * SAVE EXCELL TRACKING FOR HEXAGONAL GEOMETRY *---- CALL LCMSIX(IPTRK,'EXCELL ',1) IPARAM(:NSTATE)=0 IAUX=NV+NS+1 IAUXS=NS-NSMIN-NSMAX IAUXX=2*NCELAP IPARAM(1)=NDIM IPARAM(2)=NS IPARAM(3)=NV IPARAM(4)=NCELA IPARAM(5)=MMESH IPARAM(6)=IAUX IPARAM(7)=MAXCYL IPARAM(8)=MAXSEC IPARAM(9)=IPLANZ IPARAM(10)=IDIM1 IPARAM(11)=IDIM2 IPARAM(12)=IAUXS IPARAM(13)=IAUXX CALL LCMPUT(IPTRK,'ICEL ',IDIM1 ,1,ICEL ) CALL LCMPUT(IPTRK,'IPLAN ',IDIM2 ,1,IPLAN) IF(NDIM.EQ.3) >CALL LCMPUT(IPTRK,'IFACB ',IAUXX ,1,IFACB) CALL LCMPUT(IPTRK,'IFFV ',NCELA ,1,IFFV ) IF(IAUXS.GT.0) >CALL LCMPUT(IPTRK,'ISURB ',IAUXS ,1,ISURB) CALL LCMPUT(IPTRK,'REMESH ',MMESH ,2,MESH ) CALL LCMPUT(IPTRK,'KEYMRG ',IAUX ,1,KEYMRG ) CALL LCMPUT(IPTRK,'MATALB ',IAUX ,1,MATALB ) CALL LCMPUT(IPTRK,'VOLSUR ',IAUX ,2,VOLSUR ) CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,IPARAM ) CALL LCMSIX(IPTRK,'EXCELL ',2) DO 150 IANG=1,NANGLE KANG=IANG COS1=A1(IANG) COS2=A2(IANG) COS3=0.0 IF(NDIM.EQ.3) COS3=A3(IANG) PPASY=PASY PPASZ=PASZ POIDS=WW(IANG) IF(IPRT .GE. 500) WRITE(IOUT,6400) IANG CALL TRKHEX(IPRT,NCELA,ICEL,MESH,MMESH,PPASY, + DSIDE,COS1,COS2,COS3,POP,IPLAN, + IPLANZ,IPLAN(IPLANZ+1),NDIM, + ICEL(NCELA+1),IMAT,IFILE,KANG,POIDS, + IMAT(MACP+1),ICEL(2*NCELA+1),T0, + T1,IV0,IV0(MTT+1),PPASZ,RAYON,DZMIN, + DZMAX,IFACB,IFFV,SECTOR,NSMIN+1,ISURB,RAUX, + NSURF,KCORN,ICORN,IVOIS,NCEL2,NS) ZMIN=REAL(DZMIN) 150 CONTINUE DEALLOCATE(IVOIS,KCORN,IMAT,ISURB,RAUX,POP,T0,T1,IV0,WW,A1,A2, + IFFV) IF(NDIM.EQ.3) DEALLOCATE(IFACB,A3) ENDIF DEALLOCATE(MESH,ICEL,IPLAN) *---- * SET BC MATRIX *---- ALLOCATE(MATRT(NSOUT)) DO 200 ISS=1,NSOUT MATRT(ISS)=ISS 200 CONTINUE CALL LCMPUT(IPTRK,'BC-REFL+TRAN',NSOUT,1,MATRT) DEALLOCATE(MATRT) *---- * SET NCOR=1 *---- REWIND IFILE CALL XELCOR(IFILE,IFTEMP) IER=KDRCLS(IFILE,2) IF(IER.LT.0) GO TO 999 DEALLOCATE(KEYMRG,MATALB,VOLSUR) RETURN * 998 WRITE(IOUT,'(31H ECHO = UNABLE TO OPEN FILE FT,I4)') IFILE CALL XABORT('XHXTRK: OPEN FAILED') 999 WRITE(IOUT,'(31H ECHO = UNABLE TO CLOSE FILE FT,I4)') IFILE CALL XABORT('XHXTRK: CLOSE FAILED') 6000 FORMAT(' >>> GLOBAL GEOMETRY NAME TREATED: ',A12/ > ' >>> EXACT EXCELL *CP* TREATMENT <<'/) 6200 FORMAT(' ---> OUTPUT FROM ',A6,': NUMBER OF CELLS = ',I10/ > 9X,'X-MESH',9X,'Y-MESH') 6201 FORMAT(1P,2E15.6) 6300 FORMAT(' ---> OUTPUT FROM ',A6,': NUMBER OF CELLS = ',I10/ > 9X,'X-MESH',9X,'Y-MESH',9X,'Z-MESH') 6301 FORMAT(1P,3E15.6) 6400 FORMAT(' INTEGRATION DIRECTION = ',I10) END