ConductionTransferFunctionCalc.f90 Source File

This File Depends On

sourcefile~~conductiontransferfunctioncalc.f90~~EfferentGraph sourcefile~conductiontransferfunctioncalc.f90 ConductionTransferFunctionCalc.f90 sourcefile~general.f90 General.f90 sourcefile~general.f90->sourcefile~conductiontransferfunctioncalc.f90 sourcefile~dataheatbalance.f90 DataHeatBalance.f90 sourcefile~general.f90->sourcefile~dataheatbalance.f90 sourcefile~dataenvironment.f90 DataEnvironment.f90 sourcefile~general.f90->sourcefile~dataenvironment.f90 sourcefile~dataconversions.f90 DataConversions.f90 sourcefile~dataconversions.f90->sourcefile~conductiontransferfunctioncalc.f90 sourcefile~dataprecisionglobals.f90 DataPrecisionGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~conductiontransferfunctioncalc.f90 sourcefile~dataprecisionglobals.f90->sourcefile~general.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataconversions.f90 sourcefile~dataglobals.f90 DataGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataglobals.f90 sourcefile~datainterfaces.f90 DataInterfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datainterfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataheatbalance.f90 sourcefile~datahvacglobals.f90 DataHVACGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datahvacglobals.f90 sourcefile~datasurfaces.f90 DataSurfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasurfaces.f90 sourcefile~dataipshortcuts.f90 DataIPShortCuts.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataipshortcuts.f90 sourcefile~dataruntimelanguage.f90 DataRuntimeLanguage.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataruntimelanguage.f90 sourcefile~inputprocessor.f90 InputProcessor.f90 sourcefile~dataprecisionglobals.f90->sourcefile~inputprocessor.f90 sourcefile~datavectortypes.f90 DataVectorTypes.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datavectortypes.f90 sourcefile~databsdfwindow.f90 DataBSDFWindow.f90 sourcefile~dataprecisionglobals.f90->sourcefile~databsdfwindow.f90 sourcefile~datasizing.f90 DataSizing.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasizing.f90 sourcefile~datasystemvariables.f90 DataSystemVariables.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasystemvariables.f90 sourcefile~datacomplexfenestration.f90 DataComplexFenestration.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datacomplexfenestration.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataenvironment.f90 sourcefile~dataequivalentlayerwindow.f90 DataEquivalentLayerWindow.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataequivalentlayerwindow.f90 sourcefile~dataglobals.f90->sourcefile~conductiontransferfunctioncalc.f90 sourcefile~dataglobals.f90->sourcefile~general.f90 sourcefile~dataglobals.f90->sourcefile~dataheatbalance.f90 sourcefile~dataglobals.f90->sourcefile~datahvacglobals.f90 sourcefile~dataglobals.f90->sourcefile~datasurfaces.f90 sourcefile~dataglobals.f90->sourcefile~dataipshortcuts.f90 sourcefile~dataglobals.f90->sourcefile~dataruntimelanguage.f90 sourcefile~dataglobals.f90->sourcefile~inputprocessor.f90 sourcefile~dataglobals.f90->sourcefile~databsdfwindow.f90 sourcefile~sortandstringutilities.f90 SortAndStringUtilities.f90 sourcefile~dataglobals.f90->sourcefile~sortandstringutilities.f90 sourcefile~dataoutputs.f90 DataOutputs.f90 sourcefile~dataglobals.f90->sourcefile~dataoutputs.f90 sourcefile~dataglobals.f90->sourcefile~datasizing.f90 sourcefile~dataglobals.f90->sourcefile~datacomplexfenestration.f90 sourcefile~dataglobals.f90->sourcefile~dataenvironment.f90 sourcefile~dataglobals.f90->sourcefile~dataequivalentlayerwindow.f90 sourcefile~datainterfaces.f90->sourcefile~conductiontransferfunctioncalc.f90 sourcefile~datainterfaces.f90->sourcefile~general.f90 sourcefile~datainterfaces.f90->sourcefile~dataheatbalance.f90 sourcefile~datainterfaces.f90->sourcefile~dataruntimelanguage.f90 sourcefile~datainterfaces.f90->sourcefile~inputprocessor.f90 sourcefile~datainterfaces.f90->sourcefile~dataenvironment.f90 sourcefile~dataheatbalance.f90->sourcefile~conductiontransferfunctioncalc.f90 sourcefile~datahvacglobals.f90->sourcefile~general.f90 sourcefile~datasurfaces.f90->sourcefile~general.f90 sourcefile~datasurfaces.f90->sourcefile~dataheatbalance.f90 sourcefile~dataipshortcuts.f90->sourcefile~general.f90 sourcefile~dataipshortcuts.f90->sourcefile~inputprocessor.f90 sourcefile~dataruntimelanguage.f90->sourcefile~general.f90 sourcefile~datastringglobals.f90 DataStringGlobals.f90 sourcefile~datastringglobals.f90->sourcefile~general.f90 sourcefile~datastringglobals.f90->sourcefile~inputprocessor.f90 sourcefile~datastringglobals.f90->sourcefile~datasystemvariables.f90 sourcefile~inputprocessor.f90->sourcefile~general.f90 sourcefile~inputprocessor.f90->sourcefile~dataheatbalance.f90 sourcefile~datavectortypes.f90->sourcefile~dataheatbalance.f90 sourcefile~datavectortypes.f90->sourcefile~datasurfaces.f90 sourcefile~datavectortypes.f90->sourcefile~databsdfwindow.f90 sourcefile~databsdfwindow.f90->sourcefile~dataheatbalance.f90 sourcefile~databsdfwindow.f90->sourcefile~datasurfaces.f90 sourcefile~sortandstringutilities.f90->sourcefile~inputprocessor.f90 sourcefile~dataoutputs.f90->sourcefile~inputprocessor.f90 sourcefile~datasizing.f90->sourcefile~inputprocessor.f90 sourcefile~datasystemvariables.f90->sourcefile~inputprocessor.f90 sourcefile~datacomplexfenestration.f90->sourcefile~dataheatbalance.f90 sourcefile~dataenvironment.f90->sourcefile~dataheatbalance.f90 sourcefile~dataequivalentlayerwindow.f90->sourcefile~dataheatbalance.f90
Help

Files Dependent On This One

sourcefile~~conductiontransferfunctioncalc.f90~~AfferentGraph sourcefile~conductiontransferfunctioncalc.f90 ConductionTransferFunctionCalc.f90 sourcefile~heatbalancemanager.f90 HeatBalanceManager.f90 sourcefile~conductiontransferfunctioncalc.f90->sourcefile~heatbalancemanager.f90 sourcefile~ecoroof.f90 EcoRoof.f90 sourcefile~conductiontransferfunctioncalc.f90->sourcefile~ecoroof.f90 sourcefile~sizingmanager.f90 SizingManager.f90 sourcefile~heatbalancemanager.f90->sourcefile~sizingmanager.f90 sourcefile~simulationmanager.f90 SimulationManager.f90 sourcefile~heatbalancemanager.f90->sourcefile~simulationmanager.f90 sourcefile~heatbalancesurfacemanager.f90 HeatBalanceSurfaceManager.f90 sourcefile~ecoroof.f90->sourcefile~heatbalancesurfacemanager.f90 sourcefile~sizingmanager.f90->sourcefile~simulationmanager.f90 sourcefile~energyplus.f90 EnergyPlus.f90 sourcefile~simulationmanager.f90->sourcefile~energyplus.f90 sourcefile~utilityroutines.f90 UtilityRoutines.f90 sourcefile~simulationmanager.f90->sourcefile~utilityroutines.f90 sourcefile~heatbalancesurfacemanager.f90->sourcefile~heatbalancemanager.f90 sourcefile~heatbalancesurfacemanager.f90->sourcefile~simulationmanager.f90
Help


Source Code

