HeatBalFiniteDifferenceManager.f90 Source File

This File Depends On

sourcefile~~heatbalfinitedifferencemanager.f90~~EfferentGraph sourcefile~heatbalfinitedifferencemanager.f90 HeatBalFiniteDifferenceManager.f90 sourcefile~datamoisturebalance.f90 DataMoistureBalance.f90 sourcefile~datamoisturebalance.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~general.f90 General.f90 sourcefile~general.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~psychroutines.f90 PsychRoutines.f90 sourcefile~general.f90->sourcefile~psychroutines.f90 sourcefile~dataheatbalance.f90 DataHeatBalance.f90 sourcefile~general.f90->sourcefile~dataheatbalance.f90 sourcefile~dataenvironment.f90 DataEnvironment.f90 sourcefile~general.f90->sourcefile~dataenvironment.f90 sourcefile~schedulemanager.f90 ScheduleManager.f90 sourcefile~general.f90->sourcefile~schedulemanager.f90 sourcefile~datasurfaces.f90 DataSurfaces.f90 sourcefile~datasurfaces.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~datasurfaces.f90->sourcefile~general.f90 sourcefile~heatbalancemovableinsulation.f90 HeatBalanceMovableInsulation.f90 sourcefile~datasurfaces.f90->sourcefile~heatbalancemovableinsulation.f90 sourcefile~datasurfaces.f90->sourcefile~dataheatbalance.f90 sourcefile~dataheatbalsurface.f90 DataHeatBalSurface.f90 sourcefile~dataheatbalsurface.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~psychroutines.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataprecisionglobals.f90 DataPrecisionGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datamoisturebalance.f90 sourcefile~dataprecisionglobals.f90->sourcefile~general.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasurfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataheatbalsurface.f90 sourcefile~dataprecisionglobals.f90->sourcefile~psychroutines.f90 sourcefile~dataglobals.f90 DataGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataglobals.f90 sourcefile~dataipshortcuts.f90 DataIPShortCuts.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataipshortcuts.f90 sourcefile~dataairflownetwork.f90 DataAirflowNetwork.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataairflownetwork.f90 sourcefile~dataheatbalfansys.f90 DataHeatBalFanSys.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataheatbalfansys.f90 sourcefile~dataprecisionglobals.f90->sourcefile~heatbalancemovableinsulation.f90 sourcefile~datainterfaces.f90 DataInterfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datainterfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataheatbalance.f90 sourcefile~inputprocessor.f90 InputProcessor.f90 sourcefile~dataprecisionglobals.f90->sourcefile~inputprocessor.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataenvironment.f90 sourcefile~datahvacglobals.f90 DataHVACGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datahvacglobals.f90 sourcefile~dataruntimelanguage.f90 DataRuntimeLanguage.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataruntimelanguage.f90 sourcefile~datavectortypes.f90 DataVectorTypes.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datavectortypes.f90 sourcefile~databsdfwindow.f90 DataBSDFWindow.f90 sourcefile~dataprecisionglobals.f90->sourcefile~databsdfwindow.f90 sourcefile~dataprecisionglobals.f90->sourcefile~schedulemanager.f90 sourcefile~datasystemvariables.f90 DataSystemVariables.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasystemvariables.f90 sourcefile~datacomplexfenestration.f90 DataComplexFenestration.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datacomplexfenestration.f90 sourcefile~dataequivalentlayerwindow.f90 DataEquivalentLayerWindow.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataequivalentlayerwindow.f90 sourcefile~datasizing.f90 DataSizing.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasizing.f90 sourcefile~dataglobals.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataglobals.f90->sourcefile~general.f90 sourcefile~dataglobals.f90->sourcefile~datasurfaces.f90 sourcefile~dataglobals.f90->sourcefile~psychroutines.f90 sourcefile~dataglobals.f90->sourcefile~dataipshortcuts.f90 sourcefile~dataglobals.f90->sourcefile~dataairflownetwork.f90 sourcefile~dataglobals.f90->sourcefile~dataheatbalance.f90 sourcefile~dataglobals.f90->sourcefile~inputprocessor.f90 sourcefile~dataglobals.f90->sourcefile~dataenvironment.f90 sourcefile~dataglobals.f90->sourcefile~datahvacglobals.f90 sourcefile~dataglobals.f90->sourcefile~dataruntimelanguage.f90 sourcefile~dataglobals.f90->sourcefile~databsdfwindow.f90 sourcefile~dataglobals.f90->sourcefile~schedulemanager.f90 sourcefile~dataglobals.f90->sourcefile~datacomplexfenestration.f90 sourcefile~dataglobals.f90->sourcefile~dataequivalentlayerwindow.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~dataipshortcuts.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataipshortcuts.f90->sourcefile~general.f90 sourcefile~dataipshortcuts.f90->sourcefile~inputprocessor.f90 sourcefile~dataipshortcuts.f90->sourcefile~schedulemanager.f90 sourcefile~dataairflownetwork.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataheatbalfansys.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~heatbalancemovableinsulation.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~datainterfaces.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~datainterfaces.f90->sourcefile~general.f90 sourcefile~datainterfaces.f90->sourcefile~psychroutines.f90 sourcefile~datainterfaces.f90->sourcefile~heatbalancemovableinsulation.f90 sourcefile~datainterfaces.f90->sourcefile~dataheatbalance.f90 sourcefile~datainterfaces.f90->sourcefile~inputprocessor.f90 sourcefile~datainterfaces.f90->sourcefile~dataenvironment.f90 sourcefile~datainterfaces.f90->sourcefile~dataruntimelanguage.f90 sourcefile~datainterfaces.f90->sourcefile~schedulemanager.f90 sourcefile~dataheatbalance.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataheatbalance.f90->sourcefile~heatbalancemovableinsulation.f90 sourcefile~inputprocessor.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~inputprocessor.f90->sourcefile~general.f90 sourcefile~inputprocessor.f90->sourcefile~dataheatbalance.f90 sourcefile~inputprocessor.f90->sourcefile~schedulemanager.f90 sourcefile~dataenvironment.f90->sourcefile~heatbalfinitedifferencemanager.f90 sourcefile~dataenvironment.f90->sourcefile~psychroutines.f90 sourcefile~dataenvironment.f90->sourcefile~dataheatbalance.f90 sourcefile~dataenvironment.f90->sourcefile~schedulemanager.f90 sourcefile~datahvacglobals.f90->sourcefile~general.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~schedulemanager.f90 sourcefile~datastringglobals.f90->sourcefile~datasystemvariables.f90 sourcefile~datavectortypes.f90->sourcefile~datasurfaces.f90 sourcefile~datavectortypes.f90->sourcefile~dataheatbalance.f90 sourcefile~datavectortypes.f90->sourcefile~databsdfwindow.f90 sourcefile~databsdfwindow.f90->sourcefile~datasurfaces.f90 sourcefile~databsdfwindow.f90->sourcefile~dataheatbalance.f90 sourcefile~schedulemanager.f90->sourcefile~heatbalancemovableinsulation.f90 sourcefile~datasystemvariables.f90->sourcefile~inputprocessor.f90 sourcefile~datasystemvariables.f90->sourcefile~schedulemanager.f90 sourcefile~datacomplexfenestration.f90->sourcefile~dataheatbalance.f90 sourcefile~dataequivalentlayerwindow.f90->sourcefile~dataheatbalance.f90 sourcefile~sortandstringutilities.f90->sourcefile~inputprocessor.f90 sourcefile~dataoutputs.f90->sourcefile~inputprocessor.f90 sourcefile~datasizing.f90->sourcefile~inputprocessor.f90
Help

Files Dependent On This One

sourcefile~~heatbalfinitedifferencemanager.f90~~AfferentGraph sourcefile~heatbalfinitedifferencemanager.f90 HeatBalFiniteDifferenceManager.f90 sourcefile~heatbalancesurfacemanager.f90 HeatBalanceSurfaceManager.f90 sourcefile~heatbalfinitedifferencemanager.f90->sourcefile~heatbalancesurfacemanager.f90 sourcefile~heatbalancemanager.f90 HeatBalanceManager.f90 sourcefile~heatbalancesurfacemanager.f90->sourcefile~heatbalancemanager.f90 sourcefile~simulationmanager.f90 SimulationManager.f90 sourcefile~heatbalancesurfacemanager.f90->sourcefile~simulationmanager.f90 sourcefile~heatbalancemanager.f90->sourcefile~simulationmanager.f90 sourcefile~sizingmanager.f90 SizingManager.f90 sourcefile~heatbalancemanager.f90->sourcefile~sizingmanager.f90 sourcefile~energyplus.f90 EnergyPlus.f90 sourcefile~simulationmanager.f90->sourcefile~energyplus.f90 sourcefile~utilityroutines.f90 UtilityRoutines.f90 sourcefile~simulationmanager.f90->sourcefile~utilityroutines.f90 sourcefile~sizingmanager.f90->sourcefile~simulationmanager.f90
Help


Source Code

MODULE HeatBalFiniteDiffManager

          ! Module containing the heat balance simulation routines

          ! MODULE INFORMATION:
          !       AUTHOR         Richard J. Liesen
          !       DATE WRITTEN   October 2003
          !       RE-ENGINEERED  Curtis Pedersen, 2006, Changed to Implicit FD calc for conduction.
          !                      and included enthalpy formulations for phase change materials
          ! PURPOSE OF THIS MODULE:
          ! To encapsulate the data and algorithms required to
          ! manage the fiite difference heat balance simulation on the building.

          ! METHODOLOGY EMPLOYED:
          !

          ! REFERENCES:
          ! The MFD moisture balance method
          !  C. O. Pedersen, Enthalpy Formulation of conduction heat transfer problems
          !    involving latent heat, Simulation, Vol 18, No. 2, February 1972

          ! OTHER NOTES:
          !
          !

          ! USE STATEMENTS:
          ! Use statements for data only modules
USE DataPrecisionGlobals
USE DataGlobals, Only:MaxNameLength,KelvinConv,outputfiledebug, NumOfTimeStepInHour,DisplayExtraWarnings,   &
                      TimeStepZone,Hourofday,TimeStep,OutputFileInits,BeginEnvrnFlag,DayofSim,WarmupFlag,SecInHour
USE DataInterfaces
USE DataMoistureBalance
USE DataHeatBalance, Only: Material, TotMaterials,MaxLayersInConstruct, QRadThermInAbs,Construct,      &
                           TotConstructs,UseCondFD,RegularMaterial,Air,Zone
USE DataHeatBalSurface, Only: NetLWRadToSurf, QRadSWOutAbs, QRadSWInAbs,QRadSWOutMvIns, TempSource, &
                         OpaqSurfInsFaceConductionFlux,OpaqSurfOutsideFaceConductionFlux,  &
                         OpaqSurfOutsideFaceConduction,QsrcHist,OpaqSurfInsFaceConduction,  &
                         QdotRadOutRepPerArea,MinSurfaceTempLimit,MaxSurfaceTempLimit,QdotRadNetSurfInRep,QRadNetSurfInReport
Use DataSurfaces, Only: Surface, TotSurfaces,Ground,SurfaceClass_Window, HeatTransferModel_CondFD
Use DataHeatBalFanSys, Only: Mat, ZoneAirHumRat, QHTRadSysSurf, QHWBaseboardSurf, QSteamBaseboardSurf, QElecBaseboardSurf
Use DataEnvironment, Only: SkyTemp, IsRain
USE Psychrometrics , Only:PsyRhFnTdbRhovLBnd0C,PsyWFnTdbRhPb,PsyHgAirFnWTdb
  ! Fan system Source/Sink heat value, and source/sink location temp from CondFD
USE DataHeatBAlFanSys,  ONLY: QRadSysSource,TCondFDSourceNode,QPVSysSource
USE HeatBalanceMovableInsulation,  ONLY : EvalOutsideMovableInsulation

IMPLICIT NONE   ! Enforce explicit typing of all variables

PRIVATE   !Make everything private except what is specifically Public

          ! MODULE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: Lambda=2500000.0d0
REAL(r64), PARAMETER :: smalldiff=1.d-8  ! Used in places where "equality" tests should not be used.

INTEGER, PARAMETER :: CrankNicholsonSecondOrder = 1 ! original CondFD scheme.  semi implicit, second order in time
INTEGER, PARAMETER :: FullyImplicitFirstOrder  = 2 ! fully implicit scheme, first order in time.
CHARACTER(len=*), PARAMETER, DIMENSION(2) :: cCondFDSchemeType=  &
   (/'CrankNicholsonSecondOrder',  &
     'FullyImplicitFirstOrder  '/)

REAL(r64),PARAMETER :: TempInitValue = 23.0d0    ! Initialization value for Temperature
REAL(r64),PARAMETER :: RhovInitValue = 0.0115d0  ! Initialization value for Rhov
REAL(r64),PARAMETER :: EnthInitValue = 100.0d0   ! Initialization value for Enthalpy

          ! DERIVED TYPE DEFINITIONS:
TYPE, PUBLIC ::  ConstructionDataFD
  CHARACTER(len=MaxNameLength), ALLOCATABLE, DIMENSION(:) :: Name ! Name of construction
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: DelX                 !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TempStability
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: MoistStability

  Integer, ALLOCATABLE, DIMENSION(:) :: NodeNumPoint  !
!  INTEGER, ALLOCATABLE, DIMENSION(:) :: InterfaceNodeNums   ! Layer interfaces occur at these nodes

  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: Thickness
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: NodeXlocation ! sized to TotNode, contains X distance in m from outside face
  INTEGER :: TotNodes=0

  Integer :: DeltaTime=0

END TYPE ConstructionDataFD


TYPE, PUBLIC :: SurfaceDataFD
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: T             !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TOld          !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TT
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: Rhov          !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: RhovOld       !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: RhoT
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TD            !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TDT            !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TDTLast
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TDOld         !
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: TDreport
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: RH
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: RHreport
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: EnthOld        ! Current node enthalpy
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: EnthNew        ! Node enthalpy at new time
  REAL(r64),    ALLOCATABLE, DIMENSION(:) :: EnthLast
  INTEGER                                 :: GSloopCounter = 0 ! count of inner loop iterations
  INTEGER                                 :: GSloopErrorCount = 0 ! recurring error counter
  REAL(r64)                               :: MaxNodeDelTemp = 0.0d0 ! largest change in node temps after calc

END TYPE SurfaceDataFD

TYPE MaterialDataFD
  REAL(r64) :: tk1          =0.0d0  ! Temperature coefficient for thermal conductivity
  INTEGER   :: numTempEnth  = 0   ! number of Temperature/Enthalpy pairs
  INTEGER   :: numTempCond  = 0   ! number of Temperature/Conductivity pairs
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)  :: TempEnth  !  Temperature enthalpy Function Pairs,
                                                       !  TempEnth(1,1)= first Temp, TempEnth(1,2) = First Enthalpy,
                                                       !  TempEnth(2,1) = secomd Temp, etc.
  REAL(r64), ALLOCATABLE, DIMENSION(:,:)  :: TempCond  !  Temperature thermal conductivity Function Pairs,
                                                       !  TempCond(1,1)= first Temp, Tempcond(1,2) = First conductivity,
                                                       !  TempEnth(2,1) = secomd Temp, etc.
END TYPE
          ! MODULE VARIABLE DECLARATIONS:
