summaryrefslogtreecommitdiff
path: root/Dragon/src/MOCFFIT.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MOCFFIT.f')
-rw-r--r--Dragon/src/MOCFFIT.f129
1 files changed, 129 insertions, 0 deletions
diff --git a/Dragon/src/MOCFFIT.f b/Dragon/src/MOCFFIT.f
new file mode 100644
index 0000000..3b8f669
--- /dev/null
+++ b/Dragon/src/MOCFFIT.f
@@ -0,0 +1,129 @@
+*DECK MOCFFIT
+ SUBROUTINE MOCFFIT(SUBSCH,NR,NS,NUN,MT,LINE,SEGLEN,NRSEG,NE,
+ 1 MATALB,SIGANG,KEYFLX,YG,FLUX,EXPT,EXP2,FLM,
+ 2 FLP,CYM,CYP,IDIR,OMG2)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of transport equation on a cyclic track
+* (isotropic scattering case, 'MOCC/MCI' integration strategy).
+*
+*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): R. Roy and R. Le Tellier
+*
+*Parameters: input
+* SUBSCH track coefficients calculation subroutine.
+* NR number of volumes.
+* NS number of surfaces.
+* NUN total number of unknowns in vectors FLUX.
+* MT number of material mixtures.
+* LINE number of segments on this tracking line.
+* SEGLEN vector containing the lenght of the different segments of this
+* track.
+* NRSEG vector containing the region number of the different segments
+* of this track.
+* NE order of the polar quadrature set.
+* MATALB index-number of the mixture assigned to each volume
+* and the albedo to each surface.
+* SIGANG total cross-sections and albedos.
+* KEYFLX position of flux elements in FLUX vector.
+* YG inverse of polar quadrature cosines.
+* FLM total source vector for + direction.
+* FLP total source vector for - direction.
+* IDIR direction of fundamental current for TIBERE with MoC
+* =0,1,2,3.
+* OMG2 x, y and z components of the $3\\Omega^2$.
+*
+*Parameters: output
+* FLUX vector containing the zonal flux moments.
+*
+*Parameters: scratch
+* EXPT track coefficient.
+* EXP2 quadratic expansion of (1-exp(-a*L))/L with small argument.
+* CYM undefined.
+* CYP undefined.
+*
+*-----------------------------------------------------------------------
+*
+ IMPLICIT NONE
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NR,NS,NUN,MT,LINE,NRSEG(LINE),NE,MATALB(-NS:NR),KEYFLX(NR)
+ INTEGER IDIR
+ REAL SIGANG(-6:MT),YG(NE)
+ DOUBLE PRECISION SEGLEN(LINE),FLUX(NUN),EXPT(NE,LINE),
+ 1 EXP2(NE,LINE),FLM(NE,LINE),FLP(NE,LINE),CYM(NE,LINE),
+ 2 CYP(NE,LINE),OMG2(NE,3)
+ EXTERNAL SUBSCH
+*----
+* LOCAL VARIABLES
+*----
+ INTEGER MXE
+ PARAMETER (MXE=64)
+ DOUBLE PRECISION ONE,ZERO
+ PARAMETER (ONE=1.D0,ZERO=0.D0)
+ INTEGER IE,NOIL,IL,JL,IND,INDC
+ DOUBLE PRECISION FL0(MXE),CY0(MXE),FL1(MXE),CY1(MXE)
+*
+ FL0(:NE)= ZERO
+ FL1(:NE)= ZERO
+ CY0(:NE)= ONE
+ CY1(:NE)= ONE
+*----
+* Calculation of the coefficients for this track.
+*----
+* MOCSCAT: Step-Characteristics Scheme with Tabulated Exponentials
+* MOCDDFT: Diamond-Differencing Scheme
+* MOCSCET: Step-Characteristics Scheme with Exact Exponentials
+ CALL SUBSCH(LINE,NR,NS,MT,NRSEG,MATALB,SEGLEN,SIGANG(-6),EXPT,
+ 1 EXP2,NE,YG)
+*----
+* Summation along the track in both directions
+*----
+ DO 4226 IL= 1, LINE
+ JL= LINE + 1 - IL
+ DO 4225 IE=1,NE
+ FLM(IE,IL) = (FL0(IE) - FLM(IE,IL)) * EXP2(IE,IL)
+ CYM(IE,IL) = CY0(IE) * EXP2(IE,IL)
+ FLP(IE,JL) = (FL1(IE) - FLP(IE,JL)) * EXP2(IE,JL)
+ CYP(IE,JL) = CY1(IE) * EXP2(IE,JL)
+ FL0(IE) = FL0(IE) - FLM(IE,IL)
+ CY0(IE) = CY0(IE) - CYM(IE,IL)
+ FL1(IE) = FL1(IE) - FLP(IE,JL)
+ CY1(IE) = CY1(IE) - CYP(IE,JL)
+ 4225 CONTINUE
+ 4226 CONTINUE
+ DO 4230 IE=1,NE
+ FL0(IE)= FL0(IE)/(ONE-CY0(IE))
+ FL1(IE)= FL1(IE)/(ONE-CY1(IE))
+ 4230 CONTINUE
+ DO 4240 IL= 1, LINE
+ NOIL = NRSEG(IL)
+ IF( NOIL.GT.0 )THEN
+ IND=KEYFLX(NOIL)
+ INDC=NUN/2+IND
+ DO 4241 IE=1,NE
+ FLUX(IND)= FLUX(IND)
+ > + ((FL0(IE)*CYM(IE,IL)+FLM(IE,IL))
+ > +(FL1(IE)*CYP(IE,IL)+FLP(IE,IL)))
+* CALCULATE XI, YI OR ZI FOR TIBERE
+ IF(IDIR.GT.0) THEN
+ FLUX(INDC)= FLUX(INDC)
+ > + ((FL0(IE)*CYM(IE,IL)+FLM(IE,IL))
+ > +(FL1(IE)*CYP(IE,IL)+FLP(IE,IL)))
+ > *OMG2(IE,IDIR)
+ ENDIF
+ 4241 CONTINUE
+ ENDIF
+ 4240 CONTINUE
+*
+ RETURN
+ END