MODULE ConductionTransferFunctionCalc

          ! MODULE INFORMATION:
          !       AUTHOR         Rick Strand
          !       DATE WRITTEN   January 1997
          !       MODIFIED       June-July 2000, RKS
          !       RE-ENGINEERED  June 1996, February 1997, August 1997, RKS
          !       RE-ENGINEERED  November 1999, LKL

          ! PURPOSE OF THIS MODULE:
          ! This module calculates the conduction transfer functions (CTFs) for
          ! all building constructions.

          ! METHODOLOGY EMPLOYED:
          ! This subroutine uses the state space method of calculating CTFs.
          ! The state space method involves imposing a finite difference grid
          ! to a solution space (i.e., a building construction, inside to
          ! outside surface).  The finite difference grid is only used to
          ! derive a system of differential equations.  This first order
          ! system can then be solved using matrix algebra.  In this
          ! implementation of the state space method, a conversion from the
          ! internal units of EnergyPlus (SI) to English units is used, This
          ! is done due to observations made by Russ Taylor that the solution
          ! method was not as stable numerically when SI units were used.

          ! REFERENCES:
          ! While there are several important references on the state space
          ! method, the most definitive reference which includes details on
          ! implementing the state space method is:
          ! Seem, J.E.  1987.  Modeling of Heat Transfer in Buildings, Ph.D.
          ! Dissertation, Department of Mechanical Engineering, University of
          ! Wisconsin-Madison.

          ! OTHER NOTES:
          ! na

          ! USE STATEMENTS:
USE DataPrecisionGlobals
USE DataGlobals
USE DataHeatBalance ! This is the Heat balance super block data-only module
USE DataInterfaces


IMPLICIT NONE    ! Enforce explicit typing of all variables

PRIVATE
          ! MODULE PARAMETER DEFINITIONS:
!INTEGER, PRIVATE, PARAMETER :: MaxTotNodes = 75   ! Maximum total number of
          ! nodes per construction.  This limit is a compromise between faster,
          ! less accurate solutions and slower, possibly more accurate answers.

INTEGER, PARAMETER :: NumOfPerpendNodes = 7     ! Number of nodes in the direction
          ! perpendicular to the main direction of heat transfer.  This is only used
          ! when a two-dimensional solution has been requested for a construction
          ! with a heat source/sink.

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! MODULE VARIABLE DECLARATIONS:
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:)        :: AExp   ! Exponential of AMat
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:)        :: AInv   ! Inverse of AMat
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:)        :: AMat   ! "A" matrix from Seem's dissertation
                                                                        ! (constant coefficients of linear system)
REAL(r64), PRIVATE, DIMENSION(3)                       :: BMat   ! "B" matrix of state space method (non-zero elements)
REAL(r64), PRIVATE, DIMENSION(2)                       :: CMat   ! "C" matrix of state space method (non-zero elements)
REAL(r64), PRIVATE, DIMENSION(2)                       :: DMat   ! "D" matrix of state space method (non-zero elements)
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:)          :: e      ! Coefficients for the surface flux history term
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:)        :: Gamma1 ! Intermediate calculation array corresponding to a term
                                                                        ! in Seem's dissertation
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:)        :: Gamma2 ! Intermediate calculation array corresponding to a term
                                                                        ! in Seem's dissertation
INTEGER,   PRIVATE                                     :: NodeSource ! Node at which a source or sink is present
INTEGER,   PRIVATE                                     :: NodeUserTemp ! Node where user wishes to calculate a temperature
                                                                              ! (for constructions with sources/sinks only)
INTEGER,   PRIVATE                                     :: rcmax  ! Total number of nodes in the construct (<= MaxTotNodes)
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:,:)      :: s      ! Coefficients for the surface temperature history terms
REAL(r64), PRIVATE, DIMENSION(4,3)                     :: s0     ! Coefficients for the current surface temperature terms
REAL(r64), PRIVATE                                     :: TinyLimit
REAL(r64), PRIVATE, ALLOCATABLE, DIMENSION(:,:)        :: IdenMatrix  ! Identity Matrix

          ! SUBROUTINE SPECIFICATIONS FOR MODULE ConductionTransferFunctionCalc

PUBLIC  InitConductionTransferFunctions
PRIVATE CalculateExponentialMatrix
PRIVATE CalculateInverseMatrix
PRIVATE CalculateGammas
PRIVATE CalculateCTFs
PRIVATE ReportCTFs

CONTAINS

          ! MODULE SUBROUTINES:

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

