SUBROUTINE InitConductionTransferFunctions
! SUBROUTINE INFORMATION:
! AUTHOR Russ Taylor
! DATE WRITTEN June 1990
! MODIFIED July 1994, LKL, cosmetic and improve execution time
! Dec 1995, Apr 1996, RKS, cosmetic and clean-up changes, changes to allow proper
! handling of resistive layers
! June 2000, RKS, addition of QTFs (both 1- and 2-D solutions for constructions
! with embedded/internal heat sources/sinks)
! July 2010-August 2011, RKS, R-value only layer enhancement
! RE-ENGINEERED June 1996, February 1997, August-October 1997, RKS; Nov 1999, LKL
! PURPOSE OF THIS SUBROUTINE:
! This subroutine serves as the main drive for the
! calculation of Conduction Transfer Functions (CTFs)
! using the state space method.
! METHODOLOGY EMPLOYED:
! The basic steps of this routine (which may be a little difficult
! to decipher until another major revision is done) are:
! 1. Determine if enough material info has been entered
! 2. Determine whether construct is (a) all resistive,
! (b) the reverse of a previously calculated construct, or
! (c) neither (a) nor (b), i.e. a layer for which CTFs must
! be calculated.
! 3. If the answer to 2 is (a), calculate the overall resistance
! and use this as the CTF (steady state conduction).
! 4. If the answer to 2 is (b), transfer the CTFs for the reverse
! construction to the CTF arrays for this construct (reversing
! the inside and outside terms).
! 5. If the answer to 2 is (c), calculate the CTFs using the state
! space method described below.
! The state space method of calculating CTFs involves
! applying a finite difference grid to a multilayered
! building element and performing linear algebra on the
! resulting system of equations (in matrix form).
! CTFs must be calculated for non-reversed layers which
! have an appreciable thermal mass. A conversion from
! SI units to English units is made due to concerns
! about round off problems noted in earlier version of
! this subroutine.
! REFERENCES:
! Seem, J.E. "Modeling of Heat Transfer in Buildings",
! Department of Mechanical Engineering, University of
! Wisconsin-Madison, 1987.
! Strand, R.K. "Testing Design Description for the CTF
! Calculation Code in BEST", BSO internal document,
! May/June 1996.
! Strand, R.K. "Heat Source Transfer Functions and Their
! Applicatoin to Low Temperature Radiant Heating System",
! Ph.D. Dissertation, Department of Mechanical and
! Industrial Engineering, University of Illinois at
! Urbana-Champaign, 1995.
! USE STATEMENTS:
USE DataConversions
USE General, ONLY: RoundSigDigits
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
! na
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: PhysPropLimit = 1.0d-6 ! Physical properties limit.
! This is more or less the traditional value from BLAST.
REAL(r64), PARAMETER :: RValueLowLimit = 1.0d-3 ! Physical properties limit for R-value only layers
! This value was based on trial and error related to CR 7791 where a
! user had entered a "no insulation" layer with an R-value of 1.0E-05.
! Some trial and error established this as a potential value though
! there is no guarantee that this is a good value.
INTEGER, PARAMETER :: MinNodes = 6 ! Minimum number of state space nodes
! per layer. This value was chosen based on experience with IBLAST.
REAL(r64), PARAMETER :: MaxAllowedCTFSumError = 0.01d0 ! Allow a 1 percent
! difference between the CTF series summations. If the difference is
! greater than this, then the coefficients will not yield a valid steady
! state solution.
REAL(r64), PARAMETER :: MaxAllowedTimeStep = 4.0d0 ! Sets the maximum allowed time step
! for CTF calculations to be 4 hours. This is done in response to some
! rare situations where odd or faulty input will cause the routine to
! go off and get some huge time step (in excess of 20 hours). This value
! is a compromise that does not really solve any input problems. One run
! indicated that 2 meters of concrete will result in a time step of slightly
! more than 3 hours. So, 4 hours was arbitrarily picked as a ceiling for
! time steps so that an error message can be produced to warn the user
! that something isn't right. Note that the 4 hour limit does not guarantee
! that problems won't exist and it does not necessarily avoid any problems
! that interpolated temperature histories might cause.
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
INTEGER, &
DIMENSION(MaxLayersInConstruct) :: AdjacentResLayerNum ! Layers that are adjacent to each other which are resistive
! only can and should be combined
INTEGER :: AdjLayer ! Loop counter for adjacent resistance-only layers
REAL(r64) :: amatx ! Intermediate calculation variable
REAL(r64) :: amatxx ! Intermediate calculation variable
REAL(r64) :: amaty ! Intermediate calculation variable
REAL(r64) :: BiggestSum ! Largest CTF series summation (maximum of SumXi, SumYi, and SumZi)
REAL(r64) :: cap ! Thermal capacitance of a node (intermediate calculation)
REAL(r64) :: capavg ! Thermal capacitance of a node (average value for a node at an interface)
REAL(r64) :: cnd ! Total thermal conductance (1/Rtot) of the bldg element
INTEGER :: Constr ! Loop counter
INTEGER :: ConstrNum ! Loop counter (construct number)
REAL(r64), &
DIMENSION(MaxLayersInConstruct) :: cp ! Specific heat of a material layer
LOGICAL :: CTFConvrg ! Set after CTFs are calculated, based on whether there are too
! many CTF terms
INTEGER :: CurrentLayer ! Pointer to material number in Material derived type (current layer)
REAL(r64), &
DIMENSION(MaxLayersInConstruct) :: dl ! Thickness of a material layer
REAL(r64) :: dtn ! Intermediate calculation of the time step
REAL(r64), &
DIMENSION(MaxLayersInConstruct) :: dx ! Distance between nodes in a particular material layer
REAL(r64) :: dxn ! Intermediate calculation of nodal spacing
REAL(r64) :: dxtmp ! Intermediate calculation variable ( = 1/dx/cap)
REAL(r64) :: dyn ! Nodal spacing in the direction perpendicular to the main direction
! of heat transfer (only valid for a 2-D solution)
LOGICAL :: ErrorsFound=.false. !Flag for input error condition
INTEGER :: HistTerm ! Loop counter
INTEGER :: ipts1 ! Intermediate calculation for number of nodes per layer
INTEGER :: ir ! Loop control for constructing Identity Matrix
INTEGER :: Layer ! Loop counter
INTEGER :: Layer1 ! Loop counter
INTEGER :: LayersInConstruct ! Array containing the number of layers for each construct
! Different from TotLayers because shades are not include in local var
REAL(r64), &
DIMENSION(MaxLayersInConstruct) :: lr ! R value of a material layer
INTEGER :: Node ! Loop counter
INTEGER :: Node2 ! Node number (modification of Node and NodeInRow)
INTEGER :: NodeInLayer ! Loop counter
INTEGER :: NodeInRow ! Loop counter
INTEGER, &
DIMENSION(MaxLayersInConstruct) :: Nodes ! Array containing the number of nodes per layer
INTEGER :: NumResLayers ! Number of resistive layers in the construction
! property allowable, traditional value from BLAST
INTEGER :: NumAdjResLayers ! Number of resistive layers that are adjacent
INTEGER :: OppositeLayer ! Used for comparing constructions (to see if one is the reverse of another)
LOGICAL, &
DIMENSION(MaxLayersInConstruct) :: ResLayer ! Set true if the layer must be handled as a resistive
LOGICAL :: RevConst ! Set true if one construct is the reverse of another (CTFs already
! available)
REAL(r64), &
DIMENSION(MaxLayersInConstruct) :: rho ! Density of a material layer
REAL(r64), &
DIMENSION(MaxLayersInConstruct) :: rk ! Thermal conductivity of a material layer
REAL(r64) :: rs ! Total thermal resistance of the building element
REAL(r64) :: SumXi ! Summation of all of the Xi terms (inside CTFs) for a construction
REAL(r64) :: SumYi ! Summation of all of the Xi terms (cross CTFs) for a construction
REAL(r64) :: SumZi ! Summation of all of the Xi terms (outside CTFs) for a construction
LOGICAL :: DoCTFErrorReport
REAL(r64) :: Alpha ! thermal diffusivity in m2/s, for local check of properties
REAL(r64) :: DeltaTimestep ! zone timestep in seconds, for local check of properties
REAL(r64) :: ThicknessThreshold ! min thickness consistent with other thermal properties, for local check
! FLOW:
! Subroutine initializations
TinyLimit=rTinyValue
DoCTFErrorReport=.false.
DO ConstrNum = 1, TotConstructs ! Begin construction loop ...
Construct(ConstrNum)%CTFCross = 0.0D0
Construct(ConstrNum)%CTFFlux = 0.0D0
Construct(ConstrNum)%CTFInside = 0.0D0
Construct(ConstrNum)%CTFOutside = 0.0D0
Construct(ConstrNum)%CTFSourceIn = 0.0D0
Construct(ConstrNum)%CTFSourceOut = 0.0D0
Construct(ConstrNum)%CTFTimeStep = 0.0D0
Construct(ConstrNum)%CTFTSourceOut = 0.0D0
Construct(ConstrNum)%CTFTSourceIn = 0.0D0
Construct(ConstrNum)%CTFTSourceQ = 0.0D0
Construct(ConstrNum)%CTFTUserOut = 0.0D0
Construct(ConstrNum)%CTFTUserIn = 0.0D0
Construct(ConstrNum)%CTFTUserSource = 0.0D0
Construct(ConstrNum)%NumHistories = 0
Construct(ConstrNum)%NumCTFTerms = 0
Construct(ConstrNum)%UValue = 0.0d0
AdjacentResLayerNum = 0 ! Zero this out for each construct
IF(Construct(ConstrNum)%TypeIsWindow) CYCLE
! Initialize construct parameters
Construct(ConstrNum)%CTFTimeStep = TimeStepZone
rs = 0.0D0
LayersInConstruct = 0
NumResLayers = 0
ResLayer = .FALSE.
DO Layer = 1, Construct(ConstrNum)%TotLayers ! Begin layer loop ...
! Loop through all of the layers in the current construct. The purpose
! of this loop is to define the thermal properties necessary to
! calculate the CTFs.
CurrentLayer = Construct(ConstrNum)%LayerPoint(Layer)
LayersInConstruct = LayersInConstruct + 1
! Obtain thermal properties from the Material derived type
dl(Layer) = Material(CurrentLayer)%Thickness
rk(Layer) = Material(CurrentLayer)%Conductivity
rho(Layer) = Material(CurrentLayer)%Density
cp(Layer) = Material(CurrentLayer)%SpecHeat ! Must convert
! from kJ/kg-K to J/kg-k due to rk units
IF (Construct(ConstrNum)%SourceSinkPresent .AND. .Not. Material(CurrentLayer)%WarnedForHighDiffusivity) THEN
! check for materials that are too conductive or thin
IF ((rho(Layer) * cp(Layer)) > 0.d0) THEN
Alpha = rk(Layer) / (rho(Layer) * cp(Layer))
IF (Alpha > HighDiffusivityThreshold ) THEN
DeltaTimestep = TimeStepZone * SecInHour
ThicknessThreshold = SQRT(Alpha * DeltaTimestep * 3.d0)
IF (Material(CurrentLayer)%Thickness < ThicknessThreshold) THEN
CALL ShowSevereError('InitConductionTransferFunctions: Found Material that is too thin and/or too ' &
//'highly conductive, material name = ' // TRIM(Material(CurrentLayer)%Name))
CALL ShowContinueError('High conductivity Material layers are not well supported for internal source ' &
//'constructions, material conductivity = ' &
//TRIM(RoundSigDigits(Material(CurrentLayer)%Conductivity, 3)) //' [W/m-K]')
CALL ShowContinueError('Material thermal diffusivity = ' //TRIM(RoundSigDigits(Alpha, 3)) //' [m2/s]')
CALL ShowContinueError('Material with this thermal diffusivity should have thickness > ' &
// TRIM(RoundSigDigits(ThicknessThreshold , 5)) // ' [m]')
IF (Material(CurrentLayer)%Thickness < ThinMaterialLayerThreshold) THEN
CALL ShowContinueError('Material may be too thin to be modeled well, thickness = ' &
// TRIM(RoundSigDigits(Material(currentLayer)%Thickness, 5 ))//' [m]')
CALL ShowContinueError('Material with this thermal diffusivity should have thickness > ' &
// TRIM(RoundSigDigits(ThinMaterialLayerThreshold , 5)) // ' [m]')
ENDIF
Material(CurrentLayer)%WarnedForHighDiffusivity = .TRUE.
ENDIF
ENDIF
ENDIF
ENDIF
IF (rk(Layer) <= PhysPropLimit) THEN ! Thermal conductivity too small,
! thus this must be handled as a resistive layer
ResLayer(Layer) = .TRUE.
ELSE
lr(Layer) = dl(Layer)/rk(Layer)
IF ((dl(Layer)*SQRT(rho(Layer)*cp(Layer)/rk(Layer))) &
< PhysPropLimit) THEN
! If the time constant is smaller than some reasonable
! limit, this should be treated as a resistive layer.
ResLayer(Layer) = .TRUE.
ELSE ! Layer has significant thermal mass (non-resistive)
ResLayer(Layer) = .FALSE.
END IF
END IF
! If not a resistive layer, nothing further is required
! for this layer.
IF (ResLayer(Layer)) THEN ! Resistive layer-check for R-value, etc.
NumResLayers = NumResLayers + 1 ! Increment number of resistive layers
lr(Layer) = Material(CurrentLayer)%Resistance ! User defined thermal resistivity
IF (lr(Layer) < RValueLowLimit) THEN ! User didn't define enough
! parameters to calculate CTFs for a building element
! containing this layer.
CALL ShowSevereError('InitConductionTransferFunctions: Material='//TRIM(Material(CurrentLayer)%Name)// &
'R Value below lowest allowed value')
CALL ShowContinueError('Lowest allowed value=['//trim(RoundSigDigits(RValueLowLimit,3))// &
'], Material R Value=['//trim(RoundSigDigits(lr(Layer),3))//'].')
ErrorsFound=.true.
ELSE ! A valid user defined R-value is available.
! If this is either the first or last layer in the construction,
! then assign other properties based on air at 1 atm, 300K.
! Reference for air properties: Incropera and DeWitt,
! Introduction to Heat Transfer, Appendix A, Table A.4,
! John Wiley & Sons, New York, 1985.
! If this is not the first or last layer in the construction,
! then use the "exact" approach to model a massless layer
! based on the node equations for the state space method.
IF ( (Layer == 1) .OR. (Layer == Construct(ConstrNum)%TotLayers) .OR. &
(.NOT. Material(Construct(ConstrNum)%LayerPoint(Layer))%ROnly) ) THEN
cp(Layer) = 1.007d0
rho(Layer) = 1.1614d0
rk(Layer) = 0.0263d0
dl(Layer) = rk(Layer)*lr(Layer)
ELSE
cp(Layer) = 0.0d0
rho(Layer) = 0.0d0
rk(Layer) = 1.0d0
dl(Layer) = lr(Layer)
END IF
END IF
END IF ! ... end of resistive layer determination IF-THEN block.
END DO ! ... end of layer loop.
! If errors have been found, just cycle
IF (ErrorsFound) CYCLE
! Combine any adjacent resistive-only (no mass) layers together
! to avoid a divide by zero error in the CTF calculations below.
! Since the inner and outer layers cannot be resistive layers
! (inner and outer layer still converted to equivalent air layer)
! there can only be resistive layers adjacent to one another if
! there are more than three total layers and more than one
! resistive layer.
IF ((LayersInConstruct > 3).AND.(NumResLayers > 1)) THEN
NumAdjResLayers = 0
DO Layer = 2, LayersInConstruct-2
IF ( (ResLayer(Layer)) .AND. (ResLayer(Layer+1)) ) THEN
NumAdjResLayers = NumAdjResLayers + 1
! There is method to the next assignment statement. As the layers get shifted, the layer
! numbers will also shift. Thus, we have to also shift which layer we are dealing with.
AdjacentResLayerNum(NumAdjResLayers) = Layer + 1 - NumAdjResLayers
END IF
END DO
DO AdjLayer = 1, NumAdjResLayers
Layer = AdjacentResLayerNum(AdjLayer)
! Double check to make sure we are in the right place...
IF ( (ResLayer(Layer)) .AND. (ResLayer(Layer+1)) ) THEN
! Shift layers forward after combining two adjacent layers. Then
! restart the do loop.
cp(Layer) = 0.0d0
rho(Layer) = 0.0d0
rk(Layer) = 1.0d0
lr(Layer) = lr(Layer) + lr(Layer+1)
dl(Layer) = lr(Layer)
NumResLayers = NumResLayers - 1 ! Combining layers so decrease number of resistive layers
DO Layer1 = Layer+1, LayersInConstruct - 1
lr(Layer1) = lr(Layer1+1)
dl(Layer1) = dl(Layer1+1)
rk(Layer1) = rk(Layer1+1)
rho(Layer1) = rho(Layer1+1)
cp(Layer1) = cp(Layer1+1)
ResLayer(Layer1) = ResLayer(Layer1+1)
END DO
! Then zero out the layer that got shifted forward
cp(LayersInConstruct) = 0.0d0
rho(LayersInConstruct) = 0.0d0
rk(LayersInConstruct) = 0.0d0
lr(LayersInConstruct) = 0.0d0
dl(LayersInConstruct) = 0.0d0
! Now reduce the number of layers in construct since merger is complete
LayersInConstruct = LayersInConstruct - 1
! Also adjust layers with source/sinks if two layers are merged
IF (Construct(ConstrNum)%SourceSinkPresent) THEN
Construct(ConstrNum)%SourceAfterLayer = Construct(ConstrNum)%SourceAfterLayer - 1
Construct(ConstrNum)%TempAfterLayer = Construct(ConstrNum)%TempAfterLayer - 1
END IF
ELSE ! These are not adjacent layers and there is a logic flaw here (should not happen)
CALL ShowFatalError('Combining resistance layers failed for '//TRIM(Construct(ConstrNum)%Name))
CALL ShowContinueError('This should never happen. Contact EnergyPlus Support for further assistance.')
END IF
END DO
END IF
! Convert SI units to English. In theory, conversion to English
! units is not necessary; however, Russ Taylor noted that some
! numerical problems when SI units were used and decided to continue
! calculating CTFs in English units.
DO Layer = 1,LayersInConstruct ! Begin units conversion loop ...
lr(Layer) = lr(Layer)*CFU
dl(Layer) = dl(Layer)/CFL
rk(Layer) = rk(Layer)/CFK
rho(Layer) = rho(Layer)/CFD
cp(Layer) = cp(Layer)/(CFC*1000.0d0)
END DO ! ... end of layer loop for units conversion.
IF (Construct(ConstrNum)%SolutionDimensions == 1) THEN
dyn = 0.0d0
ELSE
dyn = (Construct(ConstrNum)%ThicknessPerpend/CFL)/REAL(NumOfPerpendNodes-1,r64)
END IF
! Compute total construct conductivity and resistivity.
DO Layer = 1,LayersInConstruct
rs = rs + lr(Layer) ! Resistances in series sum algebraically
END DO
cnd = 1.0D0/rs ! Conductivity is the inverse of resistivity
IF (LayersInConstruct > NumResLayers) THEN
! One or more are not simple resistive layers so CTFs will have to be
! calculated unless this is a reverse of a previously defined
! construction.
! Check for reversed construction of interzone surfaces by checking
! previous constructions for same number of layers as first indicator.
RevConst = .FALSE.
DO Constr = 1, ConstrNum-1 ! Constructions loop (check reversed) ...
! If a source or sink is present in this construction, do not allow any
! checks for reversed constructions, i.e., always force EnergyPlus to
! calculate CTF/QTFs. So, don't even check for reversed constructions.
IF (Construct(ConstrNum)%SourceSinkPresent) EXIT ! Constr DO loop
IF (Construct(ConstrNum)%TotLayers == & ! Same number of layers--now
Construct(Constr)%TotLayers) THEN ! check for reversed construct.
RevConst = .TRUE.
DO Layer = 1, Construct(ConstrNum)%TotLayers ! Begin layers loop ...
! RevConst is set to FALSE anytime a mismatch in materials is found.
! This will exit this DO immediately and go on to the next construct
! (if any remain).
OppositeLayer = Construct(ConstrNum)%TotLayers-Layer+1
IF (Construct(ConstrNum)%LayerPoint(Layer) /= &
Construct(Constr)%LayerPoint(OppositeLayer)) THEN
RevConst = .FALSE.
EXIT ! Layer DO loop
END IF
END DO ! ... end of layers loop.
IF (RevConst) THEN ! Curent construction is a reverse of
! construction Constr. Thus, CTFs do not need to be re-
! calculated. Copy CTF info for construction Constr to
! construction ConstrNum.
Construct(ConstrNum)%CTFTimeStep = Construct(Constr)%CTFTimeStep
Construct(ConstrNum)%NumHistories = Construct(Constr)%NumHistories
Construct(ConstrNum)%NumCTFTerms = Construct(Constr)%NumCTFTerms
! Transfer the temperature and flux history terms to CTF arrays.
! Loop through the number of CTF history terms ...
DO HistTerm = 0, Construct(ConstrNum)%NumCTFTerms
Construct(ConstrNum)%CTFInside(HistTerm) = Construct(Constr)%CTFOutside(HistTerm)
Construct(ConstrNum)%CTFCross(HistTerm) = Construct(Constr)%CTFCross(HistTerm)
Construct(ConstrNum)%CTFOutside(HistTerm) = Construct(Constr)%CTFInside(HistTerm)
IF (HistTerm /= 0) &
Construct(ConstrNum)%CTFFlux(HistTerm) = Construct(Constr)%CTFFlux(HistTerm)
END DO ! ... end of CTF history terms loop.
EXIT ! Constr DO loop
END IF ! ... end of reversed construction found block
END IF ! ... end of reversed construct (same number of layers) block.
END DO ! ... end of construct loop (check reversed--Constr)
IF (.NOT.RevConst) THEN ! Calculate CTFs (non-reversed constr)
! Estimate number of nodes each layer of the construct will require
! and calculate the nodal spacing from that
DO Layer = 1, LayersInConstruct ! Begin loop thru layers ...
! The calculation of dxn used here is based on a standard stability
! criteria for explicit finite difference solutions. This criteria
! was chosen not because it is viewed to be correct, but rather for
! lack of any better criteria at this time. The use of a Fourier
! number based criteria such as this is probably physically correct,
! though the coefficient (2.0) may not be.
! If this is a "resistive" layer, only need a single node
IF ( (ResLayer(Layer)) .AND. (Layer>1) .AND. (Layer<LayersInConstruct) ) THEN
Nodes(Layer) = 1
dx(Layer) = dl(Layer)
ELSE
dxn = sqrt(2.0d0*(rk(Layer)/rho(Layer)/cp(Layer)) &
*Construct(ConstrNum)%CTFTimeStep)
ipts1 = int(dl(Layer)/dxn) ! number of nodes=thickness/spacing
! Limit the upper and lower bounds of the number of
! nodes to MaxCTFTerms and MinNodes respectively.
IF (ipts1 > MaxCTFTerms) THEN ! Too many nodes
Nodes(Layer) = MaxCTFTerms
ELSE IF (ipts1 < MinNodes) THEN ! Too few nodes
Nodes(Layer) = MinNodes
ELSE ! Calculated number of nodes ok
Nodes(Layer) = ipts1
END IF
IF (Construct(ConstrNum)%SolutionDimensions > 1) THEN
IF (ipts1 > MaxCTFTerms/2) ipts1 = MaxCTFTerms/2
END IF
dx(Layer) = dl(Layer)/REAL(Nodes(Layer),r64) ! calc node spacing
END IF
END DO ! . .. end of layers in construction loop (calculating #nodes per layer)
! Determine the total number of nodes (rcmax)
rcmax = 0
DO Layer = 1, LayersInConstruct
rcmax = rcmax + Nodes(Layer)
END DO
! Nodes are placed throughout layers and at the interface between
! layers. As a result, the end layers share a node with the adjacent
! layer-leaving one less node total for all layers.
rcmax = rcmax-1
IF (Construct(ConstrNum)%SolutionDimensions > 1) rcmax = rcmax * NumOfPerpendNodes
! This section no longer needed as rcmax/number of total nodes is allowed to float.
! If reinstated, this node reduction section would have to be modified to account for
! the possibility that a 2-D solution is potentially being performed.
! Check to see if the maximum number of nodes for the construct has
! been exceeded. Reduce the nodes per layer if necessary, but only
! if the number of nodes in a particular layer is greater than the
! minimum node limit.
! DO WHILE (rcmax > MaxTotNodes) ! Begin total node reduction loop ...
! rcmax = 0
! DO Layer = 1, LayersInConstruct ! Begin layer node reduction ...
! ! If more nodes than the minimum limit for a layer, reduce the
! ! number of nodes.
! IF (Nodes(Layer) > MinNodes) THEN
! Nodes(Layer) = Nodes(Layer)-1
! dx(Layer) = dl(Layer)/dble(float(Nodes(Layer))) ! Recalc node spacing
! END IF
! rcmax = rcmax + Nodes(Layer) ! Recalculate total number of nodes
! END DO ! ... end of layer loop for node reduction.
! rcmax = rcmax-1 ! See note above on counting rcmax
! END DO ! ... end of total node reduction loop.
! For constructions that have sources or sinks present, determine which
! node the source/sink is applied at and also where the temperature
! calculation has been requested.
NodeSource = 0
NodeUserTemp = 0
IF (Construct(ConstrNum)%SourceSinkPresent) THEN
DO Layer = 1, Construct(ConstrNum)%SourceAfterLayer
NodeSource = NodeSource + Nodes(Layer)
END DO
IF ( (NodeSource > 0) .AND. (Construct(ConstrNum)%SolutionDimensions > 1) ) &
NodeSource = ((NodeSource-1)*NumOfPerpendNodes) + 1
DO Layer = 1, Construct(ConstrNum)%TempAfterLayer
NodeUserTemp = NodeUserTemp + Nodes(Layer)
END DO
IF ( (NodeUserTemp > 0) .AND. (Construct(ConstrNum)%SolutionDimensions > 1) ) &
NodeUserTemp = ((NodeUserTemp-1)*NumOfPerpendNodes) + 1
END IF
! "Adjust time step to ensure stability." If the time step is too
! small, it will result in too many history terms which can lead to
! solution instability. The method used here to determine whether or
! not the time step will produce a stable solution is based on a pure
! Fourier number calculation (Fo = 1) and has not proven to be
! completely effective. If too many history terms are calculated,
! the time step is adjusted and the CTFs end up being recalculated
! (see later code in this routine).
dtn = 0.0D0
Construct(ConstrNum)%CTFTimeStep = 0.0D0
DO Layer = 1, LayersInConstruct
IF (Nodes(Layer) >= MaxCTFTerms) THEN
IF (Construct(ConstrNum)%SolutionDimensions == 1) THEN
dtn = rho(Layer)*cp(Layer)*dx(Layer)**2/rk(Layer)
ELSE ! 2-D solution requested-->this changes length parameter in Fourier number calculation
dtn = rho(Layer)*cp(Layer)*((dx(Layer)**2)+(dyn**2))/rk(Layer)
END IF
IF (dtn > Construct(ConstrNum)%CTFTimeStep) &
Construct(ConstrNum)%CTFTimeStep = dtn
END IF
END DO
! If the user defined time step is significantly different than the
! calculated time step for this construct, then CTFTimeStep must be
! revised.
IF (ABS((TimeStepZone-Construct(ConstrNum)%CTFTimeStep)/TimeStepZone) > 0.1d0) THEN
IF (Construct(ConstrNum)%CTFTimeStep > TimeStepZone) THEN
! CTFTimeStep larger than TimeStepZone: Make sure TimeStepZone
! divides evenly into CTFTimeStep
Construct(ConstrNum)%NumHistories = &
INT((Construct(ConstrNum)%CTFTimeStep/TimeStepZone)+0.5D0)
Construct(ConstrNum)%CTFTimeStep = &
TimeStepZone*REAL(Construct(ConstrNum)%NumHistories,r64)
ELSE
! CTFTimeStep smaller than TimeStepZone: Set to TimeStepZone
Construct(ConstrNum)%CTFTimeStep = TimeStepZone
Construct(ConstrNum)%NumHistories = 1
END IF
END IF
! Calculate the CTFs using the state space method
! outlined in Seem's dissertation. The main matrices
! AMat, BMat, CMat, and DMat must be derived from
! applying a finite difference network to the layers of
! each bldg element.
! This section must continue looping until the CTFs
! calculated here will produce a stable solution (less
! history terms than MaxCTFTerms).
! This first subsection calculates the elements of AMat
! which characterizes the heat transfer inside the
! building element.
CTFConvrg = .FALSE. ! Initialize loop control logical
ALLOCATE(AExp(rcmax,rcmax))
AExp=0.0D0
ALLOCATE(AMat(rcmax,rcmax))
AMat=0.0D0
ALLOCATE(AInv(rcmax,rcmax))
AInv=0.0D0
ALLOCATE(IdenMatrix(rcmax,rcmax))
IdenMatrix=0.0D0
DO ir=1,rcmax
IdenMatrix(ir,ir)=1.0D0
ENDDO
ALLOCATE(e(rcmax))
e=0.0D0
ALLOCATE(Gamma1(rcmax,3))
Gamma1=0.0D0
ALLOCATE(Gamma2(rcmax,3))
Gamma2=0.0D0
ALLOCATE(s(rcmax,4,3))
s=0.0D0
DO WHILE (.NOT.CTFConvrg) ! Begin CTF calculation loop ...
BMat(3) = 0.0d0
IF (Construct(ConstrNum)%SolutionDimensions == 1) THEN
! Set up intermediate calculations for the first layer.
cap = rho(1)*cp(1)*dx(1)
cap = 1.5d0*cap ! For the first node, account for the fact that the
! half-node at the surface results in a "loss" of some
! thermal mass. Therefore, for simplicity, include it
! at this node. Same thing done at the last node...
dxtmp = 1.0D0/dx(1)/cap
AMat(1,1) = -2.0D0*rk(1)*dxtmp ! Assign the matrix values for the
AMat(1,2) = rk(1)*dxtmp ! first node.
BMat(1) = rk(1)*dxtmp ! Assign non-zero value of BMat.
Layer = 1 ! Initialize the "layer" counter
NodeInLayer = 2 ! Initialize the node (in a layer) counter (already
! on the second node for the first layer
DO Node = 2,rcmax-1 ! Begin nodes loop (includes all nodes except the
! first/last which have special equations) ...
IF ((NodeInLayer == Nodes(Layer)) .AND. &
(LayersInConstruct /= 1)) THEN ! For a node at
! the interface between two adjacent layers, the
! capacitance of the node must be calculated from the 2
! halves which may be made up of 2 different materials.
cap = ( rho(Layer)*cp(Layer)*dx(Layer) &
+rho(Layer+1)*cp(Layer+1)*dx(Layer+1) ) * 0.5D0
AMat(Node,Node-1) = rk(Layer)/dx(Layer)/cap ! Assign matrix
AMat(Node,Node) = -1.0D0 * ( rk(Layer)/dx(Layer)+ & ! values for
rk(Layer+1)/dx(Layer+1) ) / cap ! the current
AMat(Node,Node+1) = rk(Layer+1)/dx(Layer+1)/cap ! node.
NodeInLayer = 0 ! At an interface, reset nodes in layer counter
Layer = Layer+1 ! Also increment the layer counter
ELSE ! Standard node within any layer
cap = rho(Layer)*cp(Layer)*dx(Layer) ! Intermediate
dxtmp = 1.0D0/dx(Layer)/cap ! calculations.
AMat(Node,Node-1) = rk(Layer)*dxtmp ! Assign matrix
AMat(Node,Node) = -2.0D0*rk(Layer)*dxtmp ! values for the
AMat(Node,Node+1) = rk(Layer)*dxtmp ! current node.
END IF
NodeInLayer = NodeInLayer+1 ! Increment nodes in layer counter
IF (Node == NodeSource) BMat(3) = 1.0d0/cap
END DO ! ... end of nodes loop.
! Intermediate calculations for the last node.
cap = rho(LayersInConstruct)*cp(LayersInConstruct)* dx(LayersInConstruct)
cap = 1.5d0*cap ! For the last node, account for the fact that the
! half-node at the surface results in a "loss" of some
! thermal mass. Therefore, for simplicity, include it
! at this node. Same thing done at the first node...
dxtmp = 1.0D0/dx(LayersInConstruct)/cap
AMat(rcmax,rcmax) = -2.0D0*rk(LayersInConstruct)*dxtmp ! Assign matrix
AMat(rcmax,rcmax-1) = rk(LayersInConstruct)*dxtmp ! values for the
BMat(2) = rk(LayersInConstruct)*dxtmp ! last node.
CMat(1) = -rk(1)/dx(1) ! Compute the necessary elements
CMat(2) = rk(LayersInConstruct)/dx(LayersInConstruct) ! of all other
DMat(1) = rk(1)/dx(1) ! matrices for the state
DMat(2) = -rk(LayersInConstruct)/dx(LayersInConstruct) ! space method
ELSE ! 2-D solution requested (assign matrices appropriately)
! As with the 1-D solution, we are accounting for the thermal mass
! of the half-node at the surface by adding it to the first row
! of interior nodes at both sides of the construction. This is not
! exact, but it does take all of the thermal mass into account.
amatx = rk(1)/(1.5d0*rho(1)*cp(1)*dx(1)*dx(1))
amaty = rk(1)/(1.5d0*rho(1)*cp(1)*dyn*dyn)
! FIRST ROW OF NODES: This first row within the first material layer
! is special in that it is exposed to a boundary condition. Thus,
! the equations are slightly different.
! Note also that the first and last nodes in a row are slightly
! different from the rest since they are on an adiabatic plane in
! the direction perpendicular to the main direction of heat transfer.
AMat(1,1) =-2.0d0*(amatx+amaty)
AMat(1,2) = 2.0d0*amaty
AMat(1,NumOfPerpendNodes+1) = amatx
DO Node = 2, NumOfPerpendNodes-1
AMat(Node,Node-1) = amaty
AMat(Node,Node) =-2.0d0*(amatx+amaty)
AMat(Node,Node+1) = amaty
AMat(Node,Node+NumOfPerpendNodes) = amatx
END DO
AMat(NumOfPerpendNodes,NumOfPerpendNodes) =-2.0d0*(amatx+amaty)
AMat(NumOfPerpendNodes,NumOfPerpendNodes-1) = 2.0d0*amaty
AMat(NumOfPerpendNodes,NumOfPerpendNodes+NumOfPerpendNodes) = amatx
BMat(1) = amatx
Layer = 1
NodeInLayer = 2
amatx = rk(1)/(rho(1)*cp(1)*dx(1)*dx(1)) ! Reset these to the normal capacitance
amaty = rk(1)/(rho(1)*cp(1)*dyn*dyn) ! Reset these to the normal capacitance
DO Node = NumOfPerpendNodes+1, rcmax+1-2*NumOfPerpendNodes, NumOfPerpendNodes
! INTERNAL ROWS OF NODES: This is the majority of nodes which are all within
! a solid layer and not exposed to a boundary condition.
IF ((LayersInConstruct == 1).OR.(NodeInLayer /= Nodes(Layer))) THEN
! Single material row: This row of nodes are all contained within a material
! and thus there is no special considerations necessary.
IF (NodeInLayer == 1) THEN
! These intermediate variables only need to be reassigned when a new layer is started.
! When this is simply another row of the same material, these have already been assigned correctly.
amatx = rk(Layer)/(rho(Layer)*cp(Layer)*dx(Layer)*dx(Layer))
amaty = rk(Layer)/(rho(Layer)*cp(Layer)*dyn*dyn)
END IF
! Note that the first and last layers in a row are slightly different
! from the rest since they are on an adiabatic plane in the direction
! perpendicular to the main direction of heat transfer.
AMat(Node,Node) =-2.0d0*(amatx+amaty)
AMat(Node,Node+1) = 2.0d0*amaty
AMat(Node,Node-NumOfPerpendNodes) = amatx
AMat(Node,Node+NumOfPerpendNodes) = amatx
DO NodeInRow = 2, NumOfPerpendNodes-1
Node2 = Node + NodeInRow - 1
AMat(Node2,Node2-1) = amaty
AMat(Node2,Node2) =-2.0d0*(amatx+amaty)
AMat(Node2,Node2+1) = amaty
AMat(Node2,Node2-NumOfPerpendNodes) = amatx
AMat(Node2,Node2+NumOfPerpendNodes) = amatx
END DO
Node2 = Node - 1 + NumOfPerpendNodes
AMat(Node2,Node2) =-2.0d0*(amatx+amaty)
AMat(Node2,Node2-1) = 2.0d0*amaty
AMat(Node2,Node2-NumOfPerpendNodes) = amatx
AMat(Node2,Node2+NumOfPerpendNodes) = amatx
ELSE ! Row at a two-layer interface (half of node consists of one layer's materials
! and the other half consist of the next layer's materials)
capavg = 0.5d0*(rho(Layer)*cp(Layer)*dx(Layer)+rho(Layer+1)*cp(Layer+1)*dx(Layer+1))
amatx = rk(Layer)/(capavg*dx(Layer))
amatxx = rk(Layer+1)/(capavg*dx(Layer+1))
amaty = (rk(Layer)*dx(Layer)+rk(Layer+1)*dx(Layer+1))/(capavg*dyn*dyn)
AMat(Node,Node) =-amatx-amatxx-2.0d0*amaty
AMat(Node,Node+1) = 2.0d0*amaty
AMat(Node,Node-NumOfPerpendNodes) = amatx
AMat(Node,Node+NumOfPerpendNodes) = amatxx
DO NodeInRow = 2, NumOfPerpendNodes-1
Node2 = Node + NodeInRow - 1
AMat(Node2,Node2-1) = amaty
AMat(Node2,Node2) =-amatx-amatxx-2.0d0*amaty
AMat(Node2,Node2+1) = amaty
AMat(Node2,Node2-NumOfPerpendNodes) = amatx
AMat(Node2,Node2+NumOfPerpendNodes) = amatxx
END DO
Node2 = Node - 1 + NumOfPerpendNodes
AMat(Node2,Node2) =-amatx-amatxx-2.0d0*amaty
AMat(Node2,Node2-1) = 2.0d0*amaty
AMat(Node2,Node2-NumOfPerpendNodes) = amatx
AMat(Node2,Node2+NumOfPerpendNodes) = amatxx
IF (Node == NodeSource) BMat(3) = 2.0d0*REAL(NumOfPerpendNodes-1,r64)/capavg
NodeInLayer = 0
Layer = Layer + 1
END IF
NodeInLayer = NodeInLayer + 1
END DO
! LAST ROW OF NODES: Like the first row of nodes, this row is exposed to a boundary
! condition and thus has slightly modified nodal equations.
! As with the 1-D solution, we are accounting for the thermal mass
! of the half-node at the surface by adding it to the first row
! of interior nodes at both sides of the construction. This is not
! exact, but it does take all of the thermal mass into account.
amatx = amatx/1.5d0
amaty = amaty/1.5d0
Node = rcmax + 1 - NumOfPerpendNodes
AMat(Node,Node) =-2.0d0*(amatx+amaty)
AMat(Node,Node+1) = 2.0d0*amaty
AMat(Node,Node-NumOfPerpendNodes) = amatx
DO Node = rcmax+2-NumOfPerpendNodes, rcmax-1
AMat(Node,Node-1) = amaty
AMat(Node,Node) =-2.0d0*(amatx+amaty)
AMat(Node,Node+1) = amaty
AMat(Node,Node-NumOfPerpendNodes) = amatx
END DO
AMat(rcmax,rcmax) =-2.0d0*(amatx+amaty)
AMat(rcmax,rcmax-1) = 2.0d0*amaty
AMat(rcmax,rcmax-NumOfPerpendNodes) = amatx
BMat(2) = amatx
CMat(1) =-rk(1)/dx(1)/REAL(NumOfPerpendNodes-1,r64)
CMat(2) = rk(LayersInConstruct)/dx(LayersInConstruct)/REAL(NumOfPerpendNodes-1,r64)
DMat(1) = rk(1)/dx(1)/REAL(NumOfPerpendNodes-1,r64)
DMat(2) =-rk(LayersInConstruct)/dx(LayersInConstruct)/REAL(NumOfPerpendNodes-1,r64)
END IF
! Calculation of the CTFs based on the state space
! method. This process involves finding the exponential
! and inverse of AMat and using these results to
! determine the CTFs. The Gammas are an intermediate
! calculations which are necessary before the CTFs can
! be computed in TransFuncCoeffs.
CALL DisplayNumberAndString(ConstrNum,'Calculating CTFs for "'//TRIM(Construct(ConstrNum)%Name)//'", Construction #')
! CALL DisplayNumberandString(ConstrNum,'Matrix exponential for Construction #')
CALL CalculateExponentialMatrix(Construct(ConstrNum)%CTFTimeStep) ! Compute exponential of AMat
! CALL DisplayNumberandString(ConstrNum,'Invert Matrix for Construction #')
CALL CalculateInverseMatrix ! Compute inverse of AMat
! CALL DisplayNumberandString(ConstrNum,'Gamma calculation for Construction #')
CALL CalculateGammas(Construct(ConstrNum)%CTFTimeStep,Construct(ConstrNum)%SolutionDimensions)
! Compute "gamma"s from AMat, AExp, and AInv
! CALL DisplayNumberandString(ConstrNum,'Compute CTFs for Construction #')
CALL CalculateCTFs(Construct(ConstrNum)%NumCTFTerms,Construct(ConstrNum)%SolutionDimensions) ! Compute CTFs
! Now check to see if the number of transfer functions
! is greater than MaxCTFTerms. If it is, then increase the
! time step and the number of history terms and
! recalculate. Whether or not it will be necessary to
! recalculate the CTFs is controlled by this DO WHILE
! loop and the logical CTFConvrg.
CTFConvrg = .TRUE. ! Assume solution convergence
! If too many terms, then solution did not converge. Increase the
! number of histories and the time step. Reset CTFConvrg to continue
! the DO loop.
IF (Construct(ConstrNum)%NumCTFTerms > (MaxCTFTerms-1)) THEN
Construct(ConstrNum)%NumHistories = Construct(ConstrNum)%NumHistories+1
Construct(ConstrNum)%CTFTimeStep = Construct(ConstrNum)%CTFTimeStep &
+TimeStepZone
CTFConvrg = .FALSE.
END IF
! If the number of terms is okay, then do a further check on the summation of
! the various series summations. In theory, Sum(Xi) = Sum(Yi) = Sum(Zi). If
! this is not the case, then the terms have not reached a valid solution, and
! we need to increase the number of histories and the time step as above.
IF (CTFConvrg) THEN
SumXi = s0(2,2)
SumYi = s0(2,1)
SumZi = s0(1,1)
DO HistTerm = 1, Construct(ConstrNum)%NumCTFTerms
SumXi = SumXi + s(HistTerm,2,2)
SumYi = SumYi + s(HistTerm,2,1)
SumZi = SumZi + s(HistTerm,1,1)
END DO
SumXi = ABS(SumXi)
SumYi = ABS(SumYi)
SumZi = ABS(SumZi)
BiggestSum = MAX(SumXi,SumYi,SumZi)
IF (BiggestSum > 0.0d0) THEN
IF ( ((ABS(SumXi-SumYi)/BiggestSum) > MaxAllowedCTFSumError) .OR. &
((ABS(SumZi-SumYi)/BiggestSum) > MaxAllowedCTFSumError) ) THEN
Construct(ConstrNum)%NumHistories = Construct(ConstrNum)%NumHistories + 1
Construct(ConstrNum)%CTFTimeStep = Construct(ConstrNum)%CTFTimeStep + TimeStepZone
CTFConvrg = .FALSE.
END IF
ELSE ! Something terribly wrong--the surface has no CTFs, not even an R-value
CALL ShowFatalError('Illegal construction definition, no CTFs calculated for '//TRIM(Construct(ConstrNum)%Name))
END IF
END IF
! Once the time step has reached a certain point, it is highly likely that
! there is either a problem with the input or the solution. This should
! be extremely rare since other checks should flag most bad user input.
! Thus, if the time step reaches a certain point, error out and let the
! user know that something needs to be checked in the input file.
IF (Construct(ConstrNum)%CTFTimeStep >= MaxAllowedTimeStep) THEN
CALL ShowSevereError('CTF calculation convergence problem for Construction="'//TRIM(Construct(ConstrNum)%Name)//'".')
CALL ShowContinueError('...with Materials (outside layer to inside)')
CALL ShowContinueError('(outside)="'//trim(Material(Construct(ConstrNum)%LayerPoint(1))%Name)//'"')
DO Layer=2,Construct(ConstrNum)%TotLayers
IF (Layer /= Construct(ConstrNum)%TotLayers) THEN
CALL ShowContinueError('(next)="'//trim(Material(Construct(ConstrNum)%LayerPoint(Layer))%Name)//'"')
ELSE
CALL ShowContinueError('(inside)="'//trim(Material(Construct(ConstrNum)%LayerPoint(Layer))%Name)//'"')
ENDIF
ENDDO
CALL ShowContinueError('The Construction report will be produced. This will show more '// &
'details on Constructions and their materials.')
CALL ShowContinueError('Attempts will be made to complete the CTF process but the report may be incomplete.')
CALL ShowContinueError('Constructs reported after this construction may appear to have all 0 CTFs.')
CALL ShowContinueError('The potential causes of this problem are related to the input for the construction')
CALL ShowContinueError('listed in the severe error above. The CTF calculate routine is unable to come up')
CALL ShowContinueError('with a series of CTF terms that have a reasonable time step and this indicates an')
CALL ShowContinueError('error. Check the definition of this construction and the materials that make up')
CALL ShowContinueError('the construction. Very thin, highly conductive materials may cause problems.')
CALL ShowContinueError('This may be avoided by ignoring the presence of those materials since they probably')
CALL ShowContinueError('do not effect the heat transfer characteristics of the construction. Highly')
CALL ShowContinueError('conductive or highly resistive layers that are alternated with high mass layers')
CALL ShowContinueError('may also result in problems. After confirming that the input is correct and')
CALL ShowContinueError('realistic, the user should contact the EnergyPlus support team.')
DoCTFErrorReport=.true.
ErrorsFound=.true.
EXIT
! CALL ShowFatalError('Program terminated for reasons listed (InitConductionTransferFunctions) ')
END IF
END DO ! ... end of CTF calculation loop.
END IF ! ... end of IF block for non-reversed constructs.
ELSE ! Construct has only resistive layers (no thermal mass).
! CTF calculation not necessary, overall resistance
! (R-value) is all that is needed.
! Set time step for construct to user time step and the number of
! inter-time step interpolations to 1
Construct(ConstrNum)%CTFTimeStep = TimeStepZone
Construct(ConstrNum)%NumHistories = 1
Construct(ConstrNum)%NumCTFTerms = 1
s0(1,1) = cnd ! CTFs for current time
s0(1,2) = -cnd ! step are set to the
s0(2,1) = cnd ! overall conductance
s0(2,2) = -cnd ! of the construction.
ALLOCATE(e(1))
e=0.0D0
ALLOCATE(s(1,2,2))
s=0.0D0
s(1,1,1) = 0.0D0 ! CTF temperature
s(1,1,2) = 0.0D0 ! and flux
s(1,2,1) = 0.0D0 ! history terms
s(1,2,2) = 0.0D0 ! are all
e(1) = 0.0D0 ! zero.
IF (Construct(ConstrNum)%SourceSinkPresent) THEN
CALL ShowSevereError('Sources/sinks not allowed in purely resistive constructions --> '//Construct(ConstrNum)%Name)
ErrorsFound = .TRUE.
END IF
RevConst = .FALSE. ! In the code that follows, handle a resistive
! layer as a non-reversed construction.
END IF ! ... end of resistive construction IF block.
! Transfer the CTFs to the storage arrays for all non-reversed
! constructions. This transfer was done earlier in the routine for
! reversed constructions.
IF (.NOT.RevConst) THEN ! If this is either a new construction or a non-
! reversed construction, the CTFs must be stored
! in the proper arrays. If this is a reversed
! construction, nothing further needs to be done.
! Copy the CTFs into the storage arrays, converting them back to SI
! units in the process. First the "zero" terms and then the history terms...
Construct(ConstrNum)%CTFOutside(0) = s0(1,1)*CFU
Construct(ConstrNum)%CTFCross(0) = s0(2,1)*CFU
Construct(ConstrNum)%CTFInside(0) = -s0(2,2)*CFU
IF (Construct(ConstrNum)%SourceSinkPresent) THEN
! QTFs...
Construct(ConstrNum)%CTFSourceOut(0) = s0(1,3)
Construct(ConstrNum)%CTFSourceIn(0) = s0(2,3)
! QTFs for temperature calculation at source/sink location
Construct(ConstrNum)%CTFTSourceOut(0) = s0(3,1)
Construct(ConstrNum)%CTFTSourceIn(0) = s0(3,2)
Construct(ConstrNum)%CTFTSourceQ(0) = s0(3,3)/CFU
IF (Construct(ConstrNum)%TempAfterLayer /= 0) THEN
! QTFs for user specified interior temperature calculations...
Construct(ConstrNum)%CTFTUserOut(0) = s0(4,1)
Construct(ConstrNum)%CTFTUserIn(0) = s0(4,2)
Construct(ConstrNum)%CTFTUserSource(0) = s0(4,3)/CFU
END IF
END IF
DO HistTerm = 1, Construct(ConstrNum)%NumCTFTerms
! "REGULAR" CTFs...
Construct(ConstrNum)%CTFOutside(HistTerm) = s(HistTerm,1,1)*CFU
Construct(ConstrNum)%CTFCross(HistTerm) = s(HistTerm,2,1)*CFU
Construct(ConstrNum)%CTFInside(HistTerm) = -s(HistTerm,2,2)*CFU
IF (HistTerm /= 0) Construct(ConstrNum)%CTFFlux(HistTerm) = -e(HistTerm)
IF (Construct(ConstrNum)%SourceSinkPresent) THEN
! QTFs...
Construct(ConstrNum)%CTFSourceOut(HistTerm) = s(HistTerm,1,3)
Construct(ConstrNum)%CTFSourceIn(HistTerm) = s(HistTerm,2,3)
! QTFs for temperature calculation at source/sink location
Construct(ConstrNum)%CTFTSourceOut(HistTerm) = s(HistTerm,3,1)
Construct(ConstrNum)%CTFTSourceIn(HistTerm) = s(HistTerm,3,2)
Construct(ConstrNum)%CTFTSourceQ(HistTerm) = s(HistTerm,3,3)/CFU
IF (Construct(ConstrNum)%TempAfterLayer /= 0) THEN
! QTFs for user specified interior temperature calculations...
Construct(ConstrNum)%CTFTUserOut(HistTerm) = s(HistTerm,4,1)
Construct(ConstrNum)%CTFTUserIn(HistTerm) = s(HistTerm,4,2)
Construct(ConstrNum)%CTFTUserSource(HistTerm) = s(HistTerm,4,3)/CFU
END IF
END IF
END DO
END IF ! ... end of the reversed construction IF block.
Construct(ConstrNum)%UValue = cnd*CFU
IF (ALLOCATED(AExp)) DEALLOCATE(AExp)
IF (ALLOCATED(AMat)) DEALLOCATE(AMat)
IF (ALLOCATED(AInv)) DEALLOCATE(AInv)
IF (ALLOCATED(IdenMatrix)) DEALLOCATE(IdenMatrix)
IF (ALLOCATED(e)) DEALLOCATE(e)
IF (ALLOCATED(Gamma1)) DEALLOCATE(Gamma1)
IF (ALLOCATED(Gamma2)) DEALLOCATE(Gamma2)
IF (ALLOCATED(s)) DEALLOCATE(s)
END DO ! ... end of construction loop.
CALL ReportCTFs(DoCTFErrorReport)
IF (ErrorsFound) THEN
CALL ShowFatalError('Program terminated for reasons listed (InitConductionTransferFunctions)')
ENDIF
RETURN
END SUBROUTINE InitConductionTransferFunctions