TYPE (ConstructionDataFD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: ConstructFD
TYPE (SurfaceDataFD), PUBLIC,   ALLOCATABLE, DIMENSION(:) :: SurfaceFD
TYPE (MaterialDataFD), ALLOCATABLE, DIMENSION(:) :: MaterialFD

!REAL(r64) :: TFDout   =0.0d0
!REAL(r64) :: TFDin    =0.0d0
!REAL(r64) :: rhovFDout=0.0d0
!REAL(r64) :: rhovFDin =0.0d0
!REAL(r64) :: TDryout  =0.0d0
!REAL(r64) :: Tdryin   =0.0d0
!REAL(r64) :: RHOut    =0.0d0
!REAL(r64) :: RHIn     =0.0d0
REAL(r64), ALLOCATABLE, DIMENSION(:) :: SigmaR      !  Total Resistance of construction layers
REAL(r64), ALLOCATABLE, DIMENSION(:) :: SigmaC      !  Total Capacitance of construction layers

!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: WSurfIn         !Humidity Ratio of the inside surface for reporting
!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QMassInFlux     !MassFlux on Surface for reporting
!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QMassOutFlux    !MassFlux on Surface for reporting
REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QHeatInFlux     !HeatFlux on Surface for reporting
REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QHeatOutFlux    !HeatFlux on Surface for reporting
!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QFluxZoneToInSurf !sum of Heat flows at the surface to air interface,
!                                 ! zone-side boundary conditions W/m2 before CR 8280 was not reported, but was calculated.
!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QFluxOutsideToOutSurf !sum of Heat flows at the surface to air interface, Out-side boundary conditions W/m2
!                                                           ! before CR 8280 was
!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QFluxInArrivSurfCond !conduction between surface node and first node into the surface (sensible)
!                                                           ! before CR 8280 was -- Qdryin    !HeatFlux on Surface for reporting for Sensible only
!REAL(r64), ALLOCATABLE, DIMENSION(:)   :: QFluxOutArrivSurfCond  !HeatFlux on Surface for reporting for Sensible only
!                                                                 ! before CR 8280 -- Qdryout         !HeatFlux on Surface for reporting for Sensible only


INTEGER   :: CondFDSchemeType = FullyImplicitFirstOrder  ! solution scheme for CondFD - default
REAL(r64) :: SpaceDescritConstant = 3.d0 ! spatial descritization constant,
REAL(r64) :: MinTempLimit = -100.d0 ! lower limit check, degree C
REAL(r64) :: MaxTempLimit =  100.d0 ! upper limit check, degree C
!feb2012 INTEGER   :: MaxGSiter = 200  ! maximum number of Gauss Seidel iterations
INTEGER   :: MaxGSiter = 30  ! maximum number of Gauss Seidel iterations
REAL(r64) :: fracTimeStepZone_Hour=0.0d0
LOGICAL   :: GetHBFiniteDiffInputFlag=.true.
INTEGER :: WarmupSurfTemp=0
          ! Subroutine Specifications for the Heat Balance Module
          ! Driver Routines
PUBLIC  ManageHeatBalFiniteDiff

          ! Initialization routines for module
PUBLIC  InitHeatBalFiniteDiff
PRIVATE GetCondFDInput

          ! Algorithms for the module
PRIVATE CalcHeatBalFiniteDiff
PRIVATE ExteriorBCEqns
PRIVATE InteriorBCEqns
PRIVATE terpld
PRIVATE InteriorNodeEqns
PRIVATE IntInterfaceNodeEqns


          ! Reporting routines for module
PRIVATE ReportFiniteDiffInits

          ! Update Data Routine
PUBLIC  UpdateMoistureBalanceFD


CONTAINS

! MODULE SUBROUTINES:
!*************************************************************************
SUBROUTINE ManageHeatBalFiniteDiff(SurfNum,TempSurfInTmp,TempSurfOutTmp)
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard Liesen
          !       DATE WRITTEN   May 2000
          !       MODIFIED       na
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine manages the moisture balance method.  It is called
          ! from the HeatBalanceManager at the time step level.
          ! This driver manages the calls to all of
          ! the other drivers and simulation algorithms.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:

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

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
  Integer, Intent(In) :: SurfNum
  REAL(r64), Intent(InOut) :: TempSurfInTmp       !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64), Intent(InOut) :: TempSurfOutTmp      !Outside Surface Temperature of each Heat Transfer Surface

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:


          ! FLOW:

  ! Get the moisture balance input at the beginning of the simulation only
  IF (GetHBFiniteDiffInputFlag) THEN
      ! Obtains conduction FD related parameters from input file
    CALL GetCondFDInput
    GetHBFiniteDiffInputFlag=.false.
  ENDIF

  ! Condition is taken care of by calling routine:
  ! IF (Surface(SurfNum)%HeatTransSurf .and. Surface(SurfNum)%Class /= SurfaceClass_Window)

  ! Solve the zone heat & moisture balance using a finite difference solution
  CALL CalcHeatBalFiniteDiff(SurfNum,TempSurfInTmp,TempSurfOutTmp)

  RETURN

END SUBROUTINE ManageHeatBalFiniteDiff


! Get Input Section of the Module
!******************************************************************************
SUBROUTINE GetCondFDInput
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Curtis Pedersen
          !       DATE WRITTEN   July 2006
          !       MODIFIED       Brent Griffith Mar 2011, user settings
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine is the main driver for initializations for the variable property CondFD part of the
          ! MFD algorithm

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE DataIPShortCuts
  USE InputProcessor, ONLY: GetNumObjectsFound,GetObjectItem,FindItemInList
  USE DataHeatBalance, ONLY: MaxAllowedDelTempCondFD, CondFDRelaxFactor, CondFDRelaxFactorInput
  USE General, ONLY: RoundSigDigits

  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:
  INTEGER :: IOStat           ! IO Status when calling get input subroutine
  CHARACTER(len=MaxNameLength),DIMENSION(3) &
          :: MaterialNames ! Number of Material Alpha names defined
  CHARACTER(len=MaxNameLength),DIMENSION(3)  :: ConstructionName ! Name of Construction with CondFDsimplified
  INTEGER :: MaterNum         ! Counter to keep track of the material number
  INTEGER :: MaterialNumAlpha ! Number of material alpha names being passed
  INTEGER :: MaterialNumProp  ! Number of material properties being passed
  REAL(r64), DIMENSION(40) :: MaterialProps !Temporary array to transfer material properties
  LOGICAL :: ErrorsFound = .false. ! If errors detected in input