SUBROUTINE CalculateExponentialMatrix(delt)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Russ Taylor
          !       DATE WRITTEN   June 1990
          !       MODIFIED       Dec 1995, Apr 1996, RKS; June 2000 RKS
          !       RE-ENGINEERED  June 1996, RKS; Nov 1999, LKL;

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine computes the exponential matrix exp(AMat*delt) for
          ! use in the state space method for the calculation of CTFs.

          ! METHODOLOGY EMPLOYED:
          ! Uses the method of Taylor expansion combined with scaling and
          ! squaring to most efficiently compute the exponential matrix.  The
          ! steps in the procedure are outlined in Seem's dissertation in
          ! Appendix A, page 128.  Exponential matrix multiplication modified
          ! to take advantage of the characteristic form of AMat.  AMat starts
          ! out as a tri-diagonal matrix.  Each time AMat is raised to a higher
          ! power two extra non-zero diagonals are added.  ExponMatrix now
          ! recognizes this.  This should speed up the calcs somewhat.  Also, a
          ! new cut-off criteria based on the significant figures of double-
          ! precision variables has been added.  The main loop for higher powers
          ! of AMat is now stopped whenever these powers of AMat will no longer
          ! add to the summation (AExp) instead ofstopping potentially at the
          ! artifical limit of AMat**100.

          ! 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.

          ! USE STATEMENTS:
          ! none

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine


          ! SUBROUTINE ARGUMENT DEFINITIONS:
  REAL(r64) :: delt  ! Time step of the resulting CTFs

          ! SUBROUTINE PARAMETER DEFINITIONS:

  REAL(r64), PARAMETER :: DPLimit = 1.0d-20
          ! This argument is nice, but not sure it's accurate -- LKL Nov 1999.
          ! Parameter set to the significant figures limit of double
          ! precision variables plus a safety factor.- The argument for setting this parameter to 1E-20 involves the
          ! number of significant figures for REAL(r64) variables which is 16 and the largest power to which
          ! AMat will be raised which is 100.  This would be a factor of 1E-18.  A factor of "safety" of another 100
          ! arrives at the value chosen.  It is argued that if one number is 1E-16 larger than a second number, then
          ! adding the second to the first will not effect the first.  However, on the conservative side, there could
          ! be up to 100 numbers which might, added together, still could effect the original number.  Each
          ! successive power of AMat will have terms smaller than the previous power.  Thus, when the ratio between
          ! the terms of the latest power of AMat and the total (AExp) is less than DPLim, all further powers of
          ! AMat will have absolutely no effect on the REAL(r64) value of AExp.  Thus, there is no need to
          ! continue the calculation.  In effect, AExp has "converged".  In REAL(r64)ity, 1E-16 would probably guarantee
          ! convergence since AMat terms drop off quickly, but the extra powers allows for differences between
          ! computer platforms.

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64)                         :: AMatRowNorm      ! Row norm for AMat
  REAL(r64)                         :: AMatRowNormMax   ! Largest row norm for AMat
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMat1    ! AMat factored by (delt/2^k)
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMato    ! AMat raised to the previous power (power of AMat1-1)
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMatN    ! Current value of AMat raised to power n (n = 1,2...)
  LOGICAL                                       :: Backup   ! Used when numerics get to small in Exponentiation
  REAL(r64)                                     :: CheckVal ! Used to avoid possible overflow from Double->REAL(r64)->Integer
  REAL(r64)                         :: fact     ! Intermediate calculation variable (delt/2^k)
  INTEGER                                       :: i        ! Loop counter
  INTEGER                                       :: ic       ! Loop counter
  INTEGER                                       :: ict      ! Loop counter
  INTEGER                                       :: idm      ! Loop counter
  INTEGER                                       :: ir       ! Loop counter
  INTEGER                                       :: isq      ! Loop counter
  INTEGER                                       :: j        ! Loop counter
  INTEGER                                       :: k        ! Power of 2 which is used to factor AMat
  INTEGER                                       :: l        ! Theoretical power to which the A matrix must be
                                                            ! raised to accurately calculate the exponential matrix
  LOGICAL                                       :: SigFigLimit      ! Significant figure limit logical, true
                                                                    ! when exponential calculation loop can be exited (i.e.
                                                                    ! the significant figure limit for REAL(r64)
                                                                    ! variables reached)

          ! FLOW:
  ALLOCATE(AMat1(rcmax,rcmax))
  ALLOCATE(AMato(rcmax,rcmax))
  ALLOCATE(AMatN(rcmax,rcmax))

  ! Subroutine initializations.  AMat is assigned to local variable AMat1 to
  ! avoid the corruption of the original AMat 2-d array.
  AMat1=AMat

  !  Other arrays are initialized to zero.
  AExp=0.0D0
  AMato=0.0D0
  AMatN=0.0D0


          ! Step 1, page 128 (Seem's thesis):  Compute the matrix row norm.
          ! See equation (A.3) which states that the matrix row norm is the
          ! maximum summation of the elements in a row of AMat multiplied by
          ! the time step.

  AMatRowNormMax = 0.0D0            ! Start of Step 1 ...

  DO i = 1,rcmax

    AMatRowNorm = 0.0D0
    DO j = 1,rcmax
      AMatRowNorm = AMatRowNorm+ABS(AMat1(i,j))
    END DO

    AMatRowNorm = AMatRowNorm*delt

    AMatRowNormMax=MAX(AMatRowNormMax,AMatRowNorm)

  END DO                            ! ... end of Step 1.

          ! Step 2, page 128:  Find smallest integer k such that
          ! AMatRowNormMax< = 2^k

  k = INT(LOG(AMatRowNormMax)/LOG(2.0D0))+1 !Objexx:Num Handle AMatRowNormMax=0

          ! Step 3, page 128:  Divide (AMat*delt) by 2^k.  This section of code
          ! takes advantage of the fact that AMat is tridiagonal.  Thus, it
          ! only factors the elements of the AMat that are known to be non-zero.

  fact  = delt/(2.0d0**k)   ! Start of Step 3 ...
  AMat1 = AMat1*fact        ! ... end of Step 3.

          ! Step 4, page 128:  Calculate l, the highest power to which AMat
          ! must be taken theoretically to accurately calculate its exponential.
          ! This is based on a paper by Cadzow and Martens ("Discrete-Time and
          ! Computer Control Systems",Prentice-Hall, pp. 389-390, 1970).  This
          ! number is now used as the maximum power to which AMat must be
          ! raised in order to calculate the exponential matrix.  A new cut-off
          ! criteria based on the number of significant figures in a double-
          ! precision variable is used as a more practical limit on the
          ! exponentiation algorithm.

  CheckVal = MIN(3.0d0*AMatRowNormMax+6.d0,100.0d0)
  l        = INT(CheckVal)

          ! Step 5, page 128:  Calculate the exponential.  First, add the
          ! linear term to the identity matrix.
   AExp = AMat1 + IdenMatrix        ! Start of Step 5 ...

          ! Now, add successive terms to the expansion as per the standard
          ! exponential formula.  AMato contains the last "power" of AMat
          ! which saves the program from having to remultiply the entire power
          ! of AMat each time.  Since this is still the linear power of AMat,
          ! AMat1 is still tridiagonal in nature.
  AMato = AMat1

  i = 1   ! Initialize the counter for the following DO loop

          ! The following DO WHILE loop continues to raise AMat to successive
          ! powers and add it to the exponential matrix (AExp).
  DO WHILE (i < l) ! Begin power raising loop ...

    i = i+1               ! Increment the loop counter
    SigFigLimit = .true.  ! Set the significant factor limit flag

    DO ir = 1, rcmax    ! Begin matrix multiplication loop ...
          ! The following matrix multiplication could be "optimized" since
          ! for one-dimensional heat transfer AMat is 3-diagonal, AMat squared
          ! is 5-diagonal, etc.  However, the code can be much simpler if we
          ! ignore this fact and just do a generic matrix multiplication.
          ! For 2-D heat transfer, the number of off-diagonal non-zero terms
          ! is slightly more complicated as well.
      DO ic = 1, rcmax
        AMatN(ir,ic) = 0.0d0
        DO ict = 1, rcmax
          ! Make sure the next term won't cause an underflow.  If it will end up being
          ! so small as to go below TinyLimit, then ignore it since it won't add anything
          ! to AMatN anyway.
          IF (ABS(AMat1(ict,ic)) > TinyLimit) THEN
            IF (ABS(AMato(ir,ict)) > ABS(REAL(i,r64)*TinyLimit/AMat1(ict,ic))) &
              AMatN(ir,ic) = AMatN(ir,ic)+AMato(ir,ict)*AMat1(ict,ic)/REAL(i,r64)
          END IF
        END DO
      END DO
    END DO          ! ... end of matrix multiplication loop.

          ! Update AMato and AExp matrices
    AMato = AMatN
    AExp  = AExp + AMato

          ! The next DO loop tests the significant figures limit criteria to
          ! see if any values in AExp are still changing appreciably.
    DO ir = 1, rcmax
      DO ic = 1, rcmax
          ! Test of limit criteria:
        IF (ABS(AExp(ir,ic)) > TinyLimit) THEN  ! Next line divides by AExp entry so it
                                                ! must be checked to avoid dividing by zero.
          ! If the ratio between any current element in the power
          ! of AMat and its corresponding element in AExp is
          ! greater than the number which might effect the overall
          ! exponential matrix based on stability criteria, then
          ! continue raising AMat to another power (SigFigLimit = false).

          IF (ABS(AMato(ir,ic)/AExp(ir,ic)) > DPLimit) THEN
            SigFigLimit = .FALSE.
            EXIT ! DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
          END IF

        ELSE        ! There are still elements of AExp which are zero, so
                    ! the raising of AMat to higher powers should continue.

          SigFigLimit = .FALSE.
          EXIT ! DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)

        END IF
      END DO
      IF (.NOT.SigFigLimit) EXIT ! DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
    END DO

          ! Compute next term, only if necessary.  If SigFigLimit is still true,
          ! then all of the new terms being added to AExp are too small to
          ! affect it.  Thus, there is no need to continue this do loop further.

    IF (SigFigLimit) i = 100  ! SigFigLimit is still true, set i to maximum possible
                              ! value of l (100).

  END DO            ! ... end of power raising loop and Step 5.

          ! Step 6, page 128:
          ! Square AExp "k times" to obtain the actual exponential matrix
          ! (remember that AExp was scaled earlier in this routine).

  DO isq = 1,k      ! Begin squaring DO loop and Step 6 ...

          ! Use AMato to store the old values of AExp
    AMato  = AExp
    Backup = .true.
    AExp   = 0.0D0

          ! Multiply the old value of AExp (AMato) by itself and store in AExp.
    DO ir = 1,rcmax
      DO ic = 1,rcmax
        DO idm = 1,rcmax
          IF (ABS(AMato(ir,idm)*AMato(idm,ic)) > TinyLimit) THEN
            AExp(ir,ic) = AExp(ir,ic)+AMato(ir,idm)*AMato(idm,ic)
            Backup=.false.
          END IF
        END DO
      END DO
    END DO
    ! Backup is true when every item of AExp didnt pass the TinyLimit test
    IF (Backup) THEN
      AExp=AMato
      EXIT
    END IF

  END DO            ! ... end of squaring loop and Step 6.

  DEALLOCATE(AMat1)
  DEALLOCATE(AMato)
  DEALLOCATE(AMatN)

  RETURN            ! AExp is now the exponential of AMat.

END SUBROUTINE CalculateExponentialMatrix

SUBROUTINE CalculateInverseMatrix

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Rick Strand
          !       DATE WRITTEN   Dec 1995
          !       MODIFIED       June 2000 RKS (made routine generic to allow for 2-D solutions)
          !       RE-ENGINEERED  June 1996, February 1997 RKS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine computes the inverse of AMat for use
          ! in the calculation of the CTFs.

          ! METHODOLOGY EMPLOYED:
          ! Uses row elimination to zero the off-diagonal terms of
          ! AMat while performing the same operations on another
          ! matrix which starts as the identity matrix.  Once AMat
          ! has been converted to an identity matrix(I), the other
          ! matrix which started as the I will then be the inverse
          ! of A.  This algorithm has been customized for a
          ! tri-diagonal matrix.

          ! REFERENCES:
          ! Any linear algebra test (this is a generic routine).

          ! USE STATEMENTS:
          ! none

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
          ! na

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMat1  ! Intermediate calculation matrix equivalent at first to AMat
  INTEGER                                       :: ic     ! Loop counter
  INTEGER                                       :: ir     ! Loop counter
  INTEGER                                       :: irr    ! Loop counter

          ! FLOW:

          ! Subroutine initializations ...
  ALLOCATE(AMat1(rcmax,rcmax))

  AMat1 = AMat          ! Set AMat1 = AMat to avoid AMat changes
  AInv  = IdenMatrix    ! Set AInv to Identity Matrix

          ! Use Gaussian elimination to zero all of the elements of AMat left
          ! of the diagonal.
          ! This DO loop will cycle through each of the rows in AMat except the
          ! last row which is handled later because it does not have to be used
          ! to eliminate any other rows.  The index ir is the current row
          ! number and also the column of the current diagonal element.

  DO ir = 1, rcmax-1    ! Begin forward elimination loop ...

          ! Factor all of the elements of the row being used to zero the next
          ! row in both AMat and AInv by the diagonal element of this row.
          ! We should only need to factor the elements to the right of the
          ! diagonal since those to the right of it should be zero.
    DO ic = ir+1, rcmax
      AMat1(ir,ic) = AMat1(ir,ic)/AMat1(ir,ir)
    END DO

          ! In the forward elimination process, all the elements in AInv to the
          ! right of the diagonal are zero so they do not need to be factored.
    DO ic = 1, ir
      AInv(ir,ic) = AInv(ir,ic)/AMat1(ir,ir)
    END DO

    AMat1(ir,ir) = 1.0d0      ! By definition, the diagonal of AMat is now 1.

          ! Use this factored row to eliminate the off-diagonal element of the
          ! rows below the current one (ir)...

    DO irr = ir+1, rcmax    ! Start of row reduction loop...

      DO ic = ir+1, rcmax
        AMat1(irr,ic) = AMat1(irr,ic) - AMat1(irr,ir)*AMat1(ir,ic)
      END DO

          ! Now, determine the effect on the next row of AInv.  Again, all of
          ! the elements in AInv to the right of the diagonal are zero, so they
          ! can be ignored.

      DO ic = 1,ir
        AInv(irr,ic) = AInv(irr,ic) - AMat1(irr,ir)*AInv(ir,ic)
      END DO

      AMat1(irr,ir) = 0.0d0 ! By definition, the element to the left of the
                            ! diagonal in the next row of AMat is now zero.

    END DO                  ! ...end of row reduction loop

  END DO            ! ... end of the forward elimination loop.

          ! Factor the last row of AInv by the current value of the last
          ! diagonal element of AMat. After this is done, all of the diagonal
          ! elements of AMat are unity and all of the elements in AMat left of
          ! the diagonal are zero.

  DO ic = 1,rcmax
    AInv(rcmax,ic) = AInv(rcmax,ic)/AMat1(rcmax,rcmax)
  END DO
  AMat1(rcmax,rcmax) = 1.0D0

          ! Now, use back substitution to eliminate the elements to the right
          ! of the diagonal in AMat.  The procedure is similar to the forward
          ! elimination process except that we only have to operate on AInv,
          ! though now all of the columns of AInv may be non-zero.

          ! This DO loop will cycle through the remaining rows which are not
          ! yet diagonalized in reverse order.  Note that the only effect on
          ! AMat is that the off-diagonal element is zeroed.  The diagonal
          ! (which has already been set to unity) is not effected by this row
          ! elimination process.
          ! In the following code ir is the column being zeroed and irr is the
          ! row being worked on

  DO ir = rcmax, 2, -1  ! Begin reverse elimination loop ...
    DO irr = 1, ir-1
      DO ic = 1,rcmax
        AInv(irr,ic) = AInv(irr,ic)-AMat1(irr,ir)*AInv(ir,ic)
      END DO
      AMat1(irr,ir) = 0.0d0
    END DO
  END DO                ! ... end of reverse elimination loop.

          ! At this point, AMat1 is equal to the identity matrix (I)
          ! and AInv is equal to the inverse of AMat.

  DEALLOCATE(AMat1)

  RETURN

END SUBROUTINE CalculateInverseMatrix

SUBROUTINE CalculateGammas(delt,SolutionDimensions)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Russ Taylor
          !       DATE WRITTEN   June 1990
          !       MODIFIED       na
          !       RE-ENGINEERED  July 1996, RKS

          ! PURPOSE OF THIS SUBROUTINE:
          ! Compute gammas as defined in Seem's dissertation.
          ! Runs as a subroutine of the conduction transfer
          ! function solver (InitializeCTFs).

          ! METHODOLOGY EMPLOYED:
          ! Determine the Gamma1 and Gamma2 based on the results
          ! from the ExponMatrix and InvertMatrix subroutines.
          ! This routine is specialized to take advantage of the
          ! fact that most of BMat consists of zeroes.

          ! REFERENCES:
          ! The state space method of calculating CTFs is
          ! outlined in the doctoral dissertation of John Seem,
          ! "Modeling of Heat Transfer in Buildings", Department
          ! of Mechanical Engineering, University of Wisconsin-
          ! Madison, 1987.

          ! USE STATEMENTS:
          ! none

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  REAL(r64), INTENT(IN) :: delt               ! Time increment in fraction of an hour
  INTEGER,          INTENT(IN) :: SolutionDimensions ! Integer relating whether a 1- or 2-D solution is required

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: ATemp    ! Intermediate variable equal to AExp - I
  INTEGER :: i          ! Loop counter
  INTEGER :: is1        ! Loop counter
  INTEGER :: j          ! Loop counter
  INTEGER :: SurfNode   ! Loop counter


          ! FLOW:

          ! Compute Gamma1 from equation (2.1.12) in Seem's dissertation which
          ! states that:  Gamma1  =  [AInv] * ([AExp]-[I]) * [BMat]
          ! noting that BMat contains only the non-zero values of the B Matrix.

  ALLOCATE(ATemp(rcmax,rcmax))
  ATemp  = AExp - IdenMatrix
  Gamma1 = 0.0d0

  DO i = 1,rcmax

    DO is1 = 1,rcmax

      IF (SolutionDimensions == 1) THEN
        Gamma1(i,1) = Gamma1(i,1)+AInv(i,is1)*ATemp(is1,1)*BMat(1)
        Gamma1(i,2) = Gamma1(i,2)+AInv(i,is1)*ATemp(is1,rcmax)*BMat(2)
      ELSE  ! SolutionDimensions = 2
        DO SurfNode = 1, NumOfPerpendNodes
          Gamma1(i,1) = Gamma1(i,1)+AInv(i,is1)*ATemp(is1,SurfNode)*BMat(1)
          Gamma1(i,2) = Gamma1(i,2)+AInv(i,is1)*ATemp(is1,rcmax+1-SurfNode)*BMat(2)
        END DO
      END IF

      IF (NodeSource > 0) THEN
        Gamma1(i,3) = Gamma1(i,3)+AInv(i,is1)*ATemp(is1,NodeSource)*BMat(3)
      END IF

    END DO

  END DO

  DEALLOCATE(ATemp)
          ! Compute Gamma2 from equation (2.1.13) in Seem's dissertation which
          ! states that:  Gamma2  =  [AInv] * ([Gamma1]/delt - [BMat])
          ! again noting that BMat contains only the non-zero values of B.
  Gamma2 = 0.0d0

  DO i = 1,rcmax

    DO j = 1,3

      DO is1 = 1,rcmax

        IF (SolutionDimensions == 1) THEN
          IF ( (j == 1) .AND. (is1 == 1) ) THEN
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt-BMat(1))
          ELSE IF ( (j == 2) .AND. (is1 == rcmax) ) THEN
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt-BMat(2))
          ELSE IF ( (j == 3) .AND. (is1 == NodeSource) ) THEN
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt-BMat(3))
          ELSE  ! the element of the actual BMat is zero
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt)
          END IF
        ELSE    ! SolutionDimensions = 2
          IF ( (j == 1) .AND. ((is1 >= 1).AND.(is1 <= NumOfPerpendNodes)) ) THEN
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt-BMat(1))
          ELSE IF ( (j == 2) .AND. ((is1 <= rcmax).AND.(is1 >= rcmax+1-NumOfPerpendNodes)) ) THEN
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt-BMat(2))
          ELSE IF ( (j == 3) .AND. (is1 == NodeSource) ) THEN
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt-BMat(3))
          ELSE  ! the element of the actual BMat is zero
            Gamma2(i,j) = Gamma2(i,j)+AInv(i,is1)*(Gamma1(is1,j)/delt)
          END IF
        END IF

      END DO

    END DO

  END DO

  RETURN  ! The calculation of Gamma1 and Gamma2 is now complete.

