diff options
Diffstat (limited to 'Dragon/src/PIJS3D.f')
| -rw-r--r-- | Dragon/src/PIJS3D.f | 188 |
1 files changed, 188 insertions, 0 deletions
diff --git a/Dragon/src/PIJS3D.f b/Dragon/src/PIJS3D.f new file mode 100644 index 0000000..9f3ded5 --- /dev/null +++ b/Dragon/src/PIJS3D.f @@ -0,0 +1,188 @@ +*DECK PIJS3D + SUBROUTINE PIJS3D(NREG,NSOUT,NSLINE,WEIGHT,RCUTOF,SIGTAL, + > SEGLEN,NRSEG,STAYIN,GOSOUT,DPR) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Integration for general 3D specular 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, G. Marleau +* +*Parameters: input +* NREG total number of regions. +* NSOUT number of outer surface. +* NSLINE number of segemnts on line. +* WEIGHT line weight. +* RCUTOF MFP cut-off factor (truncate lines). +* SIGTAL albedo-cross section vector. +* SEGLEN length of track. +* NRSEG region crossed by track. +* +*Parameters: output +* DPR collision probabilities. +* +*Parameters: scratch +* STAYIN stay-in zone probability. +* GOSOUT goes-out zone probability. +* +*References: +* R. Roy et al., +* A Cyclic Tracking Procedure for CP Calculations in 2-D Lattices +* Conf/Advances in Math, Comp & Reactor Physics, +* Pittsburgh, V 1, P 2.2 4-1 (1991). +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* VARIABLES +*---- + INTEGER NREG,NSOUT,NSLINE,NRSEG(NSLINE) + REAL RCUTOF,SIGTAL(-NSOUT:NREG) + DOUBLE PRECISION WEIGHT,SEGLEN(NSLINE),STAYIN(NSLINE), + > GOSOUT(NSLINE) + DOUBLE PRECISION DPR(-NSOUT:NREG,-NSOUT:NREG) +*---- +* Local variables +*---- + INTEGER IL,JL,NOIL,NOJL,ISODD,JSODD,IJDEL + DOUBLE PRECISION TTOT,XSIL,OPATH,FINV,CUTOF + REAL ZERO,ONE,HALF + PARAMETER (ZERO=0.0E0, ONE=1.0E0, HALF=0.5E0 ) + REAL SIXT,CUTEXP + PARAMETER (SIXT=HALF/3.0,CUTEXP=0.02) + DOUBLE PRECISION EXSIL,XSIL2 + TTOT= ONE +* +*1.1) CHANGE PATHS => GOSOUT AND STAYIN PATHS, INCLUDING ALBEDOS +* ADD *PII* LOCAL NON-CYCLIC CONTRIBUTIONS + ISODD=0 + DO 30 IL= 1, NSLINE + NOIL = NRSEG(IL) + IF( NOIL.LT.0 )THEN + IF(ISODD .EQ. 1) THEN + ISODD=0 +*---- +* FOR SURFACES: +* OLD VERSION BEFORE SURFACE DOUBLING +* GOSOUT= ALBEDO * SURFACE WEIGHT +* WHERE ALL SURFACE WEIGHTS WERE 1.0 +* NEW VERSION WITH SURFACE DOUBLING +* GOSOUT= ALBEDO +* STAYIN = 1- ALBEDO * SURFACE WEIGHT +* TTOT = PRODUCT OF GOSOUT +*---- + GOSOUT(IL)= SIGTAL(NOIL) + STAYIN(IL)= ONE - GOSOUT(IL) + TTOT= TTOT * GOSOUT(IL) + ELSE + ISODD=1 +*---- +* FOR SURFACES: +* OLD VERSION BEFORE SURFACE DOUBLING +* GOSOUT= ALBEDO * SURFACE WEIGHT +* WHERE ALL SURFACE WEIGHTS WERE 1.0 +* NEW VERSION WITH SURFACE DOUBLING +* GOSOUT= ALBEDO +* STAYIN = 1- ALBEDO * SURFACE WEIGHT +* TTOT = PRODUCT OF GOSOUT +*---- + GOSOUT(IL)= SIGTAL(NOIL) + STAYIN(IL)= ONE + ENDIF + ELSE +*---- +* FOR REGIONS +* STAYIN = 1 - EXP[ -CROSS SECTION * LENGTH OF NSLINE] +* GOSOUT = 1 - STAYIN +* TTOT = PRODUCT OF GOSOUT +*---- + XSIL = SIGTAL(NOIL) + IF( XSIL .EQ. ZERO) THEN + GOSOUT(IL)= ONE + STAYIN(IL)= SEGLEN(IL) + DPR(NOIL,NOIL)= DPR(NOIL,NOIL) + > + HALF*WEIGHT*STAYIN(IL)*STAYIN(IL) + ELSE IF( XSIL .LT. CUTEXP) THEN + OPATH= SIGTAL(NOIL)*SEGLEN(IL) + XSIL2=OPATH*OPATH + EXSIL=XSIL2*(HALF-SIXT*OPATH+XSIL2/24.0) + STAYIN(IL)=OPATH-EXSIL + GOSOUT(IL)= ONE - STAYIN(IL) + TTOT= TTOT * GOSOUT(IL) + DPR(NOIL,NOIL)= DPR(NOIL,NOIL) + WEIGHT*EXSIL + ELSE + OPATH= SIGTAL(NOIL)*SEGLEN(IL) + EXSIL= EXP(-OPATH) + STAYIN(IL)= ONE - EXSIL + GOSOUT(IL)= EXSIL + TTOT= TTOT * GOSOUT(IL) + DPR(NOIL,NOIL)= DPR(NOIL,NOIL) + > + WEIGHT*(OPATH-STAYIN(IL)) + ENDIF + ENDIF + 30 CONTINUE +* +*1.2) COMPUTE CYCLIC FACTORS BY ANGLE +* USING GLOBAL TRACK ATTENUATION: BETA(TOT)*EXP(-MFP(TOT)) + IF(TTOT .GE. ONE )THEN + CALL XABORT( 'PIJS3D: ALBEDOS ARE NOT COMPATIBLE') + ENDIF + FINV= WEIGHT / (ONE-TTOT) +* +*1.3) ADD *PIJ* CONTRIBUTIONS FOR FORWARD SOURCES + ISODD=0 + DO 50 IL= 1, NSLINE + NOIL = NRSEG(IL) + TTOT= FINV * STAYIN(IL) + CUTOF= RCUTOF*TTOT + IF( NOIL .LT. 0) THEN + ISODD=MOD(ISODD+1,2) + JSODD=ISODD + DO 70 IJDEL= 1, NSLINE + JL= MOD(IL+IJDEL-1,NSLINE) + 1 + NOJL=NRSEG(JL) + IF( NOJL .LT. 0 ) THEN + JSODD=MOD(JSODD+1,2) + IF( ISODD.EQ.1 .AND. JSODD .EQ.0) THEN + DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL) + TTOT= TTOT * GOSOUT(JL) + IF( TTOT.LE.CUTOF ) GO TO 55 + ENDIF + ELSE IF(ISODD.EQ.1) THEN + DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL) + TTOT= TTOT * GOSOUT(JL) + IF( TTOT.LE.CUTOF ) GO TO 55 + ENDIF + 70 CONTINUE + ELSE + JSODD=ISODD + DO 80 IJDEL= 1, NSLINE + JL= MOD(IL+IJDEL-1,NSLINE) + 1 + NOJL=NRSEG(JL) + IF( NOJL .LT. 0 ) THEN + JSODD=MOD(JSODD+1,2) + IF( JSODD .EQ.0) THEN + DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL) + TTOT= TTOT * GOSOUT(JL) + IF( TTOT.LE.CUTOF ) GO TO 55 + ENDIF + ELSE + DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL) + TTOT= TTOT * GOSOUT(JL) + IF( TTOT.LE.CUTOF ) GO TO 55 + ENDIF + 80 CONTINUE + ENDIF + 55 CONTINUE + 50 CONTINUE + RETURN + END |