diff options
Diffstat (limited to 'Dragon/src/QIJI3D.f')
| -rw-r--r-- | Dragon/src/QIJI3D.f | 188 |
1 files changed, 188 insertions, 0 deletions
diff --git a/Dragon/src/QIJI3D.f b/Dragon/src/QIJI3D.f new file mode 100644 index 0000000..59dde2c --- /dev/null +++ b/Dragon/src/QIJI3D.f @@ -0,0 +1,188 @@ +*DECK QIJI3D + SUBROUTINE QIJI3D(NREG,NSOUT,NPIJ,NGRP,MXSEG,NCOR,SWVOID,LINE, + > WEIGHT,NUMERO,LENGHT,SIGTAL,NPSYS,DPR) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Integration for general 3D isotropic tracking. +* +*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): R. Roy +* +*Parameters: input +* NREG total number of regions. +* NSOUT number of outer surface. +* NPIJ number of probabilities in one group. +* NGRP number of groups. +* MXSEG number of segemnts on line. +* NCOR maximum number of corners. +* SWVOID flag to indicate if there are voids. +* LINE line number. +* WEIGHT line weight. +* NUMERO region crossed by track. +* LENGHT length of track. +* SIGTAL albedo-cross section vector. +* NPSYS non-converged energy group indices. +* +*Parameters: output +* DPR CP matrix. +* +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE + INTEGER NREG,NSOUT,NPIJ,NGRP,MXSEG,NCOR, NUMERO(*), + > ISD(6),ISF(6),LIN2C,IL,JL,IG,LINE,NOIL,IND1, + > IND2,I,J,IN0,IN1,IN2,NPSYS(NGRP),NUNK + REAL WEIGHT, SIGTAL(-NSOUT:NREG,NGRP) + DOUBLE PRECISION LENGHT(*), DPR(NPIJ,NGRP), XSIL, XSIL2 + LOGICAL SWVOID + REAL ZERO, ONE, HALF + PARAMETER ( ZERO=0.0E0, ONE=1.0E0, HALF=0.5E0) + REAL SIXT,CUTEXP + PARAMETER ( CUTEXP=0.02) +* +* Allocated arrays + INTEGER, ALLOCATABLE, DIMENSION(:) :: NRSEG + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DSCBEG, DSCEND, + > SEGLEN, PRODUC + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: STAYIN, GOSOUT + + IN0(I) = ((I+NSOUT+1)*(I+NSOUT+2))/2 + IN1(I,J) = ((I+NSOUT+1)*(I+NSOUT))/2 + (J+NSOUT+1) + IN2(I,J) = (MAX(I+NSOUT+1,J+NSOUT+1)* + > (MAX(I+NSOUT+1,J+NSOUT+1)-1))/2 + > + MIN(I+NSOUT+1,J+NSOUT+1) +* +* Scratch storage allocation + ALLOCATE(NRSEG(MXSEG)) + ALLOCATE(DSCBEG(NGRP), DSCEND(NGRP), SEGLEN(MXSEG), PRODUC(NGRP)) + ALLOCATE(STAYIN(NGRP,MXSEG),GOSOUT(NGRP,MXSEG)) +* + SIXT=HALF/3.0 + NUNK=NSOUT+NREG+1 +* +* 0) REFORMAT TRACKING LINE + LIN2C= LINE-2*NCOR + DO 10 JL= 1, NCOR + ISD(JL)= NUMERO(JL) + ISF(JL)= NUMERO(NCOR+LIN2C+JL) + 10 CONTINUE + DO 20 IL= 1, LIN2C + NRSEG(IL)= NUMERO(NCOR+IL) + SEGLEN(IL)= LENGHT(NCOR+IL) + 20 CONTINUE + DO 90 IG= 1, NGRP + PRODUC(IG)= WEIGHT + 90 CONTINUE +* + IF( SWVOID )THEN +* +* 1) VOIDS ARE POSSIBLE +* PII CALCULATION AND ESCAPE + DO 240 IL = 1,LINE-2*NCOR + NOIL = NRSEG(IL) + IND1 = IN0(NOIL) + DO 100 IG= 1, NGRP + IF(NPSYS(IG).EQ.0) GO TO 100 + XSIL=SIGTAL(NOIL,IG)*SEGLEN(IL) + IF( XSIL.EQ.ZERO )THEN + GOSOUT(IG,IL)= ONE + STAYIN(IG,IL)= SEGLEN(IL) + DPR(IND1,IG)= DPR(IND1,IG) + > + HALF*WEIGHT*SEGLEN(IL)*SEGLEN(IL) + ELSE IF(XSIL .LT. CUTEXP) THEN + XSIL2=XSIL*XSIL + STAYIN(IG,IL)=XSIL-XSIL2*(HALF-SIXT*XSIL) + GOSOUT(IG,IL)=ONE-STAYIN(IG,IL) + PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL) + DPR(IND1,IG)= DPR(IND1,IG) + > + WEIGHT*XSIL2*(HALF-SIXT*XSIL) + ELSE + GOSOUT(IG,IL)= EXP( -XSIL ) + STAYIN(IG,IL)= (ONE - GOSOUT(IG,IL)) + PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL) + DPR(IND1,IG)= DPR(IND1,IG) + > + WEIGHT*(XSIL-STAYIN(IG,IL)) + ENDIF + 100 CONTINUE + 240 CONTINUE + ELSE + DO 241 IL = 1,LINE-2*NCOR + NOIL = NRSEG(IL) + IND1 = IN0(NOIL) + DO 101 IG= 1, NGRP + IF(NPSYS(IG).EQ.0) GO TO 101 + XSIL=SIGTAL(NOIL,IG)*SEGLEN(IL) + IF(XSIL .LT. CUTEXP) THEN + XSIL2=XSIL*XSIL + STAYIN(IG,IL)=XSIL-XSIL2*(HALF-SIXT*XSIL) + GOSOUT(IG,IL)=ONE-STAYIN(IG,IL) + PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL) + DPR(IND1,IG)= DPR(IND1,IG) + > + WEIGHT*XSIL2*(HALF-SIXT*XSIL) + ELSE + GOSOUT(IG,IL)= EXP( -XSIL ) + STAYIN(IG,IL)= (ONE - GOSOUT(IG,IL)) + PRODUC(IG)= PRODUC(IG) * GOSOUT(IG,IL) + DPR(IND1,IG)= DPR(IND1,IG) + > + WEIGHT*(XSIL-STAYIN(IG,IL)) + ENDIF + 101 CONTINUE + 241 CONTINUE + ENDIF +* PIJ CALCULATION + DO 120 IG= 1, NGRP + DSCBEG(IG)= WEIGHT + 120 CONTINUE + DO 260 IL = 1, LINE-2*NCOR + NOIL = NRSEG(IL) + DO 130 IG= 1, NGRP + IF(NPSYS(IG).NE.0) DSCEND(IG)= WEIGHT*STAYIN(IG,IL) + 130 CONTINUE + DO 250 JL = IL+1, LINE-2*NCOR + IND2= IN2(NRSEG(JL),NOIL) + DO 140 IG= 1, NGRP + IF(NPSYS(IG).EQ.0) GO TO 140 + DPR(IND2,IG)= DPR(IND2,IG) + STAYIN(IG,JL)*DSCEND(IG) + DSCEND(IG)= DSCEND(IG)*GOSOUT(IG,JL) + 140 CONTINUE + 250 CONTINUE +* PIS CALCULATION + DO 261 JL = 1, NCOR + IND1= IN1(NOIL,ISD(JL)) + IND2= IN1(NOIL,ISF(JL)) + DO 150 IG= 1, NGRP + IF(NPSYS(IG).EQ.0) GO TO 150 + DPR(IND1,IG)= DPR(IND1,IG) + DSCBEG(IG)*STAYIN(IG,IL) + DPR(IND2,IG)= DPR(IND2,IG) + DSCEND(IG) + 150 CONTINUE + 261 CONTINUE + DO 160 IG= 1, NGRP + IF(NPSYS(IG).NE.0) DSCBEG(IG)= DSCBEG(IG)*GOSOUT(IG,IL) + 160 CONTINUE + 260 CONTINUE +* PSS CALCULATION + DO 265 IL = 1, NCOR + DO 264 JL = 1, NCOR + IND2= IN2(ISD(IL),ISF(JL)) + DO 170 IG = 1, NGRP + IF(NPSYS(IG).NE.0) DPR(IND2,IG)= DPR(IND2,IG) + PRODUC(IG) + 170 CONTINUE + 264 CONTINUE + 265 CONTINUE +* +* Scratch storage deallocation + DEALLOCATE(GOSOUT,STAYIN) + DEALLOCATE(PRODUC,SEGLEN,DSCEND,DSCBEG) + DEALLOCATE(NRSEG) +* + RETURN + END |