diff options
Diffstat (limited to 'Trivac/src/READ3D.f')
| -rwxr-xr-x | Trivac/src/READ3D.f | 420 |
1 files changed, 420 insertions, 0 deletions
diff --git a/Trivac/src/READ3D.f b/Trivac/src/READ3D.f new file mode 100755 index 0000000..ba0385d --- /dev/null +++ b/Trivac/src/READ3D.f @@ -0,0 +1,420 @@ +*DECK READ3D + SUBROUTINE READ3D (MAXX,MAXY,MAXZ,MAXPTS,IPGEOM,IHEX,IR,ILK,SIDE, + 1 XXX,YYY,ZZZ,IMPX,LX,LY,LZ,MAT,NMBLK,NCODE,ICODE,ZCODE,ISPLTX, + 2 ISPLTY,ISPLTZ,ISPLTH,ISPLTL) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Recover the input data for the description of a 1-D, 2-D or 3-D +* Cartesian, cylindrical, spherical or hexagonal domain. +* +*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/output +* MAXX allocated storage for arrays of dimension LX. +* MAXY allocated storage for arrays of dimension LY. +* MAXZ allocated storage for arrays of dimension LZ. +* MAXPTS allocated storage for arrays of dimension NMBLK. +* IPGEOM L_GEOM pointer to the geometry. +* IHEX type of hexagonal geometry (=0 for non-hexagonal geometry). +* IR number of mixtures. +* ILK (ILK=.true. if neutron leakage through external boundary +* is present). +* SIDE side of the hexagons. XXX and YYY arrays are not used with +* hexagonal geometry. +* XXX Cartesian coordinates of the domain along the X-axis. +* YYY Cartesian coordinates of the domain along the Y-axis. +* ZZZ Cartesian coordinates of the domain along the Z-axis. +* IMPX print flag. Minimum printing if IMPX=0. +* LX number of elements along the X-axis after mesh-splitting +* or number of hexagons in one axial plane. +* LY number of elements along the Y-axis. +* LZ number of elements along the Z-axis. +* MAT index-number of the mixture type assigned to each volume +* after mesh-splitting. +* NMBLK number of elements in the domain. +* NCODE boundary condition relative to each side of the domain: +* =1: VOID ; =2: REFL ; =3: DIAG ; =4: TRAN ; =5: SYME +* =6: ALBE ; =7: ZERO ; =20: CYLI. +* ICODE physical albedo index on each side of the domain. +* ZCODE albedo relative to each side of the domain. +* ISPLTX mesh-splitting data for parallelepipeds along the X-axis +* negative value is used for equal-volume splitting of tubes. +* ISPLTY mesh-splitting data for parallelepipeds along the Y-axis. +* ISPLTZ mesh-splitting data for parallelepipeds along the Z-axis. +* ISPLTH mesh-splitting index for hexagons into triangles. +* ISPLTL mesh-splitting index for hexagons into lozenges. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPGEOM + INTEGER MAXX,MAXY,MAXZ,MAXPTS,IHEX,IR,IMPX,LX,LY,LZ,MAT(MAXPTS), + 1 NMBLK,NCODE(6),ICODE(6),ISPLTX(MAXX),ISPLTY(MAXY),ISPLTZ(MAXZ), + 2 ISPLTH,ISPLTL + REAL SIDE,XXX(MAXX+1),YYY(MAXY+1),ZZZ(MAXZ+1),ZCODE(6) + LOGICAL ILK + INTEGER, ALLOCATABLE, DIMENSION(:) :: DPP,MX +*---- +* LOCAL VARIABLES +*---- + PARAMETER(NSTATE=40) + LOGICAL LL1,LL2,LCYL,SWCEN,EMPTY,LCM + CHARACTER HSMG*131,GEONAM*12,TEXT12*12 + INTEGER ISTATE(NSTATE) + EQUIVALENCE (ITYPE,ISTATE(1)),(LR1,ISTATE(2)),(LX1,ISTATE(3)), + 1 (LY1,ISTATE(4)),(LZ1,ISTATE(5)) +* + IHEX=0 + CALL LCMGET(IPGEOM,'STATE-VECTOR',ISTATE) + IF((ITYPE.EQ.8).OR.(ITYPE.EQ.9)) THEN + CALL LCMLEN(IPGEOM,'IHEX',ILEN,ITYLCM) + IF(ILEN.EQ.0) CALL LCMLIB(IPGEOM) + IF(ILEN.EQ.0) CALL XABORT('READ3D: MISSING IHEX RECORD.') + CALL LCMGET(IPGEOM,'IHEX',IHEX) + ENDIF + IF((ISTATE(8).NE.0).OR.(ISTATE(9).NE.0).OR.(ISTATE(10).NE.0).OR. + 1 (ISTATE(13).NE.0)) CALL XABORT('READ3D: UNABLE TO PROCESS THE G' + 2 //'EOMETRY.') + LCYL=(ITYPE.EQ.3).OR.(ITYPE.EQ.4).OR.(ITYPE.EQ.6) + IDIM=1 + IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) IDIM=2 + IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) IDIM=3 +*---- +* RECOVER THE BOUNDARY CONDITIONS +*---- + CALL LCMGET(IPGEOM,'NCODE',NCODE) + CALL LCMGET(IPGEOM,'ZCODE',ZCODE) + CALL LCMGET(IPGEOM,'ICODE',ICODE) + DO 10 I=1,6 + IF(NCODE(I).EQ.10) NCODE(I)=2 + IF(NCODE(I).EQ.2) ZCODE(I)=1.0 + IF(NCODE(I).EQ.6) NCODE(I)=1 + IF((NCODE(I).EQ.20).AND.(ITYPE.NE.5).AND.(ITYPE.NE.7)) CALL + 1 XABORT('READ3D: CYLINDRICAL CORRECTION IS LIMITED TO CARTESIAN ' + 2 //'GEOMETRIES.') + IF((NCODE(I).GE.8).AND.(NCODE(I).NE.20)) THEN + CALL XABORT('READ3D: INVALID TYPE OF B.C.') + ENDIF + 10 CONTINUE +*---- +* CHECK COHERENCE OF THE CYLINDRICAL EXTERNAL B.C. +*---- + SWCEN=.FALSE. + ALBMAX=-1.0E35 + ALBMIN=+1.0E35 + DO 15 IC=1,6 + IF(NCODE(IC).NE.20) GO TO 15 + SWCEN=.TRUE. + IF(ZCODE(IC).LT.ALBMIN) ALBMIN=ZCODE(IC) + IF(ZCODE(IC).GT.ALBMAX) ALBMAX=ZCODE(IC) + 15 CONTINUE + IF(SWCEN.AND.(ALBMIN.NE.ALBMAX)) CALL XABORT('READ3D: CYLINDRICA' + 1 //'L IMBEDDED EXTERNAL GEOMETRY: ALBEDOS ARE INCONSISTENT.') +* + IF(ITYPE.GE.10) THEN + CALL XABORT('READ3D: INVALID TYPE OF GEOMETRY.') + ELSE IF(ITYPE.GE.8) THEN + IF((NCODE(2).NE.0).OR.(NCODE(3).NE.0).OR.(NCODE(4).NE.0)) + 1 CALL XABORT('READ3D: INVALID TYPE OF HEXAGONAL B.C.') + IF(NCODE(1).EQ.5) THEN + IF(IHEX.EQ.1) THEN + IHEX=10 + ELSE IF(IHEX.EQ.2) THEN + IHEX=11 + ELSE + CALL XABORT('READ3D: BOUNDARY CONDITION HBC WITH OPTION' + 1 //' SYME IS ONLY PERMITTED WITH S30 OR SA60 SYMMETRY.') + ENDIF + ELSE IF((NCODE(1).GT.2).AND.(NCODE(1).NE.7)) THEN + CALL XABORT('READ3D: BOUNDARY CONDITION HBC CAN ONLY BE US' + 1 //'ED WITH OPTIONS VOID, REFL, SYME, ALBE OR ZERO.') + ENDIF + ENDIF +*---- +* RECOVER THE MIXTURE NUMBERS +*---- + IF(ISTATE(6).GT.MAXPTS) THEN + WRITE (HSMG,690) 'NMBLK',ISTATE(6),'MAXPTS',MAXPTS + CALL XABORT(HSMG) + ENDIF + CALL LCMGET(IPGEOM,'MIX',MAT) + IR=0 + DO 20 I=1,ISTATE(6) + IR=MAX(IR,MAT(I)) + 20 CONTINUE +*---- +* RECOVER THE MESH COORDINATES. +*---- + IF(LCYL.AND.(LR1.GT.MAXX)) THEN + WRITE (HSMG,690) 'LX',LR1,'MAXX',MAXX + CALL XABORT(HSMG) + ELSE IF(LX1.GT.MAXX) THEN + WRITE (HSMG,690) 'LX',LX1,'MAXX',MAXX + CALL XABORT(HSMG) + ENDIF + IF(LY1.GT.MAXY) THEN + WRITE (HSMG,690) 'LY',LY1,'MAXY',MAXY + CALL XABORT(HSMG) + ENDIF + IF(LZ1.GT.MAXZ) THEN + WRITE (HSMG,690) 'LZ',LZ1,'MAXZ',MAXZ + CALL XABORT(HSMG) + ENDIF + LL1=.FALSE. + LL2=.FALSE. + LY=1 + YYY(1)=0.0 + YYY(2)=1.0 + LZ=1 + ZZZ(1)=0.0 + ZZZ(2)=1.0 + IF(ITYPE.EQ.2) THEN +* 1-D CARTESIAN GEOMETRY. + LX=LX1 + NMBLK=LX + IF((NCODE(1).EQ.0).OR.(NCODE(2).EQ.0)) GO TO 610 + CALL LCMGET(IPGEOM,'MESHX',XXX) + ELSE IF((ITYPE.EQ.3).OR.(ITYPE.EQ.4)) THEN +* 1-D CYLINDRICAL/SPHERICAL GEOMETRY. + LX=LR1 + NMBLK=LX + IF(NCODE(1).NE.0) GO TO 640 + IF(NCODE(2).EQ.0) GO TO 610 + NCODE(1)=2 + CALL LCMGET(IPGEOM,'RADIUS',XXX) + ELSE IF(ITYPE.EQ.5) THEN +* 2-D CARTESIAN GEOMETRY. + LX=LX1 + LY=LY1 + NMBLK=LX*LY + I2=0 + DO 30 IC=1,4 + IF(NCODE(IC).EQ.0) GO TO 610 + IF(NCODE(IC).EQ.3) I2=I2+1 + 30 CONTINUE + IF(I2.NE.0) THEN + IF((I2.NE.2).OR.(LX.NE.LY)) GO TO 630 + NMBLK=(LX+1)*LX/2 + LL1=(NCODE(2).EQ.3).AND.(NCODE(3).EQ.3) + LL2=(NCODE(1).EQ.3).AND.(NCODE(4).EQ.3) + IF((.NOT.LL1).AND.(.NOT.LL2)) GO TO 620 + ENDIF + CALL LCMGET(IPGEOM,'MESHX',XXX) + IF(LL1.OR.LL2) THEN + CALL LCMGET(IPGEOM,'MESHX',YYY) + ELSE + CALL LCMGET(IPGEOM,'MESHY',YYY) + ENDIF + ELSE IF(ITYPE.EQ.6) THEN +* 2-D CYLINDRICAL GEOMETRY. + LX=LR1 + LZ=LZ1 + NMBLK=LX*LZ + IF(NCODE(1).NE.0) GO TO 650 + IF((NCODE(2).EQ.3).OR.(NCODE(3).EQ.3).OR.(NCODE(4).EQ.3)) + 1 GO TO 660 + IF((NCODE(2).EQ.0).OR.(NCODE(5).EQ.0).OR.(NCODE(6).EQ.0)) + 1 GO TO 610 + NCODE(1)=2 + CALL LCMGET(IPGEOM,'RADIUS',XXX) + CALL LCMGET(IPGEOM,'MESHZ',ZZZ) + ELSE IF(ITYPE.EQ.7) THEN +* 3-D CARTESIAN GEOMETRY. + LX=LX1 + LY=LY1 + LZ=LZ1 + NMBLK=LX*LY*LZ + I2=0 + DO 40 IC=1,4 + IF(NCODE(IC).EQ.0) GO TO 610 + IF(NCODE(IC).EQ.3) I2=I2+1 + 40 CONTINUE + IF(I2.NE.0) THEN + IF((I2.NE.2).OR.(LX.NE.LY)) GO TO 630 + NMBLK=((LX+1)*LX/2)*LZ + LL1=(NCODE(2).EQ.3).AND.(NCODE(3).EQ.3) + LL2=(NCODE(1).EQ.3).AND.(NCODE(4).EQ.3) + IF((.NOT.LL1).AND.(.NOT.LL2)) GO TO 620 + ENDIF + CALL LCMGET(IPGEOM,'MESHX',XXX) + IF(LL1.OR.LL2) THEN + CALL LCMGET(IPGEOM,'MESHX',YYY) + ELSE + CALL LCMGET(IPGEOM,'MESHY',YYY) + ENDIF + CALL LCMGET(IPGEOM,'MESHZ',ZZZ) + ELSE IF(ITYPE.EQ.8) THEN +* 2-D HEXAGONAL GEOMETRY. + LX=LX1 + NMBLK=LX + CALL LCMGET(IPGEOM,'SIDE',SIDE) + ELSE IF(ITYPE.EQ.9) THEN +* 3-D HEXAGONAL GEOMETRY. + LX=LX1 + LZ=LZ1 + NMBLK=LX*LZ + CALL LCMGET(IPGEOM,'SIDE',SIDE) + CALL LCMGET(IPGEOM,'MESHZ',ZZZ) + ELSE + CALL XABORT('READ3D: INVALID TYPE OF GEOMETRY.') + ENDIF + IF(NMBLK.NE.ISTATE(6)) THEN + WRITE(HSMG,'(45HREAD3D: INVALID NUMBER OF REGIONS. NUMBER OF , + 1 13HMIX ENTRIES =,I7,20H NUMBER OF REGIONS =,I7)') ISTATE(6), + 2 NMBLK + CALL XABORT(HSMG) + ENDIF + DO 50 IC=1,6,2 + IF((NCODE(IC).EQ.4).AND.(NCODE(IC+1).NE.4)) GO TO 670 + 50 CONTINUE + IF((NCODE(2).EQ.3).AND.(NCODE(3).EQ.3)) THEN + ZCODE(3)=ZCODE(1) + ZCODE(2)=ZCODE(4) + ELSE IF((NCODE(1).EQ.3).AND.(NCODE(4).EQ.3)) THEN + ZCODE(1)=ZCODE(3) + ZCODE(4)=ZCODE(2) + ENDIF +*---- +* UNFOLD GEOMETRY IF HEXAGONAL IN LOZENGES +*---- + ISPLTL=0 + ISPLTH=0 + CALL LCMLEN(IPGEOM,'SPLITL',ILEN,ITYLCM) + IF(ILEN.GT.0) CALL LCMGET(IPGEOM,'SPLITL',ISPLTL) + CALL LCMLEN(IPGEOM,'SPLITH',ILEN,ITYLCM) + IF(ILEN.GT.0) CALL LCMGET(IPGEOM,'SPLITH',ISPLTH) + IF((ISPLTL.GT.0).AND.(IHEX.NE.9)) THEN + ALLOCATE(DPP(MAXPTS),MX(LX*LZ)) + DO 60 I=1,LX*LZ + MX(I)=MAT(I) + 60 CONTINUE + LXOLD=LX + CALL BIVALL(MAXPTS,IHEX,LXOLD,LX,DPP) + DO 80 KZ=1,LZ + DO 70 KX=1,LX + KEL=DPP(KX)+(KZ-1)*LXOLD + MAT(KX+(KZ-1)*LX)=MX(KEL) + 70 CONTINUE + 80 CONTINUE + DEALLOCATE(DPP,MX) + IHEX=9 + ENDIF +*---- +* MESH-SPLITTING +*---- + IF(ISTATE(11).NE.0) THEN + CALL LCMLEN(IPGEOM,'SPLITR',ILEN1,ITYLCM) + CALL LCMLEN(IPGEOM,'SPLITX',ILEN2,ITYLCM) + IF(LCYL.AND.(ILEN1.GT.0)) THEN + CALL LCMGET(IPGEOM,'SPLITR',ISPLTX) + ELSE IF(ILEN2.GT.0) THEN + CALL LCMGET(IPGEOM,'SPLITX',ISPLTX) + ELSE IF(ITYPE.LE.7) THEN + DO 90 I=1,LX + ISPLTX(I)=1 + 90 CONTINUE + ENDIF + CALL LCMLEN(IPGEOM,'SPLITY',ILEN,ITYLCM) + IF(ILEN.GT.0) THEN + CALL LCMGET(IPGEOM,'SPLITY',ISPLTY) + ELSE IF(LL1.OR.LL2) THEN + DO 100 I=1,LX + ISPLTY(I)=ISPLTX(I) + 100 CONTINUE + ELSE + DO 110 I=1,LY + ISPLTY(I)=1 + 110 CONTINUE + ENDIF + CALL LCMLEN(IPGEOM,'SPLITZ',ILEN,ITYLCM) + IF(ILEN.GT.0) THEN + CALL LCMGET(IPGEOM,'SPLITZ',ISPLTZ) + ELSE + DO 120 I=1,LZ + ISPLTZ(I)=1 + 120 CONTINUE + ENDIF + IF((ISPLTH.GT.0).AND.(ISPLTL.GT.0)) THEN + CALL XABORT('READ3D: SPLITH AND SPLITL KEYWORDS ARE EXCLUS' + 1 //'IVE.') + ENDIF + CALL SPLIT0(MAXPTS,ITYPE,NCODE,LXOLD,LYOLD,LZOLD,ISPLTX,ISPLTY, + 1 ISPLTZ,0,ISPLTL,NMBLK,LX,LY,LZ,SIDE,XXX,YYY,ZZZ,MAT,.TRUE., + 2 IMPX) + IF(NMBLK.GT.MAXPTS) THEN + WRITE (HSMG,690) 'NMBLK',NMBLK,'MAXPTS',MAXPTS + CALL XABORT(HSMG) + ENDIF + ENDIF +* + ILK=((NCODE(1).EQ.1).AND.(ZCODE(1).NE.1.0)).OR.(NCODE(1).EQ.7).OR. + 1 ((NCODE(2).EQ.1).AND.(ZCODE(2).NE.1.0)).OR.(NCODE(2).EQ.7).OR. + 2 ((NCODE(3).EQ.1).AND.(ZCODE(3).NE.1.0)).OR.(NCODE(3).EQ.7).OR. + 3 ((NCODE(4).EQ.1).AND.(ZCODE(4).NE.1.0)).OR.(NCODE(4).EQ.7).OR. + 4 ((NCODE(5).EQ.1).AND.(ZCODE(5).NE.1.0)).OR.(NCODE(5).EQ.7).OR. + 5 ((NCODE(6).EQ.1).AND.(ZCODE(6).NE.1.0)).OR.(NCODE(6).EQ.7).OR. + 6 ((NCODE(1).EQ.8).AND.(ZCODE(1).NE.1.0)).OR. + 7 ((NCODE(2).EQ.8).AND.(ZCODE(2).NE.1.0)).OR. + 8 ((NCODE(3).EQ.8).AND.(ZCODE(3).NE.1.0)).OR. + 9 ((NCODE(4).EQ.8).AND.(ZCODE(4).NE.1.0)).OR. + 1 ((NCODE(5).EQ.8).AND.(ZCODE(5).NE.1.0)).OR. + 2 ((NCODE(6).EQ.8).AND.(ZCODE(6).NE.1.0)) + IF(IMPX.GT.0) THEN + IF(ITYPE.EQ.2) THEN + WRITE (6,'(/19H 1-D SLAB GEOMETRY.)') + ELSE IF(ITYPE.EQ.3) THEN + WRITE (6,'(/26H 1-D CYLINDRICAL GEOMETRY.)') + ELSE IF(ITYPE.EQ.4) THEN + WRITE (6,'(/24H 1-D SPHERICAL GEOMETRY.)') + ELSE IF(ITYPE.EQ.5) THEN + WRITE (6,'(/24H 2-D CARTESIAN GEOMETRY.)') + ELSE IF(ITYPE.EQ.6) THEN + WRITE (6,'(/18H 2-D R-Z GEOMETRY.)') + ELSE IF(ITYPE.EQ.7) THEN + WRITE (6,'(/24H 3-D CARTESIAN GEOMETRY.)') + ELSE IF(ITYPE.EQ.8) THEN + WRITE (6,'(/24H 2-D HEXAGONAL GEOMETRY.)') + ELSE IF(ITYPE.EQ.9) THEN + WRITE (6,'(/24H 3-D HEXAGONAL GEOMETRY.)') + ENDIF + CALL LCMINF(IPGEOM,GEONAM,TEXT12,EMPTY,ILONG,LCM) + WRITE (6,'(1H+,26X,18HBASED ON GEOMETRY ,A12,1H./)') GEONAM + WRITE (6,770) LX,MAXX,LY,MAXY,LZ,MAXZ,IR + IF(.NOT.ILK) WRITE (6,'(17H INFINITE DOMAIN./)') + ENDIF + RETURN +* + 610 CALL XABORT('READ3D: A BOUNDARY CONDITION IS MISSING.') + 620 CALL XABORT('READ3D: THE DIAGONAL CONDITIONS X+: DIAG Y-: DIAG A' + 1 //'ND X-: DIAG Y+: DIAG ARE THE ONLY PERMITTED.') + 630 CALL XABORT('READ3D: LX=LY WITH A DIAGONAL SYMMETRY.') + 640 CALL XABORT('READ3D: CYLINDRICAL GEOMETRY - ONLY THE R+: BOUNDAR' + 1 //'Y CONDITION IS REQUIRED.') + 650 CALL XABORT('READ3D: CYLINDRICAL GEOMETRY - ONLY THE R+:, Z-: AN' + 1 //'D Z+: BOUNDARY CONDITIONS ARE REQUIRED.') + 660 CALL XABORT('READ3D: CYLINDRICAL GEOMETRY : THE DIAG BOUNDARY CO' + 1 //'NDITION CANNOT BE USED.') + 670 CALL XABORT('READ3D: THE TRANSLATION CONDITIONS X-: TRAN X+: TRA' + 1 //'N, Y-: TRAN Y+: TRAN AND Z-: TRAN Z+: TRAN ARE THE ONLY PERM' + 1 //'ITTED.') +* + 690 FORMAT (29HREAD3D: INSUFFICIENT STORAGE.,5X,A6,1H=,I7,8H ; AVAIL, + 1 13HABLE STORAGE ,A6,1H=,I7) + 770 FORMAT (/44H NUMBER OF MESH INTERVALS ALONG THE X AXIS =,I5,5X, + 1 24HAVAILABLE STORAGE MAXX =,I7/26X,18HALONG THE Y AXIS =,I5,5X, + 2 24HAVAILABLE STORAGE MAXY =,I7/26X,18HALONG THE Z AXIS =,I5,5X, + 3 24HAVAILABLE STORAGE MAXZ =,I7/28H NUMBER OF DISTINCT MIXTURES, + 4 2H =,I7/) + END |