diff options
Diffstat (limited to 'Dragon/src/XHXTRK.f')
| -rw-r--r-- | Dragon/src/XHXTRK.f | 646 |
1 files changed, 646 insertions, 0 deletions
diff --git a/Dragon/src/XHXTRK.f b/Dragon/src/XHXTRK.f new file mode 100644 index 0000000..7f6f8ac --- /dev/null +++ b/Dragon/src/XHXTRK.f @@ -0,0 +1,646 @@ +*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 |