!  INTEGER :: CondFDMat                ! Number of variable property CondFD materials in input
  INTEGER :: ConstructNumber     ! Cconstruction with CondHBsimple to be overridden with CondHBdetailed

  INTEGER :: NumConstructionAlpha
  Integer :: Loop
  INTEGER :: NumAlphas
  INTEGER :: NumNumbers
  INTEGER :: propNum
  INTEGER :: pcount
  INTEGER :: pcMat
  INTEGER :: vcMat
  INTEGER :: inegptr
  LOGICAL :: nonInc

  ! user settings for numerical parameters
  cCurrentModuleObject = 'HeatBalanceSettings:ConductionFiniteDifference'

  IF (GetNumObjectsFound(cCurrentModuleObject) > 0) THEN
    CALL GetObjectItem(cCurrentModuleObject,1,cAlphaArgs,NumAlphas, &
                       rNumericArgs,NumNumbers,IOSTAT,  &
                   AlphaBlank=lAlphaFieldBlanks,NumBlank=lNumericFieldBlanks,  &
                   AlphaFieldnames=cAlphaFieldNames,NumericFieldNames=cNumericFieldNames)

    IF (.NOT. lAlphaFieldBlanks(1)) THEN

      SELECT CASE (cAlphaArgs(1))

      CASE ('CRANKNICHOLSONSECONDORDER')
        CondFDSchemeType = CrankNicholsonSecondOrder
      CASE ('FULLYIMPLICITFIRSTORDER')
        CondFDSchemeType = FullyImplicitFirstOrder
      CASE DEFAULT
        CALL ShowSevereError(TRIM(cCurrentModuleObject)//': invalid '//TRIM(cAlphaFieldNames(1))//  &
         ' entered='//TRIM(cAlphaArgs(1))//', must match CrankNicholsonSecondOrder or FullyImplicitFirstOrder.')
        ErrorsFound=.true.
      END SELECT

    ENDIF

    IF (.NOT. lNumericFieldBlanks(1)) THEN
      SpaceDescritConstant = rNumericArgs(1)
    ENDIF
    IF (.NOT. lNumericFieldBlanks(2)) THEN
      CondFDRelaxFactorInput = rNumericArgs(2)
      CondFDRelaxFactor      = CondFDRelaxFactorInput
    ENDIF
    IF (.NOT. lNumericFieldBlanks(3)) THEN
      MaxAllowedDelTempCondFD = rNumericArgs(3)
    ENDIF

  ENDIF ! settings object

  pcMat=GetNumObjectsFound('MaterialProperty:PhaseChange')
  vcMat=GetNumObjectsFound('MaterialProperty:VariableThermalConductivity')

  ALLOCATE(MaterialFD(TotMaterials))

  ! Load the additional CondFD Material properties
  cCurrentModuleObject='MaterialProperty:PhaseChange'    ! Phase Change Information First

  IF ( pcMat .NE. 0 ) Then                      !  Get Phase Change info
!    CondFDVariableProperties = .TRUE.
    DO Loop=1,pcMat

    !Call Input Get routine to retrieve material data
      CALL GetObjectItem(cCurrentModuleObject,Loop,MaterialNames,MaterialNumAlpha, &
                       MaterialProps,MaterialNumProp,IOSTAT,  &
                   AlphaBlank=lAlphaFieldBlanks,NumBlank=lNumericFieldBlanks,  &
                   AlphaFieldnames=cAlphaFieldNames,NumericFieldNames=cNumericFieldNames)


    !Load the material derived type from the input data.
      MaterNum = FindItemInList(MaterialNames(1),Material%Name,TotMaterials)
      IF (MaterNum == 0) THEN
        CALL ShowSevereError(TRIM(cCurrentModuleObject)//': invalid '//TRIM(cAlphaFieldNames(1))//  &
           ' entered='//TRIM(MaterialNames(1))//', must match to a valid Material name.')
        ErrorsFound=.true.
        Cycle
      ENDIF

      IF (Material(MaterNum)%Group /= RegularMaterial) THEN
        CALL ShowSevereError(TRIM(cCurrentModuleObject)//  &
           ': Reference Material is not appropriate type for CondFD properties, material='//  &
           TRIM(Material(MaterNum)%Name)//', must have regular properties (L,Cp,K,D)')
        ErrorsFound=.true.
      ENDIF


    ! Once the material derived type number is found then load the additional CondFD variable material properties
    !   Some or all may be zero (default).  They will be checked when calculating node temperatures
      MaterialFD(MaterNum)%tk1          = MaterialProps(1)
      MaterialFD(MaterNum)%numTempEnth  = (MaterialNumProp-1)/2
      IF (MaterialFD(MaterNum)%numTempEnth*2 /= (MaterialNumProp-1)) THEN
        CALL ShowSevereError('GetCondFDInput: '//trim(cCurrentModuleObject)//'="'//trim(MaterialNames(1))//  &
           '", mismatched pairs')
        CALL ShowContinueError('...expected '//trim(RoundSigDigits(MaterialFD(MaterNum)%numTempEnth))//  &
           ' pairs, but only entered '//trim(RoundSigDigits(MaterialNumProp-1))//' numbers.')
        ErrorsFound=.true.
      ENDIF
      ALLOCATE(MaterialFD(MaterNum)%TempEnth(MaterialFD(MaterNum)%numTempEnth,2))
      MaterialFD(MaterNum)%TempEnth     =0.0d0
      propNum=2
      ! Temperature first
      DO pcount=1,MaterialFD(MaterNum)%numTempEnth
        MaterialFD(MaterNum)%TempEnth(pcount,1)   = MaterialProps(propNum)
        propNum=propNum+2
      ENDDO
      propNum=3
      ! Then Enthalpy
      DO pcount=1,MaterialFD(MaterNum)%numTempEnth
        MaterialFD(MaterNum)%TempEnth(pcount,2)   = MaterialProps(propNum)
        propNum=propNum+2
      ENDDO
      nonInc=.false.
      inegptr=0
      DO pcount=1,MaterialFD(MaterNum)%numTempEnth-1
        IF (MaterialFD(MaterNum)%TempEnth(pcount,1) < MaterialFD(MaterNum)%TempEnth(pcount+1,1)) CYCLE
        nonInc=.true.
        inegptr=pcount+1
        EXIT
      ENDDO
      IF (nonInc) THEN
        CALL ShowSevereError('GetCondFDInput: '//trim(cCurrentModuleObject)//'="'//trim(MaterialNames(1))//  &
           '", non increasing Temperatures. Temperatures must be strictly increasing.')
        CALL ShowContinueError('...occurs first at item=['//trim(RoundSigDigits(inegptr))//'], value=['//  &
           trim(RoundSigDigits(MaterialFD(MaterNum)%TempEnth(inegptr,1),2))//'].')
        ErrorsFound=.true.
      ENDIF
      nonInc=.false.
      inegptr=0
      DO pcount=1,MaterialFD(MaterNum)%numTempEnth-1
        IF (MaterialFD(MaterNum)%TempEnth(pcount,2) <= MaterialFD(MaterNum)%TempEnth(pcount+1,2)) CYCLE
        nonInc=.true.
        inegptr=pcount+1
        EXIT
      ENDDO
      IF (nonInc) THEN
        CALL ShowSevereError('GetCondFDInput: '//trim(cCurrentModuleObject)//'="'//trim(MaterialNames(1))//  &
           '", non increasing Enthalpy.')
        CALL ShowContinueError('...occurs first at item=['//trim(RoundSigDigits(inegptr))//'], value=['//  &
           trim(RoundSigDigits(MaterialFD(MaterNum)%TempEnth(inegptr,2),2))//'].')
        CALL ShowContinueError('...These values may be Cp (Specific Heat) rather than Enthalpy.  Please correct.')
        ErrorsFound=.true.
      ENDIF
    ENDDO
  END IF
!   Get CondFD Variable Thermal Conductivity Input

  cCurrentModuleObject='MaterialProperty:VariableThermalConductivity'    ! Variable Thermal Conductivity Info next
  IF ( vcMat .NE. 0 ) Then   !  variable k info
!    CondFDVariableProperties = .TRUE.
    DO Loop=1,vcMat

      !Call Input Get routine to retrieve material data
      CALL GetObjectItem(cCurrentModuleObject,Loop,MaterialNames,MaterialNumAlpha, &
                       MaterialProps,MaterialNumProp,IOSTAT,  &
                   AlphaBlank=lAlphaFieldBlanks,NumBlank=lNumericFieldBlanks,  &
                   AlphaFieldnames=cAlphaFieldNames,NumericFieldNames=cNumericFieldNames)


      !Load the material derived type from the input data.
      MaterNum = FindItemInList(MaterialNames(1),Material%Name,TotMaterials)
      IF (MaterNum == 0) THEN
        CALL ShowSevereError(TRIM(cCurrentModuleObject)//': invalid '//TRIM(cAlphaFieldNames(1))//  &
         ' entered='//TRIM(MaterialNames(1))//', must match to a valid Material name.')
        ErrorsFound=.true.
        Cycle
      ENDIF

      IF (Material(MaterNum)%Group /= RegularMaterial) THEN
        CALL ShowSevereError(TRIM(cCurrentModuleObject)//  &
           ': Reference Material is not appropriate type for CondFD properties, material='//  &
           TRIM(Material(MaterNum)%Name)//', must have regular properties (L,Cp,K,D)')
        ErrorsFound=.true.
      ENDIF


    ! Once the material derived type number is found then load the additional CondFD variable material properties
    !   Some or all may be zero (default).  They will be checked when calculating node temperatures
      MaterialFD(MaterNum)%numTempCond  = MaterialNumProp/2
      IF (MaterialFD(MaterNum)%numTempCond*2 /= MaterialNumProp) THEN
        CALL ShowSevereError('GetCondFDInput: '//trim(cCurrentModuleObject)//'="'//trim(MaterialNames(1))//  &
           '", mismatched pairs')
        CALL ShowContinueError('...expected '//trim(RoundSigDigits(MaterialFD(MaterNum)%numTempCond))//  &
           ' pairs, but only entered '//trim(RoundSigDigits(MaterialNumProp))//' numbers.')
        ErrorsFound=.true.
      ENDIF
      ALLOCATE(MaterialFD(MaterNum)%TempCond(MaterialFD(MaterNum)%numTempCond,2))
      MaterialFD(MaterNum)%TempCond    =0.0d0
      propNum=1
      ! Temperature first
      DO pcount=1,MaterialFD(MaterNum)%numTempCond
        MaterialFD(MaterNum)%TempCond(pcount,1)   = MaterialProps(propNum)
        propNum=propNum+2
      ENDDO
      propNum=2
      ! Then Conductivity
      DO pcount=1,MaterialFD(MaterNum)%numTempCond
        MaterialFD(MaterNum)%TempCond(pcount,2)   = MaterialProps(propNum)
        propNum=propNum+2
      ENDDO
      nonInc=.false.
      inegptr=0
      DO pcount=1,MaterialFD(MaterNum)%numTempCond-1
        IF (MaterialFD(MaterNum)%TempCond(pcount,1) < MaterialFD(MaterNum)%TempCond(pcount+1,1)) CYCLE
        nonInc=.true.
        inegptr=pcount+1
        EXIT
      ENDDO
      IF (nonInc) THEN
        CALL ShowSevereError('GetCondFDInput: '//trim(cCurrentModuleObject)//'="'//trim(MaterialNames(1))//  &
           '", non increasing Temperatures. Temperatures must be strictly increasing.')
        CALL ShowContinueError('...occurs first at item=['//trim(RoundSigDigits(inegptr))//'], value=['//  &
           trim(RoundSigDigits(MaterialFD(MaterNum)%TempCond(inegptr,1),2))//'].')
        ErrorsFound=.true.
      ENDIF
    ENDDO
  END IF

  DO MaterNum=1,TotMaterials
    IF (MaterialFD(MaterNum)%numTempEnth == 0) THEN
      MaterialFD(MaterNum)%numTempEnth=3
      ALLOCATE(MaterialFD(MaterNum)%TempEnth(MaterialFD(MaterNum)%numTempEnth,2))
      MaterialFD(MaterNum)%TempEnth=-100.0d0
    ENDIF
    IF (MaterialFD(MaterNum)%numTempCond == 0) THEN
      MaterialFD(MaterNum)%numTempCond=3
      ALLOCATE(MaterialFD(MaterNum)%TempCOnd(MaterialFD(MaterNum)%numTempCond,2))
      MaterialFD(MaterNum)%TempCond=-100.0d0
    ENDIF
  ENDDO

  IF (ErrorsFound) THEN
    CALL ShowFatalError('GetCondFDInput: Errors found getting ConductionFiniteDifference properties. Program terminates.')
  ENDIF

  CALL InitialInitHeatBalFiniteDiff

  RETURN

END SUBROUTINE GetCondFDInput

SUBROUTINE InitHeatBalFiniteDiff

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard J. Liesen
          !       DATE WRITTEN   Oct 2003
          !       MODIFIED       na
          !       RE-ENGINEERED  C O Pedersen 2006
          !                      B. Griffith May 2011 move begin-environment and every-timestep inits, cleanup formatting

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine sets the initial values for the FD moisture calculation

          ! METHODOLOGY EMPLOYED:

          ! REFERENCES:
          !

          ! USE STATEMENTS:
  USE DataInterfaces, ONLY:SetupOutputVariable
  USE DataSurfaces, ONLY: HeatTransferModel_CondFD

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

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  LOGICAL,SAVE        :: MyEnvrnFlag=.TRUE.
  INTEGER :: SurfNum
  INTEGER :: ConstrNum    ! Loop counter
  LOGICAL :: errorsFound

  IF (GetHBFiniteDiffInputFlag) THEN
      ! Obtains conduction FD related parameters from input file
    CALL GetCondFDInput
    GetHBFiniteDiffInputFlag=.false.
  ENDIF

  errorsFound=.false.

  ! now do begin environment inits.
  IF (BeginEnvrnFlag .AND. MyEnvrnFlag) THEN
    DO SurfNum=1,TotSurfaces
      IF (Surface(SurfNum)%HeatTransferAlgorithm /= HeatTransferModel_CondFD)   CYCLE
      IF(Surface(SurfNum)%Construction <= 0) CYCLE  ! Shading surface, not really a heat transfer surface
      ConstrNum=Surface(SurfNum)%Construction
      IF(Construct(ConstrNum)%TypeIsWindow) CYCLE  !  Windows simulated in Window module
      SurfaceFD(SurfNum)%T              = TempInitValue
      SurfaceFD(SurfNum)%TOld           = TempInitValue
      SurfaceFD(SurfNum)%TT             = TempInitValue
      SurfaceFD(SurfNum)%Rhov           = RhovInitValue
      SurfaceFD(SurfNum)%RhovOld        = RhovInitValue
      SurfaceFD(SurfNum)%RhoT           = RhovInitValue
      SurfaceFD(SurfNum)%TD             = TempInitValue
      SurfaceFD(SurfNum)%TDT            = TempInitValue
      SurfaceFD(SurfNum)%TDTLast        = TempInitValue
      SurfaceFD(SurfNum)%TDOld          = TempInitValue
      SurfaceFD(SurfNum)%TDreport       = TempInitValue
      SurfaceFD(SurfNum)%RH             = 0.0d0
      SurfaceFD(SurfNum)%RHreport       = 0.0d0
      SurfaceFD(SurfNum)%EnthOld        = EnthInitValue
      SurfaceFD(SurfNum)%EnthNew        = EnthInitValue
      SurfaceFD(SurfNum)%EnthLast       = EnthInitValue

      TempOutsideAirFD(SurfNum)         = 0.d0
      RhoVaporAirOut(SurfNum)           = 0.d0
      RhoVaporSurfIn(SurfNum)           = 0.d0
      RhoVaporAirIn(SurfNum)            = 0.d0
      HConvExtFD(SurfNum)               = 0.d0
      HMassConvExtFD(SurfNum)           = 0.d0
      HConvInFD(SurfNum)                = 0.d0
      HMassConvInFD(SurfNum)            = 0.d0
      HSkyFD(SurfNum)                   = 0.d0
      HGrndFD(SurfNum)                  = 0.d0
      HAirFD(SurfNum)                   = 0.d0
    ENDDO
    WarmupSurfTemp=0
    MyEnvrnFlag = .FALSE.
  END IF
  IF (.NOT. BeginEnvrnFlag) THEN
    MyEnvrnFlag=.TRUE.
  ENDIF

  ! now do every timestep inits

  DO SurfNum=1,TotSurfaces
    IF (Surface(SurfNum)%HeatTransferAlgorithm /= HeatTransferModel_CondFD)  CYCLE
    IF(Surface(SurfNum)%Construction <= 0) CYCLE  ! Shading surface, not really a heat transfer surface
    ConstrNum=Surface(SurfNum)%Construction
    IF(Construct(ConstrNum)%TypeIsWindow) CYCLE  !  Windows simulated in Window module
    SurfaceFD(SurfNum)%T        = SurfaceFD(SurfNum)%TOld
    SurfaceFD(SurfNum)%Rhov     = SurfaceFD(SurfNum)%RhovOld
    SurfaceFD(SurfNum)%TD       = SurfaceFD(SurfNum)%TDOld
    SurfaceFD(SurfNum)%TDT      = SurfaceFD(SurfNum)%TDreport !PT changes from TDold to TDreport
    SurfaceFD(SurfNum)%TDTLast  = SurfaceFD(SurfNum)%TDOld
    SurfaceFD(SurfNum)%EnthOld  = SurfaceFD(SurfNum)%EnthOld
    SurfaceFD(SurfNum)%EnthNew  = SurfaceFD(SurfNum)%EnthOld
    SurfaceFD(SurfNum)%EnthLast = SurfaceFD(SurfNum)%EnthOld
  ENDDO

  RETURN

END SUBROUTINE InitHeatBalFiniteDiff


SUBROUTINE InitialInitHeatBalFiniteDiff

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Linda Lawrie
          !       DATE WRITTEN   March 2012
          !       MODIFIED       na
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This routine performs the original allocate, inits and setup output variables for the
          ! module.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General,      ONLY : TrimSigDigits, RoundSigDigits
  USE DataSurfaces, ONLY : HeatTransferModel_CondFD
  USE DataHeatBalance, Only: HighDiffusivityThreshold, ThinMaterialLayerThreshold

  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:
  INTEGER :: Lay
  INTEGER :: SurfNum
  CHARACTER(len=MaxNameLength) :: LayChar

  REAL(r64)             :: dxn        ! Intermediate calculation of nodal spacing. This is the full dx. There is
                                                  ! a half dxn thick node at each surface. dxn is the "capacitor" spacing.
  INTEGER                           :: ipts1      ! Intermediate calculation for number of full thickness nodes per layer. There
                                                  ! are always two half nodes at the layer faces.
  INTEGER :: Layer      ! Loop counter
  INTEGER :: OutwardMatLayerNum ! layer index, layer outward of the current layer
  INTEGER :: layerNode
  INTEGER :: Delt
  INTEGER :: ConstrNum    ! Loop counter
  INTEGER :: TotNodes     ! Loop counter
  INTEGER :: CurrentLayer ! Loop counter
  INTEGER :: Surf         ! Loop counter
  INTEGER :: index        ! Loop Counters

  REAL(r64) :: Alpha
  REAL(r64) :: Malpha
  REAL(r64) :: stabilitytemp
  REAL(r64) :: stabilitymoist
  REAL(r64) :: a
  REAL(r64) :: b
  REAL(r64) :: c
  REAL(r64) :: d
  REAL(r64) :: kt
  REAL(r64) :: RhoS
  REAL(r64) :: Por
  REAL(r64) :: Cp
  REAL(r64) :: Dv
  LOGICAL :: errorsFound
  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

  ALLOCATE(ConstructFD(TotConstructs))
  ALLOCATE(SigmaR(TotConstructs))
  ALLOCATE(SigmaC(TotConstructs))

  ALLOCATE(SurfaceFD(TotSurfaces))
  ALLOCATE(QHeatInFlux(TotSurfaces))
  ALLOCATE(QHeatOutFlux(TotSurfaces))
!  ALLOCATE(QFluxInArrivSurfCond(TotSurfaces))
!  ALLOCATE(QFluxOutArrivSurfCond(TotSurfaces))
!  ALLOCATE(OpaqSurfInsFaceConductionFlux(TotSurfaces))
!  ALLOCATE(OpaqSurfOutsideFaceConductionFlux(TotSurfaces))
!  ALLOCATE(QFluxZoneToInSurf(TotSurfaces))
!  ALLOCATE(QFluxOutsideToOutSurf(TotSurfaces))

 ! And then initialize
  QHeatInFlux  = 0.d0
  QHeatOutFlux  = 0.d0
  !QFluxZoneToInSurf  = 0.d0
  !QFluxOutsideToOutSurf  = 0.d0
  !QFluxInArrivSurfCond  = 0.d0
  !QFluxOutArrivSurfCond  = 0.d0
  OpaqSurfInsFaceConductionFlux = 0.d0
  OpaqSurfOutsideFaceConductionFlux = 0.d0

  ! Setup Output Variables

  !  set a Delt that fits the zone time step and keeps it below 200s.

  fracTimeStepZone_Hour=1.0d0/REAL(NumOfTimeStepInHour,r64)

  DO index = 1, 20
    Delt = (fracTimeStepZone_Hour*SecInHour)/index   ! TimeStepZone = Zone time step in fractional hours
    IF ( Delt <= 200 )  EXIT
  END DO

  DO ConstrNum = 1, TotConstructs
    !Need to skip window constructions and eventually window materials
    IF(Construct(ConstrNum)%TypeIsWindow) CYCLE

    ALLOCATE(ConstructFD(ConstrNum)%Name(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%Thickness(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%NodeNumPoint(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%DelX(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%TempStability(Construct(ConstrNum)%TotLayers))
    ALLOCATE(ConstructFD(ConstrNum)%MoistStability(Construct(ConstrNum)%TotLayers))
      ! Node Number of the Interface node following each layer
!    ALLOCATE(ConstructFD(ConstrNum)%InterfaceNodeNums(Construct(ConstrNum)%TotLayers))


    TotNodes = 0
    SigmaR(ConstrNum) =0.0d0
    SigmaC(ConstrNum) =0.0d0

    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  and to.
          ! determine the total number of full size nodes in each layer.
          ! The number of temperature points is one more than this
          ! because of the two half nodes at the layer faces.
          ! 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.
          !  Change to implicit formulation still uses explicit stability, but
          ! now there are special equations for R-only layers.

      CurrentLayer = Construct(ConstrNum)%LayerPoint(Layer)


      ConstructFD(ConstrNum)%Name(Layer) = Material(CurrentLayer)%Name
      ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness

     ! Do some quick error checks for this section.

      IF (Material(CurrentLayer)%ROnly  ) THEN  ! Rlayer

        !  These values are only needed temporarily and to calculate flux,
        !   Layer will be handled
        !  as a pure R in the temperature calc.
        ! assign other properties based on resistance

        Material(CurrentLayer)%SpecHeat  = 0.0001d0
        Material(CurrentLayer)%Density = 1.d0
        Material(currentLayer)%Thickness = 0.1d0  !  arbitrary thickness for R layer
        Material(CurrentLayer)%Conductivity= &
        Material(currentLayer)%Thickness/Material(CurrentLayer)%Resistance
        kt = Material(CurrentLayer)%Conductivity
        ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness

        SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance of R layer
        SigmaC(ConstrNum) = SigmaC(ConstrNum)  + 0.0d0                               !  no capacitance for R layer

        Alpha = kt/(Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat)

        mAlpha = 0.d0

      ELSE IF (Material(CurrentLayer)%Group == 1 ) THEN   !  Group 1 = Air

                !  Again, these values are only needed temporarily and to calculate flux,
                !   Air layer will be handled
                !  as a pure R in the temperature calc.
                !
                ! assign
                ! other properties based on resistance

        Material(CurrentLayer)%SpecHeat  = 0.0001d0
        Material(CurrentLayer)%Density = 1.d0
        Material(currentLayer)%Thickness = 0.1d0  !  arbitrary thickness for R layer
        Material(currentLayer)%Conductivity= &
        Material(currentLayer)%Thickness/Material(CurrentLayer)%Resistance
        kt = Material(CurrentLayer)%Conductivity
        ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness

        SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance of R layer
        SigmaC(ConstrNum) = SigmaC(ConstrNum)  + 0.0d0                               !  no capacitance for R layer

        Alpha = kt/(Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat)
        mAlpha = 0.d0
      ELSE IF (Construct(ConstrNum)%TypeIsIRT) THEN  ! make similar to air? (that didn't seem to work well)
        CALL ShowSevereError('InitHeatBalFiniteDiff: Construction ="'//trim(Construct(ConstrNum)%Name)//  &
             '" uses Material:InfraredTransparent. Cannot be used currently with finite difference calculations.')
        IF (Construct(ConstrNum)%IsUsed) THEN
          CALL ShowContinueError('...since this construction is used in a surface, the simulation is not allowed.')
          errorsFound=.true.
        ELSE
          CALL ShowContinueError('...if this construction were used in a surface, the simulation would be terminated.')
        ENDIF
        CYCLE
!            ! set properties to get past other initializations.
!          Material(CurrentLayer)%SpecHeat  = 0.0001d0
!          Material(CurrentLayer)%Density = 0.0001d0
!          Material(currentLayer)%Thickness = 0.1d0  !  arbitrary thickness for R layer
!          Material(currentLayer)%Conductivity= &
!    -      Material(currentLayer)%Thickness/Material(CurrentLayer)%Resistance
!          kt = Material(CurrentLayer)%Conductivity
!          ConstructFD(ConstrNum)%Thickness(Layer) = Material(CurrentLayer)%Thickness
!
!          SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance of R layer
!          SigmaC(ConstrNum) = SigmaC(ConstrNum)  + 0.0                               !  no capacitance for R layer
!
!          Alpha = kt/(Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat)
!          mAlpha = 0.d0
      ELSE
       !    Regular material Properties
        a = Material(CurrentLayer)%MoistACoeff
        b = Material(CurrentLayer)%MoistBCoeff
        c = Material(CurrentLayer)%MoistCCoeff
        d = Material(CurrentLayer)%MoistDCoeff
        kt = Material(CurrentLayer)%Conductivity
        RhoS = Material(CurrentLayer)%Density
        Por = Material(CurrentLayer)%Porosity
        Cp = Material(CurrentLayer)%SpecHeat
         ! Need Resistance for reg layer
        Material(CurrentLayer)%Resistance = Material(currentLayer)%Thickness/Material(CurrentLayer)%Conductivity
        Dv = Material(CurrentLayer)%VaporDiffus
        SigmaR(ConstrNum) = SigmaR(ConstrNum) + Material(CurrentLayer)%Resistance  ! add resistance
        SigmaC(ConstrNum) = SigmaC(ConstrNum)  + &
        Material(CurrentLayer)%Density*Material(CurrentLayer)%SpecHeat*Material(CurrentLayer)%Thickness
        Alpha = kt/(RhoS*Cp)
        !   Alpha = kt*(Por+At*RhoS)/(RhoS*(Bv*Por*Lambda+Cp*(Por+At*RhoS)))
        MAlpha = 0.d0

        !check for Material layers that are too thin and highly conductivity (not appropriate for surface models)
        IF (Alpha > HighDiffusivityThreshold .AND. .NOT. Material(CurrentLayer)%WarnedForHighDiffusivity) THEN
          DeltaTimestep      = TimeStepZone * SecInHour
          ThicknessThreshold = SQRT(Alpha * DeltaTimestep * 3.d0)
          IF (Material(CurrentLayer)%Thickness < ThicknessThreshold) THEN
            CALL ShowSevereError('InitialInitHeatBalFiniteDiff: 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 by Conduction Finite Difference, ' &
                                    //' 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

      END IF  !  R, Air  or regular material properties and parameters

     ! Proceed with setting node sizes in layers

      DXN=SQRT(Alpha*Delt*SpaceDescritConstant)  ! The Fourier number is set using user constant

       ! number of nodes=thickness/spacing.  This is number of full size node spaces across layer.
      Ipts1=INT(Material(CurrentLayer)%Thickness/DXN)
       !  set high conductivity layers to a single full size node thickness. (two half nodes)
      IF(Ipts1 <= 1) Ipts1 = 1
      IF (Material(CurrentLayer)%ROnly .or. Material(CurrentLayer)%Group == 1 ) THEN

        Ipts1 = 1  !  single full node in R layers- surfaces of adjacent material or inside/outside layer
      END If

      Dxn=Material(CurrentLayer)%Thickness/REAL(Ipts1,r64)  ! full node thickness

      StabilityTemp = Alpha*Delt/dxn**2
      StabilityMoist = malpha*Delt/dxn**2
      ConstructFD(ConstrNum)%TempStability(Layer) = StabilityTemp
      ConstructFD(ConstrNum)%MoistStability(Layer) = StabilityMoist
      ConstructFD(ConstrNum)%Delx(Layer) = DXN

      TotNodes = TotNodes + Ipts1   !  number of full size nodes
      ConstructFD(ConstrNum)%NodeNumPoint(Layer) = Ipts1   !  number of full size nodes
    END DO   !  end of layer loop.


    ConstructFD(ConstrNum)%TotNodes = TotNodes
    ConstructFD(ConstrNum)%DeltaTime = Delt



  END DO            ! End of Construction Loop.  TotNodes in each construction now set

  ! now determine x location, or distance that nodes are from the outside face in meters
  DO ConstrNum = 1, TotConstructs
     IF (ConstructFD(ConstrNum)%TotNodes > 0) THEN
       ALLOCATE(ConstructFD(ConstrNum)%NodeXlocation(ConstructFD(ConstrNum)%TotNodes + 1 ))
       ConstructFD(ConstrNum)%NodeXlocation = 0.d0 ! init them all
       Ipts1 = 0 ! init counter
       DO Layer = 1, Construct(ConstrNum)%TotLayers
         OutwardMatLayerNum = Layer - 1
         DO layerNode = 1, ConstructFD(ConstrNum)%NodeNumPoint(Layer)
           Ipts1 = Ipts1 +1
           IF (Ipts1 == 1) THEN
             ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = 0.d0 ! first node is on outside face

           ELSEIF (LayerNode == 1) THEN
             IF ( OutwardMatLayerNum > 0 .AND. OutwardMatLayerNum <= Construct(ConstrNum)%TotLayers ) THEN
             ! later nodes are Delx away from previous, but use Delx from previous layer
               ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = ConstructFD(ConstrNum)%NodeXlocation(Ipts1 - 1) &
                                                           + ConstructFD(ConstrNum)%Delx(OutwardMatLayerNum)
             ENDIF
           ELSE
             ! later nodes are Delx away from previous
             ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = ConstructFD(ConstrNum)%NodeXlocation(Ipts1 - 1) &
                                                           + ConstructFD(ConstrNum)%Delx(Layer)
           ENDIF

         ENDDO
       ENDDO
       Layer = Construct(ConstrNum)%TotLayers
       Ipts1 = Ipts1 +1
       ConstructFD(ConstrNum)%NodeXlocation(Ipts1) = ConstructFD(ConstrNum)%NodeXlocation(Ipts1 - 1) &
                                                           + ConstructFD(ConstrNum)%Delx(Layer)
     ENDIF
  ENDDO

  DO Surf = 1,TotSurfaces
    IF (.not. Surface(Surf)%HeatTransSurf) CYCLE
    IF (Surface(Surf)%Class == SurfaceClass_Window) CYCLE
    IF (Surface(Surf)%HeatTransferAlgorithm /= HeatTransferModel_CondFD) CYCLE
!    IF(Surface(Surf)%Construction <= 0) CYCLE  ! Shading surface, not really a heat transfer surface
    ConstrNum=Surface(surf)%Construction
!    IF(Construct(ConstrNum)%TypeIsWindow) CYCLE  !  Windows simulated in Window module
    TotNodes = ConstructFD(ConstrNum)%TotNodes

     !Allocate the Surface Arrays
    ALLOCATE(SurfaceFD(Surf)%T(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TT(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%Rhov(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RhovOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RhoT(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TD(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDT(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDTLast(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%TDreport(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RH(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%RHreport(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%EnthOld(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%EnthNew(TotNodes + 1))
    ALLOCATE(SurfaceFD(Surf)%EnthLast(TotNodes + 1))


    !Initialize the allocated arrays.
    SurfaceFD(Surf)%T              = TempInitValue
    SurfaceFD(Surf)%TOld           = TempInitValue
    SurfaceFD(Surf)%TT             = TempInitValue
    SurfaceFD(Surf)%Rhov           = RhovInitValue
    SurfaceFD(Surf)%RhovOld        = RhovInitValue
    SurfaceFD(Surf)%RhoT           = RhovInitValue
    SurfaceFD(Surf)%TD             = TempInitValue
    SurfaceFD(Surf)%TDT            = TempInitValue
    SurfaceFD(Surf)%TDTLast        = TempInitValue
    SurfaceFD(Surf)%TDOld          = TempInitValue
    SurfaceFD(Surf)%TDreport       = TempInitValue
    SurfaceFD(Surf)%RH             = 0.0d0
    SurfaceFD(Surf)%RHreport       = 0.0d0
    SurfaceFD(Surf)%EnthOld        = EnthInitValue
    SurfaceFD(Surf)%EnthNew        = EnthInitValue
    SurfaceFD(Surf)%EnthLast       = EnthInitValue
  END DO


  DO SurfNum=1,TotSurfaces
    IF (.not. Surface(SurfNum)%HeatTransSurf) CYCLE
    IF (Surface(SurfNum)%Class == SurfaceClass_Window) CYCLE
    IF (Surface(SurfNum)%HeatTransferAlgorithm /= HeatTransferModel_CondFD) CYCLE
 !   If(SolutionAlgo == UseCondFD .or. SolutionAlgo == UseCondFDSimple)Then
!
!    CALL SetupOutputVariable('CondFD Outside Surface Heat Flux [W/m2]',   QFluxOutArrivSurfCond(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))
!    CALL SetupOutputVariable('CondFD Inside Surface Heat Flux [W/m2]',    QFluxInArrivSurfCond(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))
!    CALL SetupOutputVariable('CondFD Outside Heat Flux to Surface [W/m2]',QFluxOutsideToOutSurf(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))
!    CALL SetupOutputVariable('CondFD Inside Heat Flux to Surface [W/m2]', QFluxZoneToInSurf(SurfNum), &
!                             'Zone','State',TRIM(Surface(SurfNum)%Name))

    CALL SetupOutputVariable('CondFD Inner Solver Loop Iteration Count [ ]', SurfaceFD(SurfNum)%GSloopCounter, &
                             'Zone','Sum',Surface(SurfNum)%Name)



    TotNodes = ConstructFD(Surface(SurfNum)%Construction)%TotNodes  ! Full size nodes, start with outside face.
    DO Lay = 1,TotNodes +1  ! include inside face node
      CALL SetupOutputVariable('CondFD Surface Temperature Node '//TRIM(TrimSigDigits(Lay))//' [C]',&
           SurfaceFD(SurfNum)%TDreport(Lay),  &
          'Zone','State', Surface(SurfNum)%Name)
    END DO

  ENDDO  ! End of the Surface Loop for Report Variable Setup

  CALL ReportFiniteDiffInits  ! Report the results from the Finite Diff Inits

  RETURN

END SUBROUTINE InitialInitHeatBalFiniteDiff

SUBROUTINE CalcHeatBalFiniteDiff(Surf,TempSurfInTmp,TempSurfOutTmp)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard J. Liesen
          !       DATE WRITTEN   Oct 2003
          !       MODIFIED       Aug 2006 by C O Pedersen to include implicit solution and variable properties with
          !                                material enthalpy added for Phase Change Materials.
          !                      Sept 2010 B. Griffith, remove allocate/deallocate, use structure variables
          !                      March 2011 P. Tabares, add relaxation factor and add surfIteration to
          !                                 update TD and TDT, correct interzone partition
          !                      May 2011  B. Griffith add logging and errors when inner GS loop does not converge
          !                      November 2011 P. Tabares fixed problems with adiabatic walls/massless walls and PCM stability problems

          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! this routine controls the calculation of the fluxes and temperatures using
          !      finite difference procedures for
          !      all building surface constructs.

          ! METHODOLOGY EMPLOYED:
          ! <description>

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General,         ONLY : RoundSigDigits
  USE DataHeatBalance, ONLY : CondFDRelaxFactor
  USE DataGlobals,     ONLY : KickOffSimulation
  USE DataSurfaces,    ONLY : HeatTransferModel_CondFD

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

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER,   INTENT(In)    :: Surf
  REAL(r64), INTENT(InOut) :: TempSurfInTmp       !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64), INTENT(InOut) :: TempSurfOutTmp      !Outside Surface Temperature of each Heat Transfer Surface

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

  REAL(r64)  :: HMovInsul           !  Equiv H for TIM layer,  Comes with call to
                                                 ! EvalOutsideMovableInsulation
  INTEGER    :: RoughIndexMovInsul   ! roughness  Movable insulation
  REAL(r64)  :: AbsExt              ! exterior absorptivity  movable insulation

  INTEGER  :: I       !  Node number in construction
  INTEGER  :: J
  INTEGER  :: Lay
  INTEGER  :: ctr
  INTEGER  :: ConstrNum
  INTEGER  :: TotLayers
  INTEGER  :: TotNodes
  INTEGER  :: delt
  INTEGER  :: GSiter !  iteration counter for implicit repeat calculation
  REAL(r64) :: MaxDelTemp = 0.0d0
  INTEGER   :: NodeNum

  ConstrNum=Surface(surf)%Construction

  TotNodes = ConstructFD(ConstrNum)%TotNodes
  TotLayers = Construct(ConstrNum)%TotLayers

  TempSurfInTmp =0.0d0
  TempSurfOutTmp=0.0d0


  Delt = ConstructFD(ConstrNum)%DeltaTime    !   (seconds)

  CALL EvalOutsideMovableInsulation(Surf,HMovInsul,RoughIndexMovInsul,AbsExt)
 ! Start stepping through the slab with time.
  DO J=1,NINT((TimeStepZone*SecInHour)/Delt)  !PT testing higher time steps

    DO GSiter = 1, MaxGSiter  !  Iterate implicit equations
      SurfaceFD(Surf)%TDTLast = SurfaceFD(Surf)%TDT    !  Save last iteration's TDT (New temperature) values
      SurfaceFD(Surf)%EnthLast = SurfaceFD(Surf)%EnthNew  ! Last iterations new enthalpy value

!  Original loop version
      I= 1   !  Node counter
      DO Lay = 1, TotLayers    ! Begin layer loop ...

        !For the exterior surface node with a convective boundary condition
        IF(I == 1 .and. Lay ==1)THEN
          CALL ExteriorBCEqns(Delt,I,Lay,Surf,SurfaceFD(Surf)%T, &
                                        SurfaceFD(Surf)%TT, &
                                        SurfaceFD(Surf)%Rhov, &
                                        SurfaceFD(Surf)%RhoT, &
                                        SurfaceFD(Surf)%RH, &
                                        SurfaceFD(Surf)%TD, &
                                        SurfaceFD(Surf)%TDT,&
                                        SurfaceFD(Surf)%EnthOld, &
                                        SurfaceFD(Surf)%EnthNew, &
                                        TotNodes,HMovInsul)
        END IF

        !For the Layer Interior nodes.  Arrive here after exterior surface node or interface node

        IF(TotNodes .ne. 1) THEN

          DO Ctr=2,ConstructFD(ConstrNum)%NodeNumPoint(Lay)
            I=I+1
            CALL InteriorNodeEqns(Delt,I,Lay,Surf,SurfaceFD(Surf)%T, &
                                  SurfaceFD(Surf)%TT, &
                                  SurfaceFD(Surf)%Rhov, &
                                  SurfaceFD(Surf)%RhoT,&
                                  SurfaceFD(Surf)%RH, &
                                  SurfaceFD(Surf)%TD, &
                                  SurfaceFD(Surf)%TDT, &
                                  SurfaceFD(Surf)%EnthOld, &
                                  SurfaceFD(Surf)%EnthNew)
          END DO
        END IF

        IF(Lay < TotLayers .and. TotNodes .ne. 1) THEN
          !Interface equations for 2 capactive materials
          I=I+1
          CALL IntInterfaceNodeEqns(Delt,I,Lay,Surf,SurfaceFD(Surf)%T, &
                                   SurfaceFD(Surf)%TT, &
                                   SurfaceFD(Surf)%Rhov, &
                                   SurfaceFD(Surf)%RhoT, &
                                   SurfaceFD(Surf)%RH, &
                                   SurfaceFD(Surf)%TD, &
                                   SurfaceFD(Surf)%TDT, &
                                   SurfaceFD(Surf)%EnthOld, &
                                   SurfaceFD(Surf)%EnthNew,GSiter)

        ELSE IF (Lay == TotLayers) THEN
          !For the Interior surface node with a convective boundary condition
          I=I+1
          CALL InteriorBCEqns(Delt,I,Lay,Surf, &
                                   SurfaceFD(Surf)%T, &
                                   SurfaceFD(Surf)%TT, &
                                   SurfaceFD(Surf)%Rhov, &
                                   SurfaceFD(Surf)%RhoT, &
                                   SurfaceFD(Surf)%RH, &
                                   SurfaceFD(Surf)%TD, &
                                   SurfaceFD(Surf)%TDT, &
                                   SurfaceFD(Surf)%EnthOld, &
                                   SurfaceFD(Surf)%EnthNew, &
                                   SurfaceFD(Surf)%TDReport)
        END IF

      END DO    !The end of the layer loop


      IF (Gsiter .gt. 5) THEN
      !apply Relaxation factor for stability, use current (TDT) and previous (TDTLast) iteration temperature values
      !to obtain the actual temperature that is going to be used for next iteration. THis would mostly happen with PCM
        SurfaceFD(Surf)%TDT = SurfaceFD(Surf)%TDTLast+ (SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDTLast) * 0.5d0
      ENDIF

      IF (Gsiter .gt. 10) THEN
      !apply Relaxation factor for stability, use current (TDT) and previous (TDTLast) iteration temperature values
      !to obtain the actual temperature that is going to be used for next iteration. THis would mostly happen with PCM
        SurfaceFD(Surf)%TDT = SurfaceFD(Surf)%TDTLast+ (SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDTLast) * 0.25d0
      ENDIF

      IF (Gsiter .gt. 15) THEN
      !apply Relaxation factor for stability, use current (TDT) and previous (TDTLast) iteration temperature values
      !to obtain the actual temperature that is going to be used for next iteration. THis would mostly happen with PCM
        SurfaceFD(Surf)%TDT = SurfaceFD(Surf)%TDTLast+ (SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDTLast) * 0.10d0
      ENDIF

      ! the following could blow up when all the node temps sum to less than 1.0.  seems poorly formulated for temperature in C.
        !PT delete one zero and decrese number of minimum iterations, from 3 (which actually requires 4 iterations) to 2.

      IF (Gsiter .gt. 2  .and.ABS(SUM(SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDTLast)/SUM(SurfaceFD(Surf)%TDT)) < 0.00001d0 )  EXIT
      !SurfaceFD(Surf)%GSloopCounter = Gsiter  !PT moved out of GSloop so it can actually count all iterations

!feb2012 the following could blow up when all the node temps sum to less than 1.0.  seems poorly formulated for temperature in C.
!feb2012      IF (Gsiter .gt. 3  .and.ABS(SUM(SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDTLast)/SUM(SurfaceFD(Surf)%TDT)) < 0.000001d0 )  EXIT
!feb2012      SurfaceFD(Surf)%GSloopCounter = Gsiter
!      IF ((GSiter == MaxGSiter) .AND. (SolutionAlgo /= UseCondFDSimple)) THEN ! didn't ever converge
!        IF (.NOT. WarmupFlag .AND. (.NOT. KickOffSimulation)) THEN
!          ErrCount=ErrCount+1
!          ErrorSignal = ABS(SUM(SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDTLast)/SUM(SurfaceFD(Surf)%TDT))
!          IF (ErrCount < 10) THEN
!            CALL ShowWarningError('ConductionFiniteDifference inner iteration loop did not converge for surface named ='// &
!                          TRIM(Surface(Surf)%Name) // &
!                          ', with error signal ='//TRIM(RoundSigDigits(ErrorSignal, 8)) // &
!                          ' vs criteria of 0.000001')
!            CALL ShowContinueErrorTimeStamp(' ')
!          ELSE
!            CALL ShowRecurringWarningErrorAtEnd('ConductionFiniteDifference convergence problem continues for surface named ='// &
!                                                TRIM(Surface(Surf)%Name) , &
!                                               SurfaceFD(Surf)%GSloopErrorCount,ReportMaxOf=ErrorSignal,ReportMinOf=ErrorSignal,  &
!                                               ReportMaxUnits='[ ]',ReportMinUnits='[ ]')
!          ENDIF
!
!        ENDIF
!
!      ENDIF

    END DO ! End of Gauss Seidell iteration loop

    SurfaceFD(Surf)%GSloopCounter = Gsiter   !outputs GSloop iterations, useful for pinpointing stability issues with condFD
    IF (CondFDRelaxFactor/=1.0d0) THEN
     !apply Relaxation factor for stability, use current (TDT) and previous (TDreport) temperature values
     !   to obtain the actual temperature that is going to be exported/use
      SurfaceFD(Surf)%TDT = SurfaceFD(Surf)%TDreport+ (SurfaceFD(Surf)%TDT-SurfaceFD(Surf)%TDreport) * CondFDRelaxFactor
      SurfaceFD(Surf)%EnthOld = SurfaceFD(Surf)%EnthNew
    ENDIF

  END DO     !The end of the Time Loop   !PT solving time steps

  TempSurfOutTmp = SurfaceFD(Surf)%TDT(1)
  TempSurfInTmp  = SurfaceFD(Surf)%TDT(TotNodes+1)
  RhoVaporSurfIn(surf) = 0.0d0

  ! determine largest change in node temps
  MaxDelTemp      = 0.0d0
  DO NodeNum = 1, TotNodes+1   !need to consider all nodes
    MaxDelTemp = MAX(ABS(SurfaceFD(Surf)%TDT(NodeNum) - SurfaceFD(Surf)%TDreport(NodeNum)), MaxDelTemp)
  ENDDO
  SurfaceFD(Surf)%MaxNodeDelTemp = MaxDelTemp
!  SurfaceFD(Surf)%TDOld          = SurfaceFD(Surf)%TDT
  SurfaceFD(Surf)%TDreport       = SurfaceFD(Surf)%TDT
  SurfaceFD(Surf)%EnthOld        = SurfaceFD(Surf)%EnthNew

 RETURN


End SUBROUTINE CalcHeatBalFiniteDiff


! Beginning of Reporting subroutines
! *****************************************************************************

SUBROUTINE UpdateMoistureBalanceFD(Surf)

    ! SUBROUTINE INFORMATION:
    !   Authors:        Richard Liesen
    !   Date writtenn:  November, 2003
    !   Modified:       na
    !   Re-engineered:  na

    ! PURPOSE OF THIS SUBROUTINE:
    ! Update the data structures after the inside surface heat balance has converged.

    ! METHODOLOGY EMPLOYED:
    !

    ! USE STATEMENTS:

    IMPLICIT NONE
    Integer, INTENT(IN)  :: Surf ! Surface number

    INTEGER   :: ConstrNum

    ConstrNum = Surface(surf)%Construction
    SurfaceFD(Surf)%TOld = SurfaceFD(Surf)%T
    SurfaceFD(Surf)%RhovOld = SurfaceFD(Surf)%Rhov
    SurfaceFD(Surf)%TDOld = SurfaceFD(Surf)%TDreport

RETURN

END SUBROUTINE UpdateMoistureBalanceFD

SUBROUTINE ReportFiniteDiffInits

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard Liesen
          !       DATE WRITTEN   November 2003
          !       MODIFIED       B. Griffith, May 2011 add reporting of node x locations
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This routine gives a detailed report to the user about
          ! the initializations for the Fintie Difference calculations
          ! of each construction.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General, ONLY: ScanForReports, RoundSigDigits
  USE DataHeatBalance, ONLY: MaxAllowedDelTempCondFD, CondFDRelaxFactorInput, HeatTransferAlgosUsed

  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:
  LOGICAL :: DoReport
  CHARACTER(len=MaxNameLength) :: InodesChar
  INTEGER :: ThisNum
  INTEGER :: Layer
  INTEGER :: OutwardMatLayerNum
  INTEGER :: layerNode
  INTEGER :: Inodes

  Write(OutputFileInits,'(A)') '! <ConductionFiniteDifference HeatBalanceSettings>,Scheme Type,Space Discretization Constant,'//  &
     'Relaxation Factor,Inside Face Surface Temperature Convergence Criteria'
  Write(OutputFileInits,'(A)') ' ConductionFiniteDifference HeatBalanceSettings,'//  &
     trim(cCondFDSchemeType(CondFDSchemeType))//','//  &
     trim(RoundSigDigits(SpaceDescritConstant,2))//','//trim(RoundSigDigits(CondFDRelaxFactorInput,2))//','//  &
     trim(RoundSigDigits(MaxAllowedDelTempCondFD,4))
  CALL ScanForReports('Constructions',DoReport,'Constructions')

  IF (DoReport) THEN

!                                      Write Descriptions
    Write(OutputFileInits,'(A)') '! <Construction CondFD>,Construction Name,Index,#Layers,#Nodes,Time Step {hours}'
    Write(OutputFileInits,'(A)') '! <Material CondFD Summary>,Material Name,Thickness {m},#Layer Elements,Layer Delta X,'//  &
                             'Layer Alpha*Delt/Delx**2,Layer Moisture Stability'
  !HT Algo issue
    IF (ANY(HeatTransferAlgosUsed == UseCondFD))  &
       WRITE(OutputFileInits,'(A)') '! <ConductionFiniteDifference Node>,Node Identifier, '// &
        ' Node Distance From Outside Face {m}, Construction Name, Outward Material Name (or Face), Inward Material Name (or Face)'
    DO ThisNum=1,TotConstructs

      IF (Construct(ThisNum)%TypeIsWindow) CYCLE
      IF (Construct(ThisNum)%TypeIsIRT) CYCLE

      Write(OutputFileInits,700) TRIM(Construct(ThisNum)%Name),trim(RoundSigDigits(ThisNum)),  &
          trim(RoundSigDigits(Construct(ThisNum)%TotLayers)),  &
          trim(RoundSigDigits(Int(ConstructFD(ThisNum)%TotNodes+1))),  &
          trim(RoundSigDigits(ConstructFD(ThisNum)%DeltaTime/SecInHour,6))
!
 700  FORMAT(' Construction CondFD,',A,2(',',A),',',A,',',A)
 701  FORMAT(' Material CondFD Summary,',A,',',A,',',A,',',A,',',A,',',A)
 702  FORMAT(' ConductionFiniteDifference Node,', A, ',', A, ',', A, ',', A, ',', A)


      DO Layer=1,Construct(ThisNum)%TotLayers
            Write(OutputFileInits,701) TRIM(ConstructFD(ThisNum)%Name(Layer)),  &
                                       trim(RoundSigDigits(ConstructFD(ThisNum)%Thickness(Layer),4)),      &
                                       trim(RoundSigDigits(ConstructFD(ThisNum)%NodeNumPoint(Layer))),  &
                                       trim(RoundSigDigits(ConstructFD(ThisNum)%Delx(Layer),8)),  &
                                       trim(RoundSigDigits(ConstructFD(ThisNum)%TempStability(Layer),8)),  &
                                       trim(RoundSigDigits(ConstructFD(ThisNum)%MoistStability(Layer),8))
      ENDDO

      !now list each CondFD Node with its X distance from outside face in m along with other identifiers
      Inodes = 0

      DO Layer=1,Construct(ThisNum)%TotLayers
        OutwardMatLayerNum = Layer - 1
        DO layerNode = 1, ConstructFD(ThisNum)%NodeNumPoint(Layer)
          Inodes = Inodes + 1
          WRITE(InodesChar,*)Inodes
          IF (Inodes == 1) THEN
            WRITE(OutputFileInits,702) TRIM('Node #'//TRIM(ADJUSTL(InodesChar))), &
                                     TRIM(RoundSigDigits(ConstructFD(ThisNum)%NodeXlocation(Inodes), 8)), &
                                     TRIM(Construct(ThisNum)%Name) ,          &
                                     TRIM('Surface Outside Face'),          &
                                     TRIM(ConstructFD(ThisNum)%Name(Layer))

          ELSEIF (layerNode== 1 ) THEN

            IF ( OutwardMatLayerNum > 0 .AND. OutwardMatLayerNum <= Construct(ThisNum)%TotLayers ) THEN
               WRITE(OutputFileInits,702) TRIM('Node #'//TRIM(ADJUSTL(InodesChar))), &
                                     TRIM(RoundSigDigits(ConstructFD(ThisNum)%NodeXlocation(Inodes), 8)), &
                                     TRIM(Construct(ThisNum)%Name) ,          &
                                     TRIM(ConstructFD(ThisNum)%Name(OutwardMatLayerNum)), &
                                     TRIM(ConstructFD(ThisNum)%Name(Layer))

            ENDIF
          ELSEIF (layerNode >1) THEN
            OutwardMatLayerNum = Layer
            WRITE(OutputFileInits,702) TRIM('Node #'//TRIM(ADJUSTL(InodesChar))), &
                                     TRIM(RoundSigDigits(ConstructFD(ThisNum)%NodeXlocation(Inodes), 8)), &
                                     TRIM(Construct(ThisNum)%Name) ,          &
                                     TRIM(ConstructFD(ThisNum)%Name(OutwardMatLayerNum)), &
                                     TRIM(ConstructFD(ThisNum)%Name(Layer))
          ENDIF

          ENDDO
        ENDDO

        Layer = Construct(ThisNum)%TotLayers
        Inodes = Inodes + 1
        WRITE(InodesChar,*)Inodes
        WRITE(OutputFileInits,702) TRIM('Node #'//TRIM(ADJUSTL(InodesChar))), &
                                   TRIM(RoundSigDigits(ConstructFD(ThisNum)%NodeXlocation(Inodes), 8)), &
                                   TRIM(Construct(ThisNum)%Name) ,          &
                                   TRIM(ConstructFD(ThisNum)%Name(Layer)), &
                                   TRIM('Surface Inside Face')

    ENDDO

  ENDIF

  RETURN

END SUBROUTINE ReportFiniteDiffInits

! Utility Interpolation Function for the Module
!******************************************************************************
REAL(r64) Function terpld(N,a,x1,nind,ndep)
    !
    !author:c. o. pedersen
    !
    !purpose:
    !   this function performs a linear interpolation
    !     on a two dimensional array containing both
    !     dependent and independent variables.



    !inputs:
    !  a = two dimensional array
    !  nind=column containing independent variable
    !  ndep=column containing the dependent variable
    !   x1 = specific independent variable value for which
    !      interpolated output is wanted
    !outputs:
    !    the value of dependent variable corresponding
    !       to x1
    !    routine returns first or last dependent variable
    !      for out of range x1.
    !
     INTEGER, intent (IN)  :: N
     REAL(r64), DIMENSION(N,2),intent (IN)     :: a
     INTEGER, intent (IN)   :: nind
     INTEGER, intent (IN)   :: ndep
     REAL(r64), intent (IN)     :: x1
     INTEGER :: npts,first,last,i1,i2,i,irange
     INTEGER, DIMENSION(2) :: MaxLocArray
     REAL(r64) ::  fract

     npts = size(a,1)
     first=lbound(a,1)
     MaxLocArray=MaxLoc(a,1)
     last = MaxLocArray(1)
     IF ( npts ==1 .or. x1 <= a(first,nind) )then
       terpld=a(first,ndep)
     ELSEIF ( x1>= a(last,nind))then
       terpld=a(last,ndep)
     ELSE
       i1=first
       i2=last
       DO while ((i2-i1) > 1)
         irange=i2-i1
         i=i1+irange/2
         IF ( x1 < a(i,nind)) then
           i2=i
         ELSE
           i1=i
         END IF
      End DO
      i=i2
      fract=(x1-a(i-1,nind))/(a(i,nind)-a(i-1,nind))
      terpld=a(i-1,ndep)+fract*(a(i,ndep)-a(i-1,ndep))
    END IF
END Function

! Equation Types of the Module
!******************************************************************************
SUBROUTINE ExteriorBCEqns(Delt,I,Lay,Surf,T,TT,Rhov,RhoT,RH,TD,TDT,EnthOld,EnthNew,TotNodes,HMovInsul)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard Liesen
          !       DATE WRITTEN   November, 2003
          !       MODIFIED       B. Griffith 2010, fix adiabatic and other side surfaces
          !                      May 2011, B. Griffith, P. Tabares
          !                      November 2011 P. Tabares fixed problems with adiabatic walls/massless walls
          !                      November 2011 P. Tabares fixed problems PCM stability problems
          !       RE-ENGINEERED  Curtis Pedersen 2006

          ! PURPOSE OF THIS SUBROUTINE:
          ! <description>

          ! METHODOLOGY EMPLOYED:
          ! <description>

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE DataSurfaces,          ONLY : OtherSideCondModeledExt, OSCM, HeatTransferModel_CondFD
  USE DataHeatBalSurface  ,  ONLY : QdotRadOutRepPerArea, QdotRadOutRep, QRadOutReport
  IMPLICIT NONE ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN)  :: Delt  ! Time Increment
  INTEGER, INTENT(IN)  :: I     ! Node Index
  INTEGER, INTENT(IN)  :: Lay   ! Layer Number for Construction
  INTEGER, INTENT(IN)  :: Surf  ! Surface number
  REAL(r64),DIMENSION(:), INTENT(In)    :: T     !Old node Temperature in MFD finite difference solution
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TT    !New node Temperature in MFD finite difference solution.
  REAL(r64),DIMENSION(:), INTENT(In)    :: Rhov  !MFD Nodal Vapor Density[kg/m3] and is the old or last time step result.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RhoT  !MFD vapor density for the new time step.
  REAL(r64),DIMENSION(:), INTENT(In)    :: TD    !The old dry Temperature at each node for the CondFD algorithm..
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TDT   !The current or new Temperature at each node location for the CondFD solution..
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RH    !Nodal relative humidity
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthOld    ! Old Nodal enthalpy
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthNew    ! New Nodal enthalpy
  INTEGER, INTENT(IN)  :: TotNodes !  Total nodes in layer
  REAL(r64), INTENT(IN) :: HMovInsul  !  Conductance of movable(transparent) insulation.

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64) :: QRadSWOutFD   !Short wave radiation absorbed on outside of opaque surface
  REAL(r64) :: Delx
  INTEGER   :: ConstrNum
  INTEGER   :: MatLay
  INTEGER   :: IndVarCol
  INTEGER   :: DepVarCol
  REAL(r64) :: hconvo
  REAL(r64) :: kt  ! temperature dependent thermal conductivity,  kt=ko +kt1(T-20)
  REAL(r64) :: kto  ! Base 20C thermal conductivity
  REAL(r64) :: kt1 ! Thermal conductivity gradient coefficient where: kt=ko +kt1(T-20)
  REAL(r64) :: Cpo ! Specific heat from idf
  REAL(r64) :: Cp  ! specific heat modified if PCM, otherwise equal to Cpo
  REAL(r64) :: RhoS

  REAL(r64) :: Toa
  REAL(r64) :: Rhovo
  REAL(r64) :: hmasso

  REAL(r64) :: hgnd
  REAL(r64) :: hrad
  REAL(r64) :: hsky
  REAL(r64) :: Tgnd
  REAL(r64) :: Tsky
  REAL(r64) :: Tia
  REAL(r64) :: SigmaRLoc
  REAL(r64) :: SigmaCLoc
  REAL(r64) :: Rlayer
  REAL(r64) :: QNetSurfFromOutside  !  Combined outside surface net heat transfer terms
  REAL(r64) :: TInsulOut  ! Temperature of outisde face of Outside Insulation
  REAL(r64) :: QRadSWOutMvInsulFD  ! SW radiation at outside of Movable Insulation
  INTEGER   :: LayIn  ! layer number for call to interior eqs
  INTEGER   :: NodeIn ! node number "I" for call to interior eqs

  ConstrNum=Surface(surf)%Construction

  !Boundary Conditions from Simulation for Exterior
  hconvo = HConvExtFD(surf)
  hmasso = HMassConvExtFD(surf)

  hrad = HAirFD(surf)
  hsky = HSkyFD(surf)
  hgnd = HGrndFD(surf)

  Toa = TempOutsideAirFD(surf)
  rhovo = RhoVaporAirOut(surf)
  Tgnd = TempOutsideAirFD(surf)
  IF (Surface(Surf)%ExtBoundCond == OtherSideCondModeledExt) THEN
   !CR8046 switch modeled rad temp for sky temp.
    TSky = OSCM(Surface(Surf)%OSCMPtr)%TRad
    QRadSWOutFD = 0.0D0 ! eliminate incident shortwave on underlying surface

  ELSE
  !Set the external conditions to local variables
    QRadSWOutFD = QRadSWOutAbs(surf)
    QRadSWOutMvInsulFD = QRadSWOutMvIns(Surf)
    TSky = SkyTemp
  ENDIF
  Tia = Mat(Surface(Surf)%Zone)
  SigmaRLoc=SigmaR(ConstrNum)
  SigmaCLoc=SigmaC(ConstrNum)

  MatLay = Construct(ConstrNum)%LayerPoint(Lay)

  IF(Surface(Surf)%ExtBoundCond == Ground .or. IsRain)THEN
    TDT(I) = Toa
    TT(I) = Toa
    RhoT(I)= rhovo
  ELSEIF (Surface(Surf)%ExtBoundCond > 0) THEN
   ! this is actually the inside face of another surface, or maybe this same surface if adiabatic
   ! switch around arguments for the other surf and call routines as for interior side BC from opposite face

    LayIn = Construct(Surface(Surface(Surf)%ExtBoundCond)%Construction)%TotLayers
    nodeIn = ConstructFD(Surface(Surface(Surf)%ExtBoundCond)%Construction)%TotNodes + 1

    IF (Surface(Surf)%ExtBoundCond == surf) THEN   !adiabatic surface, PT addded since it is not the same as interzone wall
                                                   ! as Outside Boundary Condition Object can be left blank.


      CALL InteriorBCEqns(Delt,nodeIn,LayIn,Surf,SurfaceFD(Surf)%T, &
                                          SurfaceFD(Surf)%TT, &
                                          SurfaceFD(Surf)%Rhov, &
                                          SurfaceFD(Surf)%RhoT, &
                                          SurfaceFD(Surf)%RH, &
                                          SurfaceFD(Surf)%TD, &
                                          SurfaceFD(Surf)%TDT, &
                                          SurfaceFD(Surf)%EnthOld, &
                                          SurfaceFD(Surf)%EnthNew, &
                                          SurfaceFD(Surf)%TDReport)
      TDT(I)  = SurfaceFD(Surf)%TDT(TotNodes + 1)
      TT(I)   = SurfaceFD(Surf)%TT(TotNodes + 1)
      RhoT(I) = SurfaceFD(Surf)%RhoT(TotNodes + 1)

    ELSE

      CALL InteriorBCEqns(Delt,nodeIn,LayIn,Surface(Surf)%ExtBoundCond,SurfaceFD(Surface(Surf)%ExtBoundCond)%T, &
!potential-lkl-from old      CALL InteriorBCEqns(Delt,nodeIn,LayIn,Surf,SurfaceFD(Surface(Surf)%ExtBoundCond)%T, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%TT, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%Rhov, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%RhoT, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%RH, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%TD, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%TDT, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%EnthOld, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%EnthNew, &
                                            SurfaceFD(Surface(Surf)%ExtBoundCond)%TDReport)

      TDT(I)  = SurfaceFD(Surface(Surf)%ExtBoundCond)%TDT(TotNodes + 1)
      TT(I)   = SurfaceFD(Surface(Surf)%ExtBoundCond)%TT(TotNodes + 1)
      RhoT(I) = SurfaceFD(Surface(Surf)%ExtBoundCond)%RhoT(TotNodes + 1)

    ENDIF
!    CALL InteriorBCEqns(Delt,nodeIn,Layin,Surface(Surf)%ExtBoundCond,SurfaceFD(Surf)%T, &
!                                         SurfaceFD(Surf)%TT, &
!                                         SurfaceFD(Surf)%Rhov, &
!                                         SurfaceFD(Surf)%RhoT, &
!                                         SurfaceFD(Surf)%RH, &
!                                         SurfaceFD(Surf)%TD, &
!                                         SurfaceFD(Surf)%TDT, &
!                                         SurfaceFD(Surf)%EnthOld, &
!                                         SurfaceFD(Surf)%EnthNew)

    ! now fill results from interior BC model eqns into local result for current call
!    TDT(I)  = SurfaceFD(Surface(Surf)%ExtBoundCond)%TDT(TotNodes + 1)
!    TT(I)   = SurfaceFD(Surface(Surf)%ExtBoundCond)%TT(TotNodes + 1)
!    RhoT(I) = SurfaceFD(Surface(Surf)%ExtBoundCond)%RhoT(TotNodes + 1)
!    TDT(I)  = SurfaceFD(Surf)%TDT( i)
!    TT(I)   = SurfaceFD(Surf)%TT( i)
!    RhoT(I) = SurfaceFD(Surf)%RhoT( i)

    QNetSurfFromOutside = OpaqSurfInsFaceConductionFlux(Surface(Surf)%ExtBoundCond) !filled in InteriorBCEqns
!    QFluxOutsideToOutSurf(Surf)       = QnetSurfFromOutside
    OpaqSurfOutsideFaceConductionFlux(Surf)= -QnetSurfFromOutside
    OpaqSurfOutsideFaceConduction(Surf)     = Surface(Surf)%Area * OpaqSurfOutsideFaceConductionFlux(Surf)
    QHeatOutFlux(Surf)  = QnetSurfFromOutside

  ELSEIF (Surface(Surf)%ExtBoundCond <= 0) THEN   ! regular outside conditions

!++++++++++++++++++++++++++++++++++++++++++++++++++++++

    IF ( Surface(Surf)%HeatTransferAlgorithm == HeatTransferModel_CondFD  ) THEN

    ! regular outside conditions

     !  Set Thermal Conductivity.  Can be constant, simple linear temp dep or multiple linear segment temp function dep.


      kto = Material(MatLay)%Conductivity   !  20C base conductivity
      kt1 = MaterialFD(MatLay)%tk1  !  linear coefficient (normally zero)
      kt= kto + kt1*((TDT(I)+TDT(I+1))/2.d0 - 20.d0)

      IF( SUM(MaterialFD(MatLay)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function

        DepVarCol= 2  ! thermal conductivity
        IndVarCol=1  !temperature
            !  Use average temp of surface and first node for k
        kt =terpld(MaterialFD(MatLay)%numTempCond,MaterialFD(MatLay)%TempCond,(TDT(I)+TDT(I+1))/2.d0,IndVarCol,DepVarCol)

      ENDIF


      RhoS = Material(MatLay)%Density
      Cpo = Material(MatLay)%SpecHeat
      Cp=Cpo  !  Will be changed if PCM
      Delx = ConstructFD(ConstrNum)%Delx(Lay)

          !Calculate the Dry Heat Conduction Equation

      MatLay = Construct(ConstrNum)%LayerPoint(Lay)
      IF (Material(MatLay)%ROnly .or. Material(MatLay)%Group == 1 ) THEN  ! R Layer or Air Layer  **********
               !  Use algebraic equation for TDT based on R

        Rlayer = Material(MatLay)%Resistance

        TDT(I)=(QRadSWOutFD*Rlayer + TDT(I+1) + hgnd*Rlayer*Tgnd + &
                   hconvo*Rlayer*Toa + hrad*Rlayer*Toa +  &
                   hsky*Rlayer*Tsky)/  &
                   (1 + hconvo*Rlayer + hgnd*Rlayer + hrad*Rlayer + &
                    hsky*Rlayer)
        IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
            (TDT(I) < MinSurfaceTempLimit) ) THEN
          TDT(I) = MAX(MinSurfaceTempLimit,MIN(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!          CALL CheckFDSurfaceTempLimits(I,TDT(I))
        ENDIF

      ELSE  ! Regular or phase change material layer

             !  check for phase change material
        IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) >= 0.0d0) THEN    !  phase change material,  Use TempEnth Data to generate Cp
!               CheckhT = Material(MatLay)%TempEnth       ! debug

          !       Enthalpy function used to get average specific heat.  Updated by GS so enthalpy function is followed.

          DepVarCol= 2  ! enthalpy
          IndVarCol=1  !temperature
          EnthOld(I) =terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
          EnthNew(I) =terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
          IF (EnthNew(I)==EnthOld(I))  THEN
            Cp=Cpo
          ELSE
            Cp=MAX(Cpo,(EnthNew(I) -EnthOld(I))/(TDT(I)-TD(I)))
          END IF
         ELSE
           Cp = Cpo
         END IF ! Phase Change Material option

!     Choose Regular or Transparent Insulation Case

        IF(HMovInsul <= 0.0d0 ) THEN  !  regular  case

          SELECT CASE (CondFDSchemeType)
          CASE (CrankNicholsonSecondOrder)
             !  Second Order equation
            TDT(I)= (1.0d0*QRadSWOutFD + (0.5d0*Cp*Delx*RhoS*TD(I))/DelT + (0.5d0*kt*(-1.0d0*TD(I) + TD(I+1)))/Delx  &
                     + (0.5d0*kt*TDT(I+1))/Delx + 0.5d0*hgnd*Tgnd + 0.5d0*hgnd*(-1.0d0*TD(I) + Tgnd) + 0.5d0*hconvo*Toa +   &
                        0.5d0*hrad*Toa  &
                     + 0.5d0*hconvo*(-1.d0*TD(I) + Toa) + 0.5d0*hrad*(-1.d0*TD(I) + Toa) + 0.5d0*hsky*Tsky +   &
                        0.5d0*hsky*(-1.d0*TD(I) + Tsky))/  &
                       (0.5d0*hconvo + 0.5d0*hgnd + 0.5d0*hrad + 0.5d0*hsky + (0.5d0*kt)/Delx + (0.5d0*Cp*Delx*RhoS)/DelT)
!feb2012            TDT(I)= (1.0d0*QRadSWOutFD + (0.5d0*Cp*Delx*RhoS*TD(I))/DelT + (0.5d0*kt*(-1.0d0*TD(I) + TD(I+1)))/Delx  &
!feb2012                     + (0.5d0*kt*TDT(I+1))/Delx + 0.5d0*hgnd*Tgnd + 0.5d0*hgnd*(-1.0d0*TD(I) + Tgnd) + 0.5d0*hconvo*Toa +   &
!feb2012                        0.5d0*hrad*Toa  &
!feb2012                     + 0.5d0*hconvo*(-1.d0*TD(I) + Toa) + 0.5d0*hrad*(-1.d0*TD(I) + Toa) + 0.5d0*hsky*Tsky +   &
!feb2012                        0.5d0*hsky*(-1.d0*TD(I) + Tsky))/  &
!feb2012                       (0.5d0*hconvo + 0.5d0*hgnd + 0.5d0*hrad + 0.5d0*hsky + (0.5d0*kt)/Delx + (0.5d0*Cp*Delx*RhoS)/DelT)

          CASE (FullyImplicitFirstOrder)
           !   First Order
            TDT(I)=(2.0d0*Delt*Delx*QRadSWOutFD + Cp*Delx**2*RhoS*TD(I) &
                    + 2.0d0*Delt*kt*TDT(I+1) + 2.0d0*Delt*Delx*hgnd*Tgnd  &
                    + 2.0d0*Delt*Delx*hconvo*Toa + 2.0d0*Delt*Delx*hrad*Toa + 2.0d0*Delt*Delx*hsky*Tsky)/  &
                     (2.0d0*Delt*Delx*hconvo + 2.0d0*Delt*Delx*hgnd + 2.0d0*Delt*Delx*hrad +  &
                       2.0d0*Delt*Delx*hsky + 2.0d0*Delt*kt + Cp*Delx**2*RhoS)

          END SELECT


        ELSE IF ( HMovInsul > 0.0d0) Then      !  Transparent insulation on outside
      !  Transparent insulaton additions

        !Movable Insulation Layer Outside surface temp

          TInsulOut = (QRadSWOutMvInsulFD + hgnd*Tgnd + HmovInsul*TDT(I) + &
                   hconvo*Toa + hrad*Toa + hsky*Tsky)/  &
                 (hconvo + hgnd + HmovInsul + hrad + hsky)

          !List(List(Rule(TDT,(2*Delt*Delx*QradSWOutAbs +
          !-       Cp*Delx**2*Rhos*TD + 2*Delt*kt*TDTP1 +
          !-       2*Delt*Delx*HmovInsul*Tiso)/
          !-     (2*Delt*Delx*HmovInsul + 2*Delt*kt + Cp*Delx**2*Rhos))))


         ! Wall first node temperature behind Movable insulation
          SELECT CASE (CondFDSchemeType)
          CASE (CrankNicholsonSecondOrder)
            TDT(I)= (2*Delt*Delx*QradSWOutFD +  &
                     Cp*Delx**2*Rhos*TD(I) + 2*Delt*kt*TDT(I+1) +  &
                      2*Delt*Delx*HmovInsul*TInsulOut)/  &
                       (2*Delt*Delx*HmovInsul + 2*Delt*kt + Cp*Delx**2*Rhos)

          CASE (FullyImplicitFirstOrder)
            ! Currently same as Crank Nicholson, need fully implicit formulation
            TDT(I)= (2*Delt*Delx*QradSWOutFD +  &
                     Cp*Delx**2*Rhos*TD(I) + 2*Delt*kt*TDT(I+1) +  &
                      2*Delt*Delx*HmovInsul*TInsulOut)/  &
                       (2*Delt*Delx*HmovInsul + 2*Delt*kt + Cp*Delx**2*Rhos)

          END SELECT

        End IF   !  Regular layer or Movable insulation cases

        IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
            (TDT(I) < MinSurfaceTempLimit) ) THEN
          TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!          CALL CheckFDSurfaceTempLimits(I,TDT(I))
        ENDIF


      END IF  ! R layer or Regular layer


    END IF   !End IF--ELSE SECTION (regular detailed FD part or SigmaR SigmaC part

    !  Determine net heat flux to ooutside face.
    !One formulation that works for Fully Implicit and CrankNicholson and massless wall

    QNetSurfFromOutside = QRadSWOutFD + (hgnd*(-TDT(I) + Tgnd) +  &
             hconvo*(-TDT(I) + Toa) + hrad*(-TDT(I) + Toa) + hsky*(-TDT(I) + Tsky))

    !Same sign convention as CTFs
    OpaqSurfOutsideFaceConductionFlux(Surf) = -QnetSurfFromOutside
    OpaqSurfOutsideFaceConduction(Surf)     = Surface(Surf)%Area * OpaqSurfOutsideFaceConductionFlux(Surf)

    !Report all outside BC heat fluxes
    QdotRadOutRepPerArea(Surf) = -(hgnd*(TDT(I) - Tgnd) +  hrad*(TDT(I) - Toa) + hsky*( TDT(I) - Tsky))
    QdotRadOutRep(Surf)        = Surface(Surf)%Area * QdotRadOutRepPerArea(Surf)
    QRadOutReport(Surf)     = QdotRadOutRep(Surf) * SecInHour * TimeStepZone

  END IF  !End IF --ELSE SECTION (regular BC part of the ground and Rain check)


RETURN

END SUBROUTINE ExteriorBCEqns

SUBROUTINE InteriorNodeEqns(Delt,I,Lay,Surf,T,TT,Rhov,RhoT,RH,TD,TDT,EnthOld,EnthNew)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard Liesen
          !       DATE WRITTEN   November, 2003
          !       MODIFIED       May 2011, B. Griffith and P. Tabares
          !       RE-ENGINEERED  C. O. Pedersen, 2006

          ! PURPOSE OF THIS SUBROUTINE:
          ! <description>

          ! METHODOLOGY EMPLOYED:
          ! <description>

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN)  :: Delt  ! Time Increment
  INTEGER, INTENT(IN)  :: I     ! Node Index
  INTEGER, INTENT(IN)  :: Lay   ! Layer Number for Construction
  INTEGER, INTENT(IN)  :: Surf  ! Surface number
  REAL(r64),DIMENSION(:), INTENT(In)    :: T     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TT    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(In)    :: Rhov     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RhoT  !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(In)    :: TD     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TDT    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RH    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthOld    ! Old Nodal enthalpy
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthNew    ! New Nodal enthalpy


          ! SUBROUTINE PARAMETER DEFINITIONS:
!  REAL(r64), PARAMETER :: NinetyNine=99.0d0

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          ! na

  REAL(r64) :: Delx

  INTEGER :: ConstrNum
  INTEGER :: MatLay
  INTEGER :: DepVarCol
  INTEGER :: IndVarCol

  REAL(r64) :: kt  !  Thermal conductivity in temperature equation
  REAL(r64) :: ktA1  !  Variable Outer Thermal conductivity in temperature equation
  REAL(r64) :: ktA2  !  Thermal Inner conductivity in temperature equation
  REAL(r64) :: kto !  base 20 C thermal conductivity
  REAL(r64) :: kt1 !  temperature coefficient for simple temp dep k.
  REAL(r64) :: Cp   !  Cp used
  REAL(r64) :: Cpo  !  Const Cp from input
  REAL(r64) :: RhoS

  ConstrNum=Surface(surf)%Construction

  MatLay = Construct(ConstrNum)%LayerPoint(Lay)

   !  Set Thermal Conductivity.  Can be constant, simple linear temp dep or multiple linear segment temp function dep.
  kto = Material(MatLay)%Conductivity   !  20C base conductivity
  kt1 = MaterialFD(MatLay)%tk1  !  linear coefficient (normally zero)
  kt  = kto + kt1*((TDT(I)+TDT(I-1))/2.d0 - 20.0d0)
  ktA1 = kto + kt1*((TDT(I)+TDT(I+1))/2.d0 - 20.0d0)   ! Will be overridden if variable k
  ktA2 = kto + kt1*((TDT(I)+TDT(I-1))/2.d0 - 20.0d0)   ! Will be overridden if variable k

  IF( SUM(MaterialFD(MatLay)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function

    DepVarCol = 2  ! thermal conductivity
    IndVarCol = 1  !temperature
    ktA1 = terpld(MaterialFD(MatLay)%numTempCond,MaterialFD(MatLay)%TempCond,(TDT(I)+TDT(I+1))/2.0d0 ,IndVarCol,DepVarCol)
    ktA2 = terpld(MaterialFD(MatLay)%numTempCond,MaterialFD(MatLay)%TempCond,(TDT(I)+TDT(I-1))/2.0d0 ,IndVarCol,DepVarCol)
  ENDIF

  RhoS = Material(MatLay)%Density
  Cpo  = Material(MatLay)%SpecHeat
  Cp   = Cpo  ! Will be changed if PCM
  Delx = ConstructFD(ConstrNum)%Delx(Lay)

  IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) >= 0.0d0) THEN    !  phase change material,  Use TempEnth Data

    DepVarCol = 2  ! enthalpy
    IndVarCol = 1  !temperature
    EnthOld(I) =terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
    EnthNew(I) =terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
    IF (EnthNew(I)==EnthOld(I))  THEN
      Cp = Cpo
    ELSE
      Cp = MAX(Cpo,(EnthNew(I) -EnthOld(I))/(TDT(I)-TD(I)))
    END IF

  ELSE  ! No phase change

    Cp = Cpo

  END IF ! Phase Change case

  SELECT CASE (CondFDSchemeType)

  CASE (CrankNicholsonSecondOrder)
    ! Adams-Moulton second order
    TDT(I)= ((Cp*Delx*RhoS*TD(I))/Delt +   &
              0.5d0*((ktA2*(-1.0d0*TD(I) + TD(I-1)))/Delx + (ktA1*(-1.d0*TD(I) + TD(I+1)))/Delx) +   &
              (0.5d0*ktA2*TDT(I-1))/Delx + (0.5d0*ktA1*TDT(I+1))/Delx)/  &
              ((0.5d0*(ktA1+ktA2))/Delx + (Cp*Delx*RhoS)/Delt)


  CASE (FullyImplicitFirstOrder)
    ! Adams-Moulton First order
        TDT(I)= ((Cp*Delx*RhoS*TD(I))/Delt +   &
                 (1.d0*ktA2*TDT(I-1))/Delx + (1.d0*ktA1*TDT(I+1))/Delx)/  &
                 ((2.d0*(ktA1+ktA2)/2.d0)/Delx + (Cp*Delx*RhoS)/Delt)


  END SELECT

 IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
     (TDT(I) < MinSurfaceTempLimit) ) THEN
   TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!   CALL CheckFDSurfaceTempLimits(I,TDT(I))
  ENDIF

RETURN

END SUBROUTINE InteriorNodeEqns


SUBROUTINE IntInterfaceNodeEqns(Delt,I,Lay,Surf,T,TT,Rhov,RhoT,RH,TD,TDT,EnthOld,EnthNew,GSiter)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard Liesen
          !       DATE WRITTEN   November, 2003
          !       MODIFIED       May 2011, B. Griffith, P. Tabares,  add first order fully implicit, bug fixes, cleanup
          !       RE-ENGINEERED  Curtis Pedersen, Changed to Implit mode and included enthalpy.  FY2006

          ! PURPOSE OF THIS SUBROUTINE:
          ! calculate finite difference heat transfer for nodes that interface two different material layers inside construction

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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

          ! SUBROUTINE ARGUMENT DEFINITIONS:

  INTEGER, INTENT(IN)  :: Delt  ! Time Increment
  INTEGER, INTENT(IN)  :: I     ! Node Index
  INTEGER, INTENT(IN)  :: Lay   ! Layer Number for Construction
  INTEGER, INTENT(IN)  :: Surf  ! Surface number
  INTEGER, INTENT(IN)  :: GSiter  ! Iteration number of Gauss Seidell iteration
  REAL(r64),DIMENSION(:), INTENT(In)    :: T     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TT    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(In)    :: Rhov     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RhoT  !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(In)    :: TD     !OLD NODE TEMPERATURES OF EACH HEAT TRANSFER SURF IN CONDFD.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TDT    !NEW NODE TEMPERATURES OF EACH HEAT TRANSFER SURF IN CONDFD.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RH    !RELATIVE HUMIDITY.
  REAL(r64),DIMENSION(:), INTENT(In)    :: EnthOld    ! Old Nodal enthalpy
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthNew    ! New Nodal enthalpy

          ! SUBROUTINE PARAMETER DEFINITIONS:
!  REAL(r64), PARAMETER :: NinetyNine=99.0d0

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

  INTEGER :: ConstrNum
  INTEGER :: MatLay
  INTEGER :: MatLay2
  INTEGER :: IndVarCol
  INTEGER :: DepVarCol

  REAL(r64) :: kt1
  REAL(r64) :: kt2
  REAL(r64) :: kt1o
  REAL(r64) :: kt2o
  REAL(r64) :: kt11
  REAL(r64) :: kt21
  REAL(r64) :: Cp1
  REAL(r64) :: Cp2
  REAL(r64) :: Cpo1
  REAL(r64) :: Cpo2
  REAL(r64) :: RhoS1
  REAL(r64) :: RhoS2

  REAL(r64) :: Delx1
  REAL(r64) :: Delx2
  REAL(r64) :: Enth1New
  REAL(r64) :: Enth2New
  REAL(r64) :: Enth1Old
  REAL(r64) :: Enth2Old

  REAL(r64) :: Rlayer   !  resistance value of R Layer
  REAL(r64) :: Rlayer2  !  resistance value of next layer to inside
  REAL(r64) :: QSSFlux  !  Source/Sink flux value at a layer interface
  LOGICAL   :: RlayerPresent  = .FALSE.
  LOGICAL   :: RLayer2Present = .FALSE.

  ConstrNum=Surface(surf)%Construction

  MatLay = Construct(ConstrNum)%LayerPoint(Lay)
  MatLay2 = Construct(ConstrNum)%LayerPoint(Lay+1)
  !  Set Thermal Conductivity.  Can be constant, simple linear temp dep or multiple linear segment temp function dep.

  kt1o = Material(MatLay)%Conductivity
  kt11=  MaterialFD(MatLay)%tk1
  kt1= kt1o + kt11*((TDT(I)+TDT(I-1))/2.d0 -20.d0)

  IF( SUM(MaterialFD(MatLay)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function

    IndVarCol = 1  ! temperature
    DepVarCol = 2  ! thermal conductivity
    kt1       = terpld(MaterialFD(MatLay)%numTempCond,MaterialFD(MatLay)%TempCond,(TDT(I)+TDT(I-1))/2.d0,IndVarCol,DepVarCol)

  ENDIF

  RhoS1 = Material(MatLay)%Density
  Cpo1 = Material(MatLay)%SpecHeat   ! constant Cp from input file
  Delx1 = ConstructFD(ConstrNum)%Delx(Lay)
  Rlayer = Material(MatLay)%Resistance

  kt2o = Material(MatLay2)%Conductivity
  kt21=  MaterialFD(MatLay2)%tk1
  kt2= kt2o + kt21*((TDT(I)+TDT(I+1))/2.d0 -20.d0)

  IF( SUM(MaterialFD(MatLay2)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function

    IndVarCol = 1  ! temperature
    DepVarCol = 2  ! thermal conductivity
    kt2 =terpld(MaterialFD(MatLay2)%numTempCond,MaterialFD(MatLay2)%TempCond,(TDT(I)+TDT(I+1))/2.d0,IndVarCol,DepVarCol)

  ENDIF

  RhoS2 = Material(MatLay2)%Density
  Cpo2 = Material(MatLay2)%SpecHeat
  Delx2 = ConstructFD(ConstrNum)%Delx(Lay+1)
  Rlayer2 = Material(MatLay2)%Resistance
  Cp1=Cpo1   !  Will be reset if PCM
  Cp2=Cpo2   !  will be reset if PCM

  IF ( Surface(Surf)%HeatTransferAlgorithm ==  HeatTransferModel_CondFD )THEN  !HT Algo issue
              !Calculate the Dry Heat Conduction Equation
    RLayerPresent = .FALSE.
    RLayer2Present = .false.

!     Source/Sink Flux Capability ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    QSSFlux = 0.0d0
    IF (Surface(Surf)%Area >0.0d0 .and.   &
        Construct(ConstrNum)%SourceSinkPresent .and. Lay == Construct(ConstrNum)%SourceAfterLayer) then

      QSSFlux = QRadSysSource(Surf)/Surface(Surf)%Area   &
                 + QPVSysSource(Surf)/Surface(Surf)%Area     ! Includes QPV Source

    END IF

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    IF (Material(MatLay)%ROnly .or. Material(MatLay)%Group == 1)  RLayerPresent = .TRUE.

    IF (Material(MatLay2)%ROnly .or. Material(MatLay2)%Group ==1) Rlayer2Present =.TRUE.

    IF ( RLayerPresent .and. RLayer2Present) THEN
      TDT(I)=(Rlayer2*TDT(I-1) + Rlayer*TDT(I+1))/(Rlayer + Rlayer2)   ! two adjacent R layers

    ELSE IF ( RLayerPresent .and. .not. RLayer2Present  ) THEN  ! R-layer first

     !  Check for PCM second layer
     IndVarCol = 1  ! temperature
     DepVarCol = 2  ! thermal conductivity

      IF( Sum(MaterialFD(MatLay)%TempEnth(1:3,2)) < 0.0d0  &
             .and. Sum(MaterialFD(MatLay2)%TempEnth(1:3,2)) > 0.0d0) THEN    !  phase change material Layer2,  Use TempEnth Data

        Enth2Old = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TD(I),IndVarCol,DepVarCol)
        Enth2New = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TDT(I),IndVarCol,DepVarCol)
        EnthNew(I) = Enth2New  !  This node really doesn't have an enthalpy, this gives it a value

        IF (ABS(Enth2New-Enth2Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp2 = Cpo2
        ELSE
          Cp2 = MAX(Cpo2,(Enth2New -Enth2Old)/(TDT(I)-TD(I)))
        END IF
      ENDIF

         ! R layer first, then PCM or regular layer.
      SELECT CASE (CondFDSchemeType)
        CASE (CrankNicholsonSecondOrder)

          TDT(I) =( 2.d0*delt*Delx2*QSSFlux*Rlayer - delt*Delx2*TD(I) - delt*kt2*Rlayer*TD(I) +  &
                Cp2*Delx2**2*RhoS2*Rlayer*TD(I) + delt*Delx2*TD(I-1) + delt*kt2*Rlayer*TD(I+1) + delt*Delx2*TDT(I-1) +  &
                delt*kt2*Rlayer*TDT(I+1))/(delt*Delx2 + delt*kt2*Rlayer + Cp2*Delx2**2.d0*RhoS2*Rlayer)

        CASE (FullyImplicitFirstOrder)

          TDT(I) =( 2.d0*delt*Delx2*QSSFlux*Rlayer + Cp2*Delx2**2*RhoS2*Rlayer*TD(I) + 2.d0*delt*Delx2*TDT(I-1) +  &
                 2.d0*delt*kt2*Rlayer*TDT(I+1))/(2.d0*delt*Delx2 + 2.d0*delt*kt2*Rlayer + Cp2*Delx2**2.d0*RhoS2*Rlayer)

      END SELECT

      IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
          (TDT(I) < MinSurfaceTempLimit) ) THEN
        TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!        CALL CheckFDSurfaceTempLimits(I,TDT(I))
      ENDIF

    ELSE IF (.not. RLayerPresent .and.  RLayer2Present) THEN  ! R-layer second

      !  check for PCM layer before R layer
     IndVarCol = 1  ! temperature
     DepVarCol = 2  ! thermal conductivity

      IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) > 0.0d0  &
          .and. SUM(MaterialFD(MatLay2)%TempEnth(1:3,2)) < 0.0d0) THEN    !  phase change material Layer1,  Use TempEnth Data

        Enth1Old = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
        Enth1New = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
        EnthNew(I) = Enth1New    !  This node really doesn't have an enthalpy, this gives it a value

        IF (ABS(Enth1New-Enth1Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp1=Cpo1
        ELSE
          Cp1=Max(Cpo1,(Enth1New -Enth1Old)/(TDT(I)-TD(I)))
        END IF

      ENDIF

      SELECT CASE (CondFDSchemeType)

        CASE (CrankNicholsonSecondOrder)
          TDT(I)=(2.d0*delt*Delx1*QSSFlux*Rlayer2 - delt*Delx1*TD(I) -  &
                delt*kt1*Rlayer2*TD(I) + Cp1*Delx1**2*RhoS1*Rlayer2*TD(I) + delt*kt1*Rlayer2*TD(I-1) +  &
                delt*Delx1*TD(I+1) + delt*kt1*Rlayer2*TDT(I-1) + delt*Delx1*TDT(I+1))/  &
               (delt*Delx1 + delt*kt1*Rlayer2 + Cp1*Delx1**2*RhoS1*Rlayer2)
        CASE (FullyImplicitFirstOrder)
          TDT(I)=(2.d0*delt*Delx1*QSSFlux*Rlayer2  + Cp1*Delx1**2*RhoS1*Rlayer2*TD(I) + 2.d0*delt*kt1*Rlayer2*TDT(I-1) + &
            2.d0*delt*Delx1*TDT(I+1))/ (2.d0*delt*Delx1 + 2.d0*delt*kt1*Rlayer2 + Cp1*Delx1**2*RhoS1*Rlayer2)

      END SELECT

      IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
          (TDT(I) < MinSurfaceTempLimit) ) THEN
        TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!        CALL CheckFDSurfaceTempLimits(I,TDT(I))
      ENDIF

    ELSE  !   Regular or Phase Change on both sides of interface
          !   Consider the various PCM material location cases
      Cp1 = Cpo1    !  Will be changed if PCM
      Cp2 = Cpo2    !  Will be changed if PCM
      IndVarCol = 1  ! temperature
      DepVarCol = 2  ! thermal conductivity

      IF( Sum(MaterialFD(MatLay)%TempEnth(1:3,2)) > 0.0d0  &
          .and. Sum(MaterialFD(MatLay2)%TempEnth(1:3,2)) > 0.0d0) THEN    !  phase change material both layers,  Use TempEnth Data

        Enth1Old = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
        Enth2Old = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TD(I),IndVarCol,DepVarCol)
        Enth1New = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
        Enth2New = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TDT(I),IndVarCol,DepVarCol)

        EnthNew(I) = Enth1New    !  This node really doesn't have an enthalpy, this gives it a value

        IF (ABS(Enth1New-Enth1Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp1 = Cpo1
        ELSE
          Cp1 = MAX(Cpo1,(Enth1New -Enth1Old)/(TDT(I)-TD(I)))
        END IF

        IF (ABS(Enth2New-Enth2Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp2 = Cpo2
        ELSE
          Cp2 = Max(Cpo2,(Enth2New -Enth2Old)/(TDT(I)-TD(I)))
        END IF

      ELSE IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) > 0.0d0  &
                .and. SUM(MaterialFD(MatLay2)%TempEnth(1:3,2)) < 0.0d0) THEN    !  phase change material Layer1,  Use TempEnth Data

        Enth1Old = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)
        Enth1New = terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
        EnthNew(I) = Enth1New    !  This node really doesn't have an enthalpy, this gives it a value

        IF (ABS(Enth1New-Enth1Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp1 = Cpo1
        ELSE
          Cp1 = Max(Cpo1,(Enth1New -Enth1Old)/(TDT(I)-TD(I)))
        END IF

        Cp2 = Cpo2

      ELSE IF( SUM(MaterialFD(MatLay)%TempEnth(1:3,2)) < 0.0d0  &
              .and. SUM(MaterialFD(MatLay2)%TempEnth(1:3,2)) > 0.0d0) THEN    !  phase change material Layer2,  Use TempEnth Data

        Enth2Old = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TD(I),IndVarCol,DepVarCol)
        Enth2New = terpld(MaterialFD(MatLay2)%numTempEnth,MaterialFD(MatLay2)%TempEnth,TDT(I),IndVarCol,DepVarCol)
        EnthNew(I) = Enth2New  !  This node really doesn't have an enthalpy, this gives it a value

        IF (ABS(Enth2New-Enth2Old) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp2 = Cpo2
        ELSE
          Cp2 = MAX(Cpo2,(Enth2New -Enth2Old)/(TDT(I)-TD(I)))
        END IF

        Cp1 = Cpo1

      END If  ! Phase change material check

      SELECT CASE (CondFDSchemeType)

        CASE (CrankNicholsonSecondOrder)
        !     Regular Internal Interface Node with Source/sink using Adams Moulton second order
          TDT(I) = (2.d0*delt*Delx1*Delx2*QSSFlux - delt*Delx2*kt1*TD(I) - delt*Delx1*kt2*TD(I) +  &
                 Cp1*Delx1**2.d0*Delx2*RhoS1*TD(I) + Cp2*Delx1*Delx2**2.d0*RhoS2*TD(I) + delt*Delx2*kt1*TD(I-1) +  &
                 delt*Delx1*kt2*TD(I+1) + delt*Delx2*kt1*TDT(I-1) + delt*Delx1*kt2*TDT(I+1))/   &
                 (delt*Delx2*kt1 + delt*Delx1*kt2 + Cp1*Delx1**2.d0*Delx2*RhoS1 + Cp2*Delx1*Delx2**2.d0*RhoS2)

        CASE (FullyImplicitFirstOrder)
        ! first order adams moulton
          TDT(I) = (2.d0*delt*Delx1*Delx2*QSSFlux + &
                 Cp1*Delx1**2.d0*Delx2*RhoS1*TD(I) + Cp2*Delx1*Delx2**2.d0*RhoS2*TD(I) +  &
                 2.d0*delt*Delx2*kt1*TDT(I-1) + 2.d0*delt*Delx1*kt2*TDT(I+1))/   &
                 (2.d0*delt*Delx2*kt1 + 2.d0*delt*Delx1*kt2 + Cp1*Delx1**2.d0*Delx2*RhoS1 + Cp2*Delx1*Delx2**2.d0*RhoS2)

      END SELECT

      IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
          (TDT(I) < MinSurfaceTempLimit) ) THEN
        TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!        CALL CheckFDSurfaceTempLimits(I,TDT(I))
      ENDIF


      IF (Construct(ConstrNum)%SourceSinkPresent .and. Lay == Construct(ConstrNum)%SourceAfterLayer) Then
        TcondFDSourceNode(Surf)= TDT(I) ! transfer node temp to Radiant System
        TempSource(Surf) = TDT(I)  !  Transfer node temp to DataHeatBalSurface  module.

      ENDIF
!+++++++++++++++++++++++++++++++
    END IF   !  end of R-layer and Regular check

  END IF  !End of the CondFD if block


RETURN

END SUBROUTINE IntInterfaceNodeEqns


SUBROUTINE InteriorBCEqns(Delt,I,Lay,Surf,T,TT,Rhov,RhoT,RH,TD,TDT,EnthOld,EnthNew,TDReport)
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Richard Liesen
          !       DATE WRITTEN   November, 2003
          !       MODIFIED       B. Griffith, P. Tabares, May 2011, add first order fully implicit, bug fixes, cleanup
          !                      November 2011 P. Tabares fixed problems with adiabatic walls/massless walls
          !                      November 2011 P. Tabares fixed problems PCM stability problems
          !       RE-ENGINEERED  C. O. Pedersen 2006

          ! PURPOSE OF THIS SUBROUTINE:
          ! Calculate the heat transfer at the node on the surfaces inside face (facing zone)

          ! METHODOLOGY EMPLOYED:
          ! <description>

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE DataHeatBalFanSys, ONLY: Mat, ZoneAirHumRat, QHTRadSysSurf, QHWBaseboardSurf, QSteamBaseboardSurf, QElecBaseboardSurf
  USE DataSurfaces,      ONLY: HeatTransferModel_CondFD
  IMPLICIT NONE ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN)  :: Delt  ! Time Increment
  INTEGER, INTENT(IN)  :: I     ! Node Index
  INTEGER, INTENT(IN)  :: Lay   ! Layer Number for Construction
  INTEGER, INTENT(IN)  :: Surf  ! Surface number
  REAL(r64),DIMENSION(:), INTENT(In)    :: T     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF (Old).
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TT    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF (New).
  REAL(r64),DIMENSION(:), INTENT(In)    :: Rhov     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RhoT  !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(In)    :: TD     !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TDT    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: RH    !INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthOld    ! Old Nodal enthalpy
  REAL(r64),DIMENSION(:), INTENT(InOut) :: EnthNew    ! New Nodal enthalpy
  REAL(r64),DIMENSION(:), INTENT(InOut) :: TDreport    ! Temperature value from previous HeatSurfaceHeatManager titeration's value

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64) :: NetLWRadToSurfFD !Net interior long wavelength radiation to surface from other surfaces
  REAL(r64) :: QRadSWInFD    !Short wave radiation absorbed on inside of opaque surface
  REAL(r64) :: QHtRadSysSurfFD  ! Current radiant heat flux at a surface due to the presence of high temperature radiant heaters
  Real(r64) :: QHWBaseboardSurfFD  ! Current radiant heat flux at a surface due to the presence of hot water baseboard heaters
  Real(r64) :: QSteamBaseboardSurfFD  ! Current radiant heat flux at a surface due to the presence of steam baseboard heaters
  Real(r64) :: QElecBaseboardSurfFD  ! Current radiant heat flux at a surface due to the presence of electric baseboard heaters
  REAL(r64) :: QRadThermInFD !Thermal radiation absorbed on inside surfaces
  REAL(r64) :: Delx
  REAL(r64), PARAMETER :: IterDampConst = 5.0d0  ! Damping constant for inside surface temperature iterations. Only used for massless (R-value only) Walls

  INTEGER    :: ConstrNum
  INTEGER    :: MatLay
  INTEGER    :: IndVarCol
  INTEGER    :: DepVarCol

  REAL(r64)  :: kto
  REAL(r64)  :: kt1
  REAL(r64)  :: kt
  REAL(r64)  :: Cp
  REAL(r64)  :: Cpo
  REAL(r64)  :: RhoS
  REAL(r64)  :: Tia
  REAL(r64)  :: Rhovi
  REAL(r64)  :: hmassi
  REAL(r64)  :: hconvi

  REAL(r64)  :: Rlayer
  REAL(r64)  :: SigmaRLoc
  REAL(r64)  :: SigmaCLoc
  REAL(r64)  :: QNetSurfInside


  ConstrNum = Surface(surf)%Construction
  SigmaRLoc = SigmaR(ConstrNum)
  SigmaCLoc = SigmaC(ConstrNum)

  !Set the internal conditions to local variables
  NetLWRadToSurfFD      = NetLWRadToSurf(surf)
  QRadSWInFD            = QRadSWInAbs(surf)
  QHtRadSysSurfFD       = QHtRadSysSurf(surf)
  QHWBaseboardSurfFD    = QHWBaseboardSurf(surf)
  QSteamBaseboardSurfFD = QSteamBaseboardSurf(surf)
  QElecBaseboardSurfFD  = QElecBaseboardSurf(surf)
  QRadThermInFD         = QRadThermInAbs(Surf)

  !Boundary Conditions from Simulation for Interior
  hconvi = HConvInFD(surf)
  hmassi = HMassConvInFD(surf)

  Tia    = Mat(Surface(Surf)%Zone)
  rhovi  = RhoVaporAirIn(surf)

!++++++++++++++++++++++++++++++++++++++++++++++++++++++
   !    Do all the nodes in the surface   Else will switch to SigmaR,SigmaC
  IF (Surface(Surf)%HeatTransferAlgorithm == HeatTransferModel_CondFD) THEN

    MatLay = Construct(ConstrNum)%LayerPoint(Lay)
     !  Set Thermal Conductivity.  Can be constant, simple linear temp dep or multiple linear segment temp function dep.
    kto = Material(MatLay)%Conductivity   !  20C base conductivity
    kt1 = MaterialFD(MatLay)%tk1  !  linear coefficient (normally zero)
    kt  = kto + kt1*((TDT(I)+TDT(i-1))/2.d0 - 20.d0)


    IF( SUM(MaterialFD(MatLay)%TempCond(1:3,2)) >= 0.0d0) THEN ! Multiple Linear Segment Function

      DepVarCol= 2  ! thermal conductivity
      IndVarCol=1  !temperature
          !  Use average  of surface and first node temp for determining k
      kt  = terpld(MaterialFD(MatLay)%numTempCond,MaterialFD(MatLay)%TempCond,(TDT(I)+TDT(I-1))/2.0d0 ,IndVarCol,DepVarCol)

    ENDIF


    RhoS = Material(MatLay)%Density
    Cpo = Material(MatLay)%SpecHeat
    Cp  = Cpo  !  Will be changed if PCM
    Delx = ConstructFD(ConstrNum)%Delx(Lay)

    !Calculate the Dry Heat Conduction Equation

    Rlayer = Material(MatLay)%Resistance
    If (Material(MatLay)%ROnly .or. Material(MatLay)%Group == 1 ) THEN  ! R Layer or Air Layer
            !  Use algebraic equation for TDT based on R

      IF (Surface(Surf)%ExtBoundCond > 0 .and. i==1) THEN  !this is for an adiabatic partition

        TDT(I)=(NetLWRadToSurfFD*Rlayer+QHtRadSysSurfFD*Rlayer + QHWBaseboardSurfFD*Rlayer + QSteamBaseboardSurfFD*Rlayer + &
                    QElecBaseboardSurfFD*Rlayer +  QRadSWInFD*Rlayer + QRadThermInFD*Rlayer +   &
                       TDT(I+1) + hconvi*Rlayer*Tia+TDreport(I)*IterDampConst*Rlayer)/(1.0d0 + hconvi*Rlayer+IterDampConst*Rlayer)


      ELSE ! regular wall
        TDT(I)=(NetLWRadToSurfFD*Rlayer+QHtRadSysSurfFD*Rlayer + QHWBaseboardSurfFD*Rlayer + QSteamBaseboardSurfFD*Rlayer + &
                    QElecBaseboardSurfFD*Rlayer + QRadSWInFD*Rlayer + QRadThermInFD*Rlayer +   &
                       TDT(I-1) + hconvi*Rlayer*Tia+TDreport(I)*IterDampConst*Rlayer)/(1.0d0 + hconvi*Rlayer+IterDampConst*Rlayer)
      ENDIF

      IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
          (TDT(I) < MinSurfaceTempLimit) ) THEN
        TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!        CALL CheckFDSurfaceTempLimits(I,TDT(I))
      ENDIF

    ELSE  !  Regular or PCM


      IF( Sum(MaterialFD(MatLay)%TempEnth(1:3,2)) >= 0.0d0) THEN    !  phase change material,  Use TempEnth Data

        DepVarCol= 2  ! enthalpy
        IndVarCol=1  !temperature

        EnthOld(I) =terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TD(I),IndVarCol,DepVarCol)

        EnthNew(I) =terpld(MaterialFD(MatLay)%numTempEnth,MaterialFD(MatLay)%TempEnth,TDT(I),IndVarCol,DepVarCol)
        IF (ABS(EnthNew(I)-EnthOld(I)) <= smalldiff .or. ABS(TDT(I)-TD(I)) <= smalldiff)  THEN
          Cp = Cpo
        ELSE
          Cp = MAX(Cpo,(EnthNew(I) -EnthOld(I))/(TDT(I)-TD(I)))
        END IF

      ELSE  ! Not phase change material
        Cp= Cpo

      END IF  ! Phase change material check

      IF (Surface(Surf)%ExtBoundCond > 0 .and. i==1) THEN !this is for an adiabatic or interzone partition
        SELECT CASE (CondFDSchemeType)

        CASE (CrankNicholsonSecondOrder)
          ! Adams-Moulton second order
          TDT(I)=(2.d0*Delt*Delx*NetLWRadToSurfFD + 2.d0*Delt*Delx*QHtRadSysSurfFD + 2.d0*Delt*Delx*QHWBaseboardSurfFD + &
                 2.d0*Delt*Delx*QSteamBaseboardSurfFD + 2.d0*Delt*Delx*QElecBaseboardSurfFD + &
                 2.d0*Delt*Delx*QRadSWInFD + 2.d0*Delt*Delx*QRadThermInFD -  &
                  Delt*Delx*hconvi*TD(I) - Delt*kt*TD(I) + Cp*Delx**2*RhoS*TD(I) +  &
                  Delt*kt*TD(I+1) + Delt*kt*TDT(I+1) + 2.d0*Delt*Delx*hconvi*Tia)/  &
                  (Delt*Delx*hconvi + Delt*kt + Cp*Delx**2*RhoS)

        CASE (FullyImplicitFirstOrder)
        ! Adams-Moulton First order
          TDT(I)=(2.d0*Delt*Delx*NetLWRadToSurfFD + 2.d0*Delt*Delx*QHtRadSysSurfFD + 2.d0*Delt*Delx*QHWBaseboardSurfFD + &
                 2.d0*Delt*Delx*QSteamBaseboardSurfFD + 2.d0*Delt*Delx*QElecBaseboardSurfFD + &
                 2.d0*Delt*Delx*QRadSWInFD + 2.d0*Delt*Delx*QRadThermInFD +  &
                  Cp*Delx**2*RhoS*TD(I) + 2.d0*Delt*kt*TDT(I+1) + 2.d0*Delt*Delx*hconvi*Tia)/  &
                  (2.d0*Delt*Delx*hconvi + 2.d0*Delt*kt + Cp*Delx**2*RhoS)
        END SELECT

      ELSE ! for regular or interzone walls
        SELECT CASE (CondFDSchemeType)

        CASE (CrankNicholsonSecondOrder)
          TDT(I)=(2.d0*Delt*Delx*NetLWRadToSurfFD + 2.d0*Delt*Delx*QHtRadSysSurfFD + 2.d0*Delt*Delx*QHWBaseboardSurfFD + &
               2.d0*Delt*Delx*QSteamBaseboardSurfFD + 2.d0*Delt*Delx*QElecBaseboardSurfFD + &
               2.d0*Delt*Delx*QRadSWInFD + 2.d0*Delt*Delx*QRadThermInFD -  &
                Delt*Delx*hconvi*TD(I) - Delt*kt*TD(I) + Cp*Delx**2*RhoS*TD(I) +  &
                Delt*kt*TD(I-1) + Delt*kt*TDT(I-1) + 2.d0*Delt*Delx*hconvi*Tia)/  &
                (Delt*Delx*hconvi + Delt*kt + Cp*Delx**2*RhoS)
        CASE (FullyImplicitFirstOrder)
          TDT(I)=(2.d0*Delt*Delx*NetLWRadToSurfFD + 2.d0*Delt*Delx*QHtRadSysSurfFD + 2.d0*Delt*Delx*QHWBaseboardSurfFD + &
               2.d0*Delt*Delx*QSteamBaseboardSurfFD + 2.d0*Delt*Delx*QElecBaseboardSurfFD + &
               2.d0*Delt*Delx*QRadSWInFD + 2.d0*Delt*Delx*QRadThermInFD +  &
                Cp*Delx**2*RhoS*TD(I) + 2.d0*Delt*kt*TDT(I-1) + 2.d0*Delt*Delx*hconvi*Tia)/  &
                (2.d0*Delt*Delx*hconvi + 2.d0*Delt*kt + Cp*Delx**2*RhoS)
        END SELECT
      ENDIF

      IF ((TDT(I) > MaxSurfaceTempLimit) .OR. &
          (TDT(I) < MinSurfaceTempLimit) ) THEN
        TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check
!        CALL CheckFDSurfaceTempLimits(I,TDT(I))
      ENDIF
!      TDT(I) = Max(MinSurfaceTempLimit,Min(MaxSurfaceTempLimit,TDT(I)))  !  +++++ Limit Check

!  Pass inside conduction Flux [W/m2] to DataHeatBalanceSurface array
!          OpaqSurfInsFaceConductionFlux(Surf)= (TDT(I-1)-TDT(I))*kt/Delx
    END IF  ! Regular or R layer

  END IF   !  End of Regular node or SigmaR SigmaC option


  QNetSurfInside= - (NetLWRadToSurfFD + QHtRadSysSurfFD + QRadSWInFD + QRadThermInFD + QHWBaseboardSurfFD  + &
             QSteamBaseboardSurfFD+QElecBaseboardSurfFD+hconvi*(-TDT(I) + Tia))
! note -- no change ref: CR8575
!feb2012  QNetSurfInside=NetLWRadToSurfFD + QHtRadSysSurfFD + QRadSWInFD + QRadThermInFD + QHWBaseboardSurfFD  + &
!feb2012             QSteamBaseboardSurfFD+QElecBaseboardSurfFD+hconvi*(-TDT(I) + Tia)

    !  Pass inside conduction Flux [W/m2] to DataHeatBalanceSurface array
  OpaqSurfInsFaceConductionFlux(Surf)= QNetSurfInside
!  QFluxZoneToInSurf(Surf) = QNetSurfInside
  OpaqSurfInsFaceConduction(Surf)=QNetSurfInside*Surface(Surf)%Area   !for reporting as in CTF, PT

RETURN

END SUBROUTINE InteriorBCEqns

SUBROUTINE CheckFDSurfaceTempLimits(SurfNum,CheckTemperature)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Linda Lawrie
          !       DATE WRITTEN   August 2012
          !       MODIFIED       na
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! Provides a single entry point for checking surface temperature limits as well as
          ! setting up for recurring errors if too low or too high.

          ! METHODOLOGY EMPLOYED:
          ! Use methodology similar to HBSurfaceManager

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General, ONLY: RoundSigDigits
  USE DataAirFlowNetwork, ONLY: SimulateAirFlowNetwork,AirflowNetworkControlSimple

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

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN)   :: SurfNum  ! surface number
  REAL(R64), INTENT(IN) :: CheckTemperature  ! calculated temperature, not reset

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  INTEGER :: ZoneNum

  ZoneNum=Surface(SurfNum)%Zone

!      IF ((TH(SurfNum,1,2) > MaxSurfaceTempLimit) .OR. &
!          (TH(SurfNum,1,2) < MinSurfaceTempLimit) ) THEN
  IF (WarmupFlag) WarmupSurfTemp=WarmupSurfTemp+1
  IF (.not. WarmupFlag  .or. (WarmupFlag .and. WarmupSurfTemp > 10) .or. DisplayExtraWarnings) THEN
    IF (CheckTemperature < MinSurfaceTempLimit) THEN
      IF (Surface(SurfNum)%LowTempErrCount == 0) THEN
        CALL ShowSevereMessage('Temperature (low) out of bounds ['//TRIM(RoundSigDigits(CheckTemperature,2))//  &
                           '] for zone="'//trim(Zone(ZoneNum)%Name)//'", for surface="'//TRIM(Surface(SurfNum)%Name)//'"')
        CALL ShowContinueErrorTimeStamp(' ')
        IF (.not. Zone(ZoneNum)%TempOutOfBoundsReported) THEN
          CALL ShowContinueError('Zone="'//trim(Zone(ZoneNum)%Name)//'", Diagnostic Details:')
          IF (Zone(ZoneNum)%FloorArea > 0.0d0) THEN
            CALL ShowContinueError('...Internal Heat Gain ['//  &
                 trim(RoundSigDigits(Zone(ZoneNum)%InternalHeatGains/Zone(ZoneNum)%FloorArea,3))//'] W/m2')
          ELSE
            CALL ShowContinueError('...Internal Heat Gain (no floor) ['//  &
                 trim(RoundSigDigits(Zone(ZoneNum)%InternalHeatGains,3))//'] W')
          ENDIF
          IF (SimulateAirflowNetwork <= AirflowNetworkControlSimple) THEN
            CALL ShowContinueError('...Infiltration/Ventilation ['//  &
                   trim(RoundSigDigits(Zone(ZoneNum)%NominalInfilVent,3))//'] m3/s')
            CALL ShowContinueError('...Mixing/Cross Mixing ['//  &
                   trim(RoundSigDigits(Zone(ZoneNum)%NominalMixing,3))//'] m3/s')
          ELSE
            CALL ShowContinueError('...Airflow Network Simulation: Nominal Infiltration/Ventilation/Mixing not available.')
          ENDIF
          IF (Zone(ZoneNum)%isControlled) THEN
            CALL ShowContinueError('...Zone is part of HVAC controlled system.')
          ELSE
            CALL ShowContinueError('...Zone is not part of HVAC controlled system.')
          ENDIF
          Zone(ZoneNum)%TempOutOfBoundsReported=.true.
        ENDIF
        CALL ShowRecurringSevereErrorAtEnd('Temperature (low) out of bounds for zone='//trim(Zone(ZoneNum)%Name)//  &
                       ' for surface='//TRIM(Surface(SurfNum)%Name),  &
                       Surface(SurfNum)%LowTempErrCount,ReportMaxOf=CheckTemperature,ReportMaxUnits='C',  &
                       ReportMinOf=CheckTemperature,ReportMinUnits='C')
      ELSE
        CALL ShowRecurringSevereErrorAtEnd('Temperature (low) out of bounds for zone='//trim(Zone(ZoneNum)%Name)//  &
                       ' for surface='//TRIM(Surface(SurfNum)%Name),  &
                       Surface(SurfNum)%LowTempErrCount,ReportMaxOf=CheckTemperature,ReportMaxUnits='C',  &
                       ReportMinOf=CheckTemperature,ReportMinUnits='C')
      ENDIF
    ELSE
      IF (Surface(SurfNum)%HighTempErrCount == 0) THEN
        CALL ShowSevereMessage('Temperature (high) out of bounds ('//TRIM(RoundSigDigits(CheckTemperature,2))//  &
                           '] for zone="'//trim(Zone(ZoneNum)%Name)//'", for surface="'//TRIM(Surface(SurfNum)%Name)//'"')
        CALL ShowContinueErrorTimeStamp(' ')
        IF (.not. Zone(ZoneNum)%TempOutOfBoundsReported) THEN
          CALL ShowContinueError('Zone="'//trim(Zone(ZoneNum)%Name)//'", Diagnostic Details:')
          IF (Zone(ZoneNum)%FloorArea > 0.0d0) THEN
            CALL ShowContinueError('...Internal Heat Gain ['//  &
                 trim(RoundSigDigits(Zone(ZoneNum)%InternalHeatGains/Zone(ZoneNum)%FloorArea,3))//'] W/m2')
          ELSE
            CALL ShowContinueError('...Internal Heat Gain (no floor) ['//  &
                 trim(RoundSigDigits(Zone(ZoneNum)%InternalHeatGains,3))//'] W')
          ENDIF
          IF (SimulateAirflowNetwork <= AirflowNetworkControlSimple) THEN
            CALL ShowContinueError('...Infiltration/Ventilation ['//  &
                   trim(RoundSigDigits(Zone(ZoneNum)%NominalInfilVent,3))//'] m3/s')
            CALL ShowContinueError('...Mixing/Cross Mixing ['//  &
                   trim(RoundSigDigits(Zone(ZoneNum)%NominalMixing,3))//'] m3/s')
          ELSE
            CALL ShowContinueError('...Airflow Network Simulation: Nominal Infiltration/Ventilation/Mixing not available.')
          ENDIF
          IF (Zone(ZoneNum)%isControlled) THEN
            CALL ShowContinueError('...Zone is part of HVAC controlled system.')
          ELSE
            CALL ShowContinueError('...Zone is not part of HVAC controlled system.')
          ENDIF
          Zone(ZoneNum)%TempOutOfBoundsReported=.true.
        ENDIF
        CALL ShowRecurringSevereErrorAtEnd('Temperature (high) out of bounds for zone='//trim(Zone(ZoneNum)%Name)//  &
                       ' for surface='//TRIM(Surface(SurfNum)%Name),  &
                       Surface(SurfNum)%HighTempErrCount,ReportMaxOf=CheckTemperature,ReportMaxUnits='C',  &
                       ReportMinOf=CheckTemperature,ReportMinUnits='C')
      ELSE
        CALL ShowRecurringSevereErrorAtEnd('Temperature (high) out of bounds for zone='//trim(Zone(ZoneNum)%Name)//  &
                      ' for surface='//TRIM(Surface(SurfNum)%Name),  &
                       Surface(SurfNum)%HighTempErrCount,ReportMaxOf=CheckTemperature,ReportMaxUnits='C',  &
                       ReportMinOf=CheckTemperature,ReportMinUnits='C')
      ENDIF
    ENDIF
  ENDIF


  RETURN

END SUBROUTINE CheckFDSurfaceTempLimits

! *****************************************************************************

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

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