END SUBROUTINE CalculateGammas

SUBROUTINE CalculateCTFs(nrf,SolutionDimensions)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Russ Taylor
          !       DATE WRITTEN   June 1990
          !       MODIFIED       Apr 96, RKS, cosmetic, algorithm neutral changes
          !       RE-ENGINEERED  July 1996, RKS; Nov 1999, LKL (allocatable arrays)

          ! PURPOSE OF THIS SUBROUTINE:
          ! Subprogram to calculate the Sj and ej coefficients of Seem's
          ! dissertation.  Follows Seem's technique to compute the coefficients
          ! in order with minimum storage requirements.

          ! METHODOLOGY EMPLOYED:
          ! Combine the results of the ExponMatrix, InvertMatrix, and
          ! CalculateGammas routines together to arrive at the temperature
          ! coefficients (s, s0) and the heat flux history coefficients (e) of
          ! the CTFs.  The outline of this subroutine is based on step 5 of
          ! Seem's suggested implementation of the state space method found on
          ! pages 26+27 of his dissertation.

          ! REFERENCES:
          ! The state space method of calculating CTFs is outlined in the
          ! doctoral dissertation of John Seem, "Modeling of Heat Transfer in
          ! Buildings", Department of Mechanical Engineering, University of
          ! Wisconsin-Madison, 1987.  In particular, the equations used for
          ! these calculations are equations (2.1.24) through (2.1.26) in Seem's
          ! dissertation.

          ! USE STATEMENTS:
          ! none

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(OUT) :: nrf                ! Number of response factor terms
  INTEGER, INTENT(IN)  :: SolutionDimensions ! Integer relating whether a 1- or 2-D solution is required

          ! SUBROUTINE PARAMETER DEFINITIONS:
  REAL(r64), PARAMETER :: ConvrgLim = 1.0d-13   ! Convergence limit (ratio) for cutting off the calculation of further
                                                       ! CTFs.  This value was found to give suitable accuracy in IBLAST.

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE DUMMY VARIABLE DECLARATIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

  REAL(r64)                                :: avg       ! Intermediate calculation variable (average)
  LOGICAL                                              :: CTFConvrg ! Set after CTFs are calculated, based on whether there are
                                                                    ! too many CTFs terms
  INTEGER                                              :: i         ! Loop counter
  INTEGER                                              :: ic        ! Loop counter
  INTEGER                                              :: inum      ! Loop counter
  INTEGER                                              :: ir        ! Loop counter
  INTEGER                                              :: is        ! Loop counter
  INTEGER                                              :: is2       ! Loop counter
  INTEGER                                              :: j         ! Loop counter
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)        :: PhiR0     ! Product of Phi( = AExp) and R0 matrices from the state
                                                                    ! space method
  REAL(r64)                                :: rat       ! Intermediate calculation variable (ratio of flux history
                                                                    ! terms)
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)        :: Rnew      ! Current R matrix
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)        :: Rold      ! R matrix from the last iteration
  INTEGER                                              :: SurfNode  ! Loop counter (for nodes at a surface)
  REAL(r64)                                :: SurfNodeFac ! Multiplying factor applied to various surface nodes
  REAL(r64)                                :: trace     ! Trace of the product of Phi( = AExp) and R0

          ! FLOW:

          ! Subroutine initializations
  ALLOCATE(PhiR0(rcmax,rcmax))
  ALLOCATE(Rnew(rcmax,rcmax))
  ALLOCATE(Rold(rcmax,rcmax))
  PhiR0 = 0.0D0
  Rold  = 0.0D0

  s0   = 0.0D0
  s    = 0.0D0
  e    = 0.0D0
  Rnew = IdenMatrix ! Rnew initialized to the identity matrix

          ! Calculate Gamma1-Gamma2.  Gamma1 is not used by itself in the
          ! equations, only Gamma1-Gamma2.  Thus, reset Gamma1 to:
          ! Gamma1-Gamma2
  DO i = 1,rcmax
    DO j = 1,3
      Gamma1(i,j) = Gamma1(i,j)-Gamma2(i,j)
    END DO
  END DO

          ! Compute s0.  See Seem's thesis equation (2.1.24) which states that:
          ! s0  =  (CMat*R0*Gamma2) + (DMat)
          ! Note that for a two-dimensional solution, there is more than one
          ! node at the surface and the effect of each of these must be added
          ! together.
  IF (SolutionDimensions == 1) THEN
    s0(1,1) = CMat(1)*Gamma2(1,1) + DMat(1)
    s0(1,2) = CMat(1)*Gamma2(1,2)
    s0(1,3) = CMat(1)*Gamma2(1,3)
    s0(2,1) = CMat(2)*Gamma2(rcmax,1)
    s0(2,2) = CMat(2)*Gamma2(rcmax,2) + DMat(2)
    s0(2,3) = CMat(2)*Gamma2(rcmax,3)
  ELSE  ! SolutionDimensions = 2
    DO SurfNode = 1, NumOfPerpendNodes
      IF ( (SurfNode == 1) .OR. (SurfNode == NumOfPerpendNodes) ) THEN
        SurfNodeFac = 0.5d0
      ELSE
        SurfNodeFac = 1.0d0
      END IF
      s0(1,1) = s0(1,1) + SurfNodeFac*CMat(1)*Gamma2(SurfNode,1)
      s0(1,2) = s0(1,2) + SurfNodeFac*CMat(1)*Gamma2(SurfNode,2)
      s0(1,3) = s0(1,3) + SurfNodeFac*CMat(1)*Gamma2(SurfNode,3)
      s0(2,1) = s0(2,1) + SurfNodeFac*CMat(2)*Gamma2(rcmax+1-SurfNode,1)
      s0(2,2) = s0(2,2) + SurfNodeFac*CMat(2)*Gamma2(rcmax+1-SurfNode,2)
      s0(2,3) = s0(2,3) + SurfNodeFac*CMat(2)*Gamma2(rcmax+1-SurfNode,3)
    END DO
    s0(1,1) = s0(1,1) + REAL(NumOfPerpendNodes-1,r64)*DMat(1)
    s0(2,2) = s0(2,2) + REAL(NumOfPerpendNodes-1,r64)*DMat(2)
  END IF

  IF (NodeSource > 0) THEN
    s0(3,1) = Gamma2(NodeSource,1)
    s0(3,2) = Gamma2(NodeSource,2)
    s0(3,3) = Gamma2(NodeSource,3)
  END IF
  IF (NodeUserTemp > 0) THEN
    s0(4,1) = Gamma2(NodeUserTemp,1)
    s0(4,2) = Gamma2(NodeUserTemp,2)
    s0(4,3) = Gamma2(NodeUserTemp,3)
  END IF

          ! Check for and enforce symmetry in the cross term (Y)
  IF (ABS(s0(1,2)) /= ABS(s0(2,1))) THEN
    avg = (ABS(s0(1,2))+ABS(s0(2,1))) * 0.5D0
    s0(1,2) = avg*s0(1,2)/ABS(s0(1,2))
    s0(2,1) = avg*s0(2,1)/ABS(s0(2,1))
  END IF

          ! Compute S's and e's from 1 to n-1.  See equations (2.1.25) and
          ! (2.1.26) and Appendix C.
  inum = 1            ! Set history term counter
  CTFConvrg = .FALSE. ! Set the convergence logical to false

          ! The following DO WHILE loop calculates each successive set of time
          ! history terms until there are rcmax number of history terms or the
          ! latest flux history term is negligibly small compared to the first
          ! flux history term.
  DO WHILE ((.NOT.CTFConvrg) .AND. (inum < rcmax))   ! Begin CTF calculation loop ...

          ! Compute e(inum) based on Appendix C (Seem's dissertation). First,
          ! compute the new PhiR0 and its trace.

    trace = 0.0D0

    DO ir = 1,rcmax

      DO ic = 1,rcmax
        PhiR0(ir,ic) = 0.0D0
        DO is = 1,rcmax
          ! Make sure the next term won't cause an underflow.  If it will end up being
          ! so small as to go below TinyLimit, then ignore it since it won't add anything
          ! to PhiR0 anyway.
          IF (ABS(Rnew(is,ic)) > TinyLimit) THEN
            IF (ABS(AExp(ir,is)) > ABS(TinyLimit/Rnew(is,ic))) &
              PhiR0(ir,ic) = PhiR0(ir,ic)+AExp(ir,is)*Rnew(is,ic)
          END IF
        END DO
      END DO

      trace = trace + PhiR0(ir,ir)

    END DO

          ! Now calculate ej from the trace.  According to Appendix C:
          ! e(j) = -Trace[AExp*R(j-1)]/j

    e(inum) = -trace/REAL(inum,r64)

          ! Update Rold and compute Rnew.  Note:  PhiR0 = AExp*R(j-1) here.
          ! According to Appendix C:  R(j) = AExp*R(j-1) + e(j-1)

    DO ir = 1,rcmax
      DO ic = 1,rcmax
        Rold(ir,ic) = Rnew(ir,ic)
        Rnew(ir,ic) = PhiR0(ir,ic)
      END DO
      Rnew(ir,ir) = Rnew(ir,ir)+e(inum)
    END DO

          ! Compute S(inum) based on eq.(2.1.25) which states:
          ! S(j)  =  CMat*[R(j-1)*(Gamma1-Gamma2)+R(j)*Gamma2]
          !          + e(j)*DMat
    IF (SolutionDimensions == 1) THEN
      DO j = 1, 3
        DO is2 = 1,rcmax
          s(inum,1,j) = s(inum,1,j) + CMat(1)*( Rold(1,is2)*Gamma1(is2,j)  &
                                               +Rnew(1,is2)*Gamma2(is2,j) )
          s(inum,2,j) = s(inum,2,j) + CMat(2)*( Rold(rcmax,is2)*Gamma1(is2,j)  &
                                               +Rnew(rcmax,is2)*Gamma2(is2,j) )
          IF (NodeSource > 0) THEN
            s(inum,3,j) = s(inum,3,j) + ( Rold(NodeSource,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeSource,is2)*Gamma2(is2,j) )
          END IF
          IF (NodeUserTemp > 0) THEN
            s(inum,4,j) = s(inum,4,j) + ( Rold(NodeUserTemp,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeUserTemp,is2)*Gamma2(is2,j) )
          END IF

        END DO
        IF (j /= 3) s(inum,j,j) = s(inum,j,j) + e(inum)*DMat(j)
      END DO
    ELSE ! SolutionDimensions = 2
      DO j = 1, 3
        DO is2 = 1,rcmax
          DO SurfNode = 1, NumOfPerpendNodes
            IF ( (SurfNode == 1) .OR. (SurfNode == NumOfPerpendNodes) ) THEN
              SurfNodeFac = 0.5d0
            ELSE
              SurfNodeFac = 1.0d0
            END IF
            s(inum,1,j) = s(inum,1,j) + SurfNodeFac*CMat(1)*( Rold(SurfNode,is2)*Gamma1(is2,j)  &
                                                             +Rnew(SurfNode,is2)*Gamma2(is2,j) )
            s(inum,2,j) = s(inum,2,j) + SurfNodeFac*CMat(2)*( Rold(rcmax+1-SurfNode,is2)*Gamma1(is2,j)  &
                                                             +Rnew(rcmax+1-SurfNode,is2)*Gamma2(is2,j) )
          END DO
          IF (NodeSource > 0) THEN
            s(inum,3,j) = s(inum,3,j) + ( Rold(NodeSource,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeSource,is2)*Gamma2(is2,j) )
          END IF
          IF (NodeUserTemp > 0) THEN
            s(inum,4,j) = s(inum,4,j) + ( Rold(NodeUserTemp,is2)*Gamma1(is2,j)  &
                                         +Rnew(NodeUserTemp,is2)*Gamma2(is2,j) )
          END IF
        END DO
      END DO
      s(inum,1,1) = s(inum,1,1) + e(inum)*DMat(1)*REAL(NumOfPerpendNodes-1,r64)
      s(inum,2,2) = s(inum,2,2) + e(inum)*DMat(2)*REAL(NumOfPerpendNodes-1,r64)
    END IF

          ! Check for and enforce symmetry in the cross term (Y)
    IF (ABS(s(inum,1,2)) /= ABS(s(inum,2,1))) THEN
      avg = (ABS(s(inum,1,2))+ABS(s(inum,2,1))) * 0.5D0
      s(inum,1,2) = avg*s(inum,1,2)/ABS(s(inum,1,2))
      s(inum,2,1) = avg*s(inum,2,1)/ABS(s(inum,2,1))
    END IF

          ! Check for convergence of the CTFs.
    IF (e(1) == 0.0D0) THEN

      nrf = 1                   ! e(1) is zero, so there are no history terms.
      CTFConvrg = .TRUE.        ! CTF calculations have converged--set logical.

    ELSE
          ! e(1) is non-zero -- Calculate and compare the ratio of the flux
          ! terms to the convergence limit.
      rat = ABS(e(inum)/e(1))

      IF (rat < ConvrgLim) THEN

          ! If the ratio is less than the convergence limit, then any other
          ! terms would have a neglible impact on the CTF-based energy balances.
        nrf = inum
        CTFConvrg = .TRUE.    ! CTF calculations have converged--set logical.

      END IF
    END IF  ! ... end of convergence check block.

    inum = inum+1

  END DO    ! ... end of CTF calculation loop.
            ! Continue to the next coefficient if the solution has not converged
                    !
  IF (.NOT.CTFConvrg) THEN  ! Compute last e and S, if still unconverged.

          ! Compute e(inum) based on Appendix C (Seem's dissertation) or see
          ! equation above.  First compute the new PhiR0 and its trace.

    trace = 0.0D0

    DO ir = 1,rcmax
      DO is = 1,rcmax
        trace = trace+AExp(ir,is)*Rnew(is,ir)
      END DO
    END DO

    e(rcmax) = -trace/REAL(rcmax,r64)     ! Now calculate ej from the trace.

          ! Compute S(inum) based on eq.(2.1.25) which states:
          !   S(last) = CMat*R(last-1)*(Gamma1-Gamma2)+e(last)*DMat

    IF (SolutionDimensions == 1) THEN
      DO j = 1, 3
        DO is2 = 1,rcmax
          s(rcmax,1,j) = s(rcmax,1,j) + CMat(1)*Rnew(1,is2)*Gamma1(is2,j)
          s(rcmax,2,j) = s(rcmax,2,j) + CMat(2)*Rnew(rcmax,is2)*Gamma1(is2,j)
          IF (NodeSource > 0) THEN
            s(rcmax,3,j) = s(rcmax,3,j) + Rnew(NodeSource,is2)*Gamma1(is2,j)
          END IF
          IF (NodeUserTemp > 0) THEN
            s(rcmax,4,j) = s(rcmax,4,j) + Rnew(NodeUserTemp,is2)*Gamma1(is2,j)
          END IF
        END DO
      END DO
      s(rcmax,1,1) = s(rcmax,1,1)+e(rcmax)*DMat(1)
      s(rcmax,2,2) = s(rcmax,2,2)+e(rcmax)*DMat(2)
      nrf = rcmax
    ELSE ! SolutionDimensions = 2
      DO j = 1, 3
        DO is2 = 1,rcmax
          DO SurfNode = 1, NumOfPerpendNodes
            IF ( (SurfNode == 1) .OR. (SurfNode == NumOfPerpendNodes) ) THEN
              SurfNodeFac = 0.5d0
            ELSE
              SurfNodeFac = 1.0d0
            END IF
            s(rcmax,1,j) = s(rcmax,1,j) + SurfNodeFac*CMat(1)*Rnew(SurfNode,is2)*Gamma1(is2,j)
            s(rcmax,2,j) = s(rcmax,2,j) + SurfNodeFac*CMat(2)*Rnew(rcmax+1-SurfNode,is2)*Gamma1(is2,j)
          END DO
          IF (NodeSource > 0) THEN
            s(rcmax,3,j) = s(rcmax,3,j) + Rnew(NodeSource,is2)*Gamma1(is2,j)
          END IF
          IF (NodeUserTemp > 0) THEN
            s(rcmax,4,j) = s(rcmax,4,j) + Rnew(NodeUserTemp,is2)*Gamma1(is2,j)
          END IF
        END DO
      END DO
      s(rcmax,1,1) = s(rcmax,1,1) + e(rcmax)*DMat(1)*REAL(NumOfPerpendNodes-1,r64)
      s(rcmax,2,2) = s(rcmax,2,2) + e(rcmax)*DMat(2)*REAL(NumOfPerpendNodes-1,r64)
    END IF

          ! Check for and enforce symmetry in the cross term (Y)

    IF (ABS(s(rcmax,1,2)) /= ABS(s(rcmax,2,1))) THEN
      avg = (ABS(s(rcmax,1,2))+ABS(s(rcmax,2,1)))*0.5D0
      s(rcmax,1,2) = avg*s(rcmax,1,2)/ABS(s(rcmax,1,2))
      s(rcmax,2,1) = avg*s(rcmax,2,1)/ABS(s(rcmax,2,1))
    END IF

  END IF    ! ... end of IF block for calculation of last e and S.

  DEALLOCATE(PhiR0)
  DEALLOCATE(Rnew)
  DEALLOCATE(Rold)

  RETURN    ! The array e and the matrices s and s0 now contain the conduction
            ! transfer functions for this construction.

END SUBROUTINE CalculateCTFs

SUBROUTINE ReportCTFs(DoReportBecauseError)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Linda K. Lawrie
          !       DATE WRITTEN   July 1999
          !       MODIFIED       na
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This routine gives a detailed report to the user about
          ! the conduction transfer functions and other thermal data
          ! of each construction.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General, ONLY: ScanForReports

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  LOGICAL, INTENT(IN) :: DoReportBecauseError

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
!  CHARACTER(len=12),DIMENSION(6)  :: Roughness = (/'VeryRough   ','Rough       ','MediumRough ', &
!                                                   'MediumSmooth','Smooth      ','VerySmooth  '/)
  LOGICAL :: DoReport

  INTEGER ThisNum
  INTEGER Layer
  INTEGER I

  CALL ScanForReports('Constructions',DoReport,'Constructions')

  IF (DoReport .or. DoReportBecauseError) THEN
!
!                                      Write Descriptions
    Write(OutputFileInits,'(A)') '! <Construction CTF>,Construction Name,Index,#Layers,#CTFs,Time Step {hours},'// &
                             'ThermalConductance {w/m2-K},'// &
                             'OuterThermalAbsorptance,InnerThermalAbsorptance,OuterSolarAbsorptance,'// &
                             'InnerSolarAbsorptance,Roughness'
    Write(OutputFileInits,'(A)') '! <Material CTF Summary>,Material Name,Thickness {m},Conductivity {w/m-K},Density {kg/m3},'//  &
                             'Specific Heat {J/kg-K},ThermalResistance {m2-K/w}'
    Write(OutputFileInits,'(A)') '! <Material:Air>,Material Name,ThermalResistance {m2-K/w}'
    Write(OutputFileInits,'(A)') '! <CTF>,Time,Outside,Cross,Inside,Flux (except final one)'


    DO ThisNum=1,TotConstructs

      IF (Construct(ThisNum)%TypeIsWindow) CYCLE

      Write(OutputFileInits,700) TRIM(Construct(ThisNum)%Name),ThisNum, Construct(ThisNum)%TotLayers,  &
                                 Construct(ThisNum)%NumCTFTerms, &
                                 Construct(ThisNum)%CTFTimeStep,                                                          &
                                 Construct(ThisNum)%UValue,Construct(ThisNum)%OutsideAbsorpThermal,                         &
                                 Construct(ThisNum)%InsideAbsorpThermal,Construct(ThisNum)%OutsideAbsorpSolar,              &
                                 Construct(ThisNum)%InsideAbsorpSolar,  &
                                 TRIM(DisplayMaterialRoughness(Construct(ThisNum)%OutsideRoughness))
!
 700  FORMAT(' Construction CTF,',A,3(',',I4),',',F8.3,',',G15.4,4(',',F8.3),',',A)
 701  FORMAT(' Material CTF Summary,',A,',',F8.4,',',F14.3,',',F11.3,',',F13.3,',',G12.4)
 702  FORMAT(' Material:Air,',A,',',G12.4)
 703  FORMAT(' CTF,',I4,4(',',G20.8))
 704  FORMAT(' CTF,',I4,3(',',G20.8))
 705  FORMAT(' QTF,',I4,2(',',G20.8))
 706  FORMAT(' Source/Sink Loc Internal Temp QTF,',I4,3(',',G20.8))
 707  FORMAT(' User Loc Internal Temp QTF,',I4,3(',',G20.8))

      DO I=1,Construct(ThisNum)%TotLayers
        Layer=Construct(ThisNum)%LayerPoint(I)
        SELECT CASE (Material(Layer)%Group)
          CASE (Air)
            Write(OutputFileInits,702) TRIM(Material(Layer)%Name),Material(Layer)%Resistance
          CASE DEFAULT
            Write(OutputFileInits,701) TRIM(Material(Layer)%Name),Material(Layer)%Thickness,Material(Layer)%Conductivity,  &
                                       Material(Layer)%Density,Material(Layer)%SpecHeat,Material(Layer)%Resistance
        END SELECT
      ENDDO

      DO I=Construct(ThisNum)%NumCTFTerms,0,-1
        IF (I /= 0) THEN
          Write(OutputFileInits,703) I,Construct(ThisNum)%CTFOutside(I),Construct(ThisNum)%CTFCross(I),   &
                                     Construct(ThisNum)%CTFInside(I),Construct(ThisNum)%CTFFlux(I)
        ELSE
          Write(OutputFileInits,704) I,Construct(ThisNum)%CTFOutside(I),Construct(ThisNum)%CTFCross(I),   &
                                     Construct(ThisNum)%CTFInside(I)
        ENDIF
      ENDDO

      IF (Construct(ThisNum)%SourceSinkPresent) THEN
          ! QTFs...
        DO I=Construct(ThisNum)%NumCTFTerms,0,-1
          WRITE(OutputFileInits,705) I,Construct(ThisNum)%CTFSourceOut(I),Construct(ThisNum)%CTFSourceIn(I)
        END DO
          ! QTFs for source/sink location temperature calculation...
        DO I=Construct(ThisNum)%NumCTFTerms,0,-1
            WRITE(OutputFileInits,706) I,Construct(ThisNum)%CTFTSourceOut(I), &
                                         Construct(ThisNum)%CTFTSourceIn(I),  &
                                         Construct(ThisNum)%CTFTSourceQ(I)
        END DO
        IF (Construct(ThisNum)%TempAfterLayer /= 0) THEN
          ! QTFs for user specified interior temperature calculation...
          DO I=Construct(ThisNum)%NumCTFTerms,0,-1
            WRITE(OutputFileInits,707) I,Construct(ThisNum)%CTFTUserOut(I),  &
                                         Construct(ThisNum)%CTFTUserIn(I),   &
                                         Construct(ThisNum)%CTFTUserSource(I)
          END DO
        END IF
      END IF

    ENDDO

  ENDIF

  RETURN

END SUBROUTINE ReportCTFs

!     NOTICE
!
!     Copyright © 1996-2013 The Board of Trustees of the University of Illinois
!     and The Regents of the University of California through Ernest Orlando Lawrence
!     Berkeley National Laboratory.  All rights reserved.
!
!     Portions of the EnergyPlus software package have been developed and copyrighted
!     by other individuals, companies and institutions.  These portions have been
!     incorporated into the EnergyPlus software package under license.   For a complete
!     list of contributors, see "Notice" located in EnergyPlus.f90.
!
!     NOTICE: The U.S. Government is granted for itself and others acting on its
!     behalf a paid-up, nonexclusive, irrevocable, worldwide license in this data to
!     reproduce, prepare derivative works, and perform publicly and display publicly.
!     Beginning five (5) years after permission to assert copyright is granted,
!     subject to two possible five year renewals, the U.S. Government is granted for
!     itself and others acting on its behalf a paid-up, non-exclusive, irrevocable
!     worldwide license in this data to reproduce, prepare derivative works,
!     distribute copies to the public, perform publicly and display publicly, and to
!     permit others to do so.
!
!     TRADEMARKS: EnergyPlus is a trademark of the US Department of Energy.
!

END MODULE ConductionTransferFunctionCalc

AirflowNetworkBalanceManager.f90 AirflowNetworkSolver.f90 BaseboardRadiator.f90 BaseboardRadiatorElectric.f90 BaseboardRadiatorSteam.f90 BaseboardRadiatorWater.f90 BranchInputManager.f90 BranchNodeConnections.f90 ConductionTransferFunctionCalc.f90 CoolTower.f90 CostEstimateManager.f90 CurveManager.f90 CVFOnlyRoutines.f90 DataAirflowNetwork.f90 DataAirLoop.f90 DataAirSystems.f90 DataBranchAirLoopPlant.f90 DataBranchNodeConnections.f90 DataBSDFWindow.f90 DataComplexFenestration.f90 DataContaminantBalance.f90 DataConvergParams.f90 DataConversions.f90 DataCostEstimate.f90 DataDaylighting.f90 DataDaylightingDevices.f90 Datadefineequip.f90 DataDElight.f90 DataEnvironment.f90 DataEquivalentLayerWindow.f90 DataErrorTracking.f90 DataGenerators.f90 DataGlobalConstants.f90 DataGlobals.f90 DataHeatBalance.f90 DataHeatBalFanSys.f90 DataHeatBalSurface.f90 DataHVACControllers.f90 DataHVACGlobals.f90 DataInterfaces.f90 DataIPShortCuts.f90 DataLoopNode.f90 DataMoistureBalance.f90 DataMoistureBalanceEMPD.f90 DataOutputs.f90 DataPhotovoltaics.f90 DataPlant.f90 DataPlantPipingSystems.f90 DataPrecisionGlobals.f90 DataReportingFlags.f90 DataRoomAir.f90 DataRootFinder.f90 DataRuntimeLanguage.f90 DataShadowingCombinations.f90 DataSizing.f90 DataStringGlobals.f90 DataSurfaceColors.f90 DataSurfaceLists.f90 DataSurfaces.f90 DataSystemVariables.f90 DataTimings.f90 DataUCSDSharedData.f90 DataVectorTypes.f90 DataViewFactorInformation.f90 DataWater.f90 DataZoneControls.f90 DataZoneEnergyDemands.f90 DataZoneEquipment.f90 DaylightingDevices.f90 DaylightingManager.f90 DElightManagerF.f90 DElightManagerF_NO.f90 DemandManager.f90 DesiccantDehumidifiers.f90 DirectAir.f90 DisplayRoutines.f90 DXCoil.f90 EarthTube.f90 EconomicLifeCycleCost.f90 EconomicTariff.f90 EcoRoof.f90 ElectricPowerGenerators.f90 ElectricPowerManager.f90 EMSManager.f90 EnergyPlus.f90 ExteriorEnergyUseManager.f90 ExternalInterface_NO.f90 FanCoilUnits.f90 FaultsManager.f90 FluidProperties.f90 General.f90 GeneralRoutines.f90 GlobalNames.f90 HeatBalanceAirManager.f90 HeatBalanceConvectionCoeffs.f90 HeatBalanceHAMTManager.f90 HeatBalanceInternalHeatGains.f90 HeatBalanceIntRadExchange.f90 HeatBalanceManager.f90 HeatBalanceMovableInsulation.f90 HeatBalanceSurfaceManager.f90 HeatBalFiniteDifferenceManager.f90 HeatRecovery.f90 Humidifiers.f90 HVACControllers.f90 HVACCooledBeam.f90 HVACDualDuctSystem.f90 HVACDuct.f90 HVACDXSystem.f90 HVACEvapComponent.f90 HVACFanComponent.f90 HVACFurnace.f90 HVACHeatingCoils.f90 HVACHXAssistedCoolingCoil.f90 HVACInterfaceManager.f90 HVACManager.f90 HVACMixerComponent.f90 HVACMultiSpeedHeatPump.f90 HVACSingleDuctInduc.f90 HVACSingleDuctSystem.f90 HVACSplitterComponent.f90 HVACStandAloneERV.f90 HVACSteamCoilComponent.f90 HVACTranspiredCollector.f90 HVACUnitaryBypassVAV.f90 HVACUnitarySystem.f90 HVACVariableRefrigerantFlow.f90 HVACWaterCoilComponent.f90 HVACWatertoAir.f90 HVACWatertoAirMultiSpeedHP.f90 InputProcessor.f90 MatrixDataManager.f90 MixedAir.f90 MoistureBalanceEMPDManager.f90 NodeInputManager.f90 NonZoneEquipmentManager.f90 OutAirNodeManager.f90 OutdoorAirUnit.f90 OutputProcessor.f90 OutputReportPredefined.f90 OutputReports.f90 OutputReportTabular.f90 PackagedTerminalHeatPump.f90 PackagedThermalStorageCoil.f90 Photovoltaics.f90 PhotovoltaicThermalCollectors.f90 PlantAbsorptionChillers.f90 PlantBoilers.f90 PlantBoilersSteam.f90 PlantCentralGSHP.f90 PlantChillers.f90 PlantCondLoopOperation.f90 PlantCondLoopTowers.f90 PlantEIRChillers.f90 PlantEvapFluidCoolers.f90 PlantExhaustAbsorptionChiller.f90 PlantFluidCoolers.f90 PlantGasAbsorptionChiller.f90 PlantGroundHeatExchangers.f90 PlantHeatExchanger.f90 PlantIceThermalStorage.f90 PlantLoadProfile.f90 PlantLoopEquipment.f90 PlantLoopSolver.f90 PlantManager.f90 PlantOutsideEnergySources.f90 PlantPipeHeatTransfer.f90 PlantPipes.f90 PlantPipingSystemManager.f90 PlantPondGroundHeatExchanger.f90 PlantPressureSystem.f90 PlantPumps.f90 PlantSolarCollectors.f90 PlantSurfaceGroundHeatExchanger.f90 PlantUtilities.f90 PlantValves.f90 PlantWaterSources.f90 PlantWaterThermalTank.f90 PlantWatertoWaterGSHP.f90 PlantWaterUse.f90 PollutionAnalysisModule.f90 PoweredInductionUnits.f90 PsychRoutines.f90 Purchasedairmanager.f90 RadiantSystemHighTemp.f90 RadiantSystemLowTemp.f90 RefrigeratedCase.f90 ReportSizingManager.f90 ReturnAirPath.f90 RoomAirManager.f90 RoomAirModelCrossVent.f90 RoomAirModelDisplacementVent.f90 RoomAirModelMundt.f90 RoomAirModelUFAD.f90 RoomAirModelUserTempPattern.f90 RootFinder.f90 RuntimeLanguageProcessor.f90 ScheduleManager.f90 SetPointManager.f90 SimAirServingZones.f90 SimulationManager.f90 SizingManager.f90 SolarReflectionManager.f90 SolarShading.f90 SortAndStringUtilities.f90 sqlite3.c SQLiteCRoutines.c SQLiteFortranRoutines.f90 SQLiteFortranRoutines_NO.f90 StandardRatings.f90 SurfaceGeometry.f90 SystemAvailabilityManager.f90 SystemReports.f90 TarcogComplexFenestration.f90 ThermalChimney.f90 ThermalComfort.f90 UnitHeater.f90 UnitVentilator.f90 UserDefinedComponents.f90 UtilityRoutines.f90 VectorUtilities.f90 VentilatedSlab.f90 WaterManager.f90 WeatherManager.f90 WindowAC.f90 WindowComplexManager.f90 WindowEquivalentLayer.f90 WindowManager.f90 WindTurbine.f90 Zoneairloopequipmentmanager.f90 ZoneContaminantPredictorCorrector.f90 ZoneDehumidifier.f90 Zoneequipmentmanager.f90 ZonePlenumComponent.f90 ZoneTempPredictorCorrector.f90