SolveForWindowTemperatures Subroutine

SUBROUTINE SolveForWindowTemperatures(SurfNum)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         F. Winkelmann
          !       DATE WRITTEN   July 2000
          !       MODIFIED       Oct 2000, FW: modify edge-of-glass correction to account
          !                       for gap radiative conductance affects
          !                      Feb 2001, FW: add interior or exterior shade to layer
          !                       heat balance calculation.
          !                      Mar 2001, FW: relax error tolerance if MaxIterations reached.
          !                      Jun 2001, FW: add interior or exterior blind
          !                      Nov 2002, FW: add relaxation on face temperatures
          !                       to improve convergence for multipane cases where outer pane
          !                       has high solar absorptance: temp --> 0.5*(temp + previous temp);
          !                       also, increase MaxIterations from 50 to 100.
          !                      Dec 2002, FW: add between-glass shade/blind for double and triple glazing.
          !                      Mar 2003, FW: remove redundant relaxation on radiative conductances
          !                      Mar 2003, FW: increase convergence tolerance from 0.01 to 0.02 to enhance
          !                                    convergence in difficult cases.
          !                      June 2003, FW: correct the expression for convective gain to zone air
          !                       from airflow windows with airflow destination = InsideAir. Previously
          !                       convective gain of air as it passed through gap was used, which is correct
          !                       for airflow source = InsideAir but incorrect for airflow source = OutsideAir.
          !                       Save SurfaceWindow%TAirflowGapOutlet for use in calculating convective heat
          !                       gain to return air when airflow source = InsideAir, destination = ReturnAir.
          !                      Dec 2003, FW: enhance converge for difficult cases by increasing relaxation
          !                       in layer surface temperatures for iterations > MaxIterations/4
          !                      May 2006, RR: add exterior window screen
          !                      January 2009, BG: inserted call to recalc inside face convection inside iteration loop
          !                        per ISO 15099 Section 8.3.2.2
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! Evaluates the coefficients Aface and Bface in the system of linear
          ! algebraic equations
          !
          !     Sum [Aface(i,j)*thetas(j)] = Bface(i), i = 1,nglfacep, j=1,nglfacep
          !
          ! where
          !
          ! nglface  = number of glass faces (= 2 * number of glass layers), or, if shade or blind is present,
          ! nglgacep = number of glass faces + 2
          !
          ! thetas(j) = temperature of face j
          !
          ! If an interior, exterior or between-glass shade or blind, or exterior screen is present
          ! the face numbering is as follows:
          !   1 to 2*nglface are glass faces, from outside to inside;
          !   2*nglface+1 and 2*nglface+2 are the shade or blind faces, from outside to inside
          ! For example, the following diagram shows the face number for an exterior shade, screen or blind
          ! on double glazing:
          !
          !     ||   ||   ||
          !    5||6 1||2 3||4
          !     ||   ||   ||
          ! bl/sh/sc gl   gl

          ! And for a between-glass shade/blind in triple glazing:
          !
          !     ||   ||   ||   ||
          !    1||2 3||4 7||8 5||6
          !     ||   ||   ||   ||
          !     gl   gl  bl/sh gl

          ! METHODOLOGY EMPLOYED:
          ! The Aface and Bface coefficients are determined by the equations for
          ! heat balance at the glass and shade/blind faces. The system of linear equations is solved
          ! by LU decomposition.

          ! REFERENCES:
          ! na
  USE General, ONLY: InterpSw, InterpSlatAng, TrimSigDigits, RoundSigDigits
  USE Psychrometrics, ONLY:PsyCpAirFnWTdb,PsyRhoAirFnPbTdbW,PsyHFnTdbW,PsyTdbFnHW
  USE InputProcessor, ONLY: SameString
  USE ConvectionCoefficients, ONLY: CalcISO15099WindowIntConvCoeff
  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER           :: SurfNum               ! Surface number

          ! SUBROUTINE PARAMETER DEFINITIONS:
  INTEGER, PARAMETER :: MaxIterations = 100   ! Maximum allowed number of iterations (increased 9/01 from 15 to 50,
                                              !   increased 11/02 from 50 to 100)
  REAL(r64), PARAMETER    :: errtemptol = 0.02d0     ! Tolerance on errtemp for convergence (increased from 0.01, 3/4/03)

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

INTEGER           :: ZoneNum               ! Zone number corresponding to SurfNum
INTEGER           :: i                     ! Counter
REAL(r64)         :: hgap(5)               ! Gap gas conductance (W/m2-K)
REAL(r64)  :: gr                    ! Grashof number of gas in a gap
REAL(r64)  :: con                   ! Gap gas conductivity
REAL(r64)  :: pr                    ! Gap gas Prandtl number
REAL(r64)  :: nu                    ! Gap gas Nusselt number
REAL(r64)         :: hr(10)                ! Radiative conductance (W/m2-K)
REAL(r64)         :: d                     ! +1 if number of row interchanges is even,
                                           ! -1 if odd (in LU decomposition)
INTEGER           :: indx(10)              ! Vector of row permutations in LU decomposition
REAL(r64)  :: Aface(10,10)          ! Coefficient in equation Aface*thetas = Bface
REAL(r64)  :: Bface(10)             ! Coefficient in equation Aface*thetas = Bface

INTEGER           :: iter                  ! Iteration number
REAL(r64)         :: hrprev(10)            ! Value of hr from previous iteration
REAL(r64)         :: errtemp               ! Absolute value of sum of face temperature differences
                                           !   between iterations, divided by number of faces
REAL(r64)         :: VGap                  ! Air velocity in gap between glass and shade/blind (m/s)
REAL(r64)         :: VAirflowGap           ! Air velocity in airflow gap between glass panes (m/s)
REAL(r64)         :: VGapPrev              ! Value of VGap from previous iteration
REAL(r64)         :: TGapNew               ! Average air temp in gap between glass and shade/blind (K)
REAL(r64)         :: TAirFlowGapNew        ! Average air temp in airflow gap between glass panes (K)
REAL(r64)         :: TGapOutlet            ! Temperature of air leaving gap between glass and shade/blind (K)
REAL(r64)         :: TAirflowGapOutlet     ! Temperature of air leaving airflow gap between glass panes (K)
REAL(r64)         :: TAirflowGapOutletC    ! Temperature of air leaving airflow gap between glass panes (C)
REAL(r64)         :: TGapNewBG(2)          ! For between-glass shade/blind, average gas temp in gaps on either
                                           !  side of shade/blind (K)
REAL(r64)         :: hcv                   ! Convection coefficient from gap glass or shade/blind to gap air (W/m2-K)
REAL(r64)         :: hcvAirflowGap         ! Convection coefficient from airflow gap glass to airflow gap air (W/m2-K)
REAL(r64)         :: hcvPrev               ! Value of hcv from previous iteration
REAL(r64)         :: hcvBG(2)              ! For between-glass shade/blind, convection coefficient from gap glass or
                                           !  shade/blind to gap gas on either side of shade/blind (W/m2-K)
REAL(r64)         :: ConvHeatFlowNatural   ! Convective heat flow from gap between glass and interior shade or blind (W)
REAL(r64)         :: ConvHeatFlowForced    ! Convective heat flow from forced airflow gap (W)
REAL(r64)         :: ShGlReflFacIR         ! Factor for long-wave inter-reflection between shade/blind and adjacent glass
REAL(r64)         :: RhoGlIR1,RhoGlIR2     ! Long-wave reflectance of glass surface facing shade/blind; 1=exterior shade/blind,
                                           !  2=interior shade/blind
REAL(r64)         :: RhoShIR1,RhoShIR2     ! Long-wave reflectance of shade/blind surface facing glass; 1=interior shade/blind,
                                           !  2=exterior shade/blind
REAL(r64)         :: EpsShIR1,EpsShIR2     ! Long-wave emissivity of shade/blind surface facing glass; 1=interior shade/blind,
                                           !  2=exterior shade/blind
REAL(r64)         :: TauShIR               ! Long-wave transmittance of isolated shade/blind
REAL(r64)         :: sconsh                ! shade/blind conductance (W/m2-K)
INTEGER           :: ShadeFlag             ! Shading flag
REAL(r64)         :: ShadeAbsFac1,ShadeAbsFac2 ! Fractions for apportioning absorbed radiation to shade/blind faces
REAL(r64)         :: AbsRadShadeFace(2)    ! Solar radiation, short-wave radiation from lights, and long-wave
                                           !  radiation from lights and zone equipment absorbed by faces of shade/blind (W/m2)
REAL(r64)         :: ShadeArea             ! shade/blind area (m2)
REAL(r64)         :: CondHeatGainGlass     ! Conduction through inner glass layer, outside to inside (W)
REAL(r64)         :: CondHeatGainShade     ! Conduction through shade/blind, outside to inside (W)
REAL(r64)         :: NetIRHeatGainGlass    ! Net IR heat gain to zone from shade/blind side of glass when interior
                                           !  shade/blind is present. Zero if shade/blind has zero IR transmittance (W)
REAL(r64)         :: NetIRHeatGainShade    ! Net IR heat gain to zone from interior shade/blind (W)
REAL(r64)         :: ConvHeatGainFrZoneSideOfShade ! Convective heat gain to zone from side of interior shade facing zone (W)
REAL(r64)         :: ConvHeatGainFrZoneSideOfGlass ! Convective heat gain to zone from side of glass facing zone when
                                                   !  no interior shade/blind is present (W)
REAL(r64)         :: IncidentSolar         ! Solar incident on outside of window (W)
INTEGER          :: ConstrNum, ConstrNumSh ! Construction number, bare and with shading device
REAL(r64)         :: TransDiff             ! Diffuse shortwave transmittance
REAL(r64)         :: RhoIR(10)             ! Face IR reflectance
REAL(r64)         :: FacRhoIR25            ! Intermediate variable
REAL(r64)         :: FacRhoIR63            ! Intermediate variable
REAL(r64)         :: RhoIRfp               ! Intermediate variable
REAL(r64)         :: RhoIRbp               ! Intermediate variable
REAL(r64)         :: FacRhoIR2fp           ! Intermediate variable
REAL(r64)         :: FacRhoIR3bp           ! Intermediate variable
REAL(r64)         :: FacRhoIR2fpRhoIR63    ! Intermediate variable
REAL(r64)         :: FacRhoIR3bpRhoIR25    ! Intermediate variable
REAL(r64)         :: FacRhoIR47            ! Intermediate variable
REAL(r64)         :: FacRhoIR85            ! Intermediate variable
REAL(r64)         :: FacRhoIR4fp           ! Intermediate variable
REAL(r64)         :: FacRhoIR5bp           ! Intermediate variable
REAL(r64)         :: FacRhoIR4fpRhoIR85    ! Intermediate variable
REAL(r64)         :: FacRhoIR5bpRhoIR47    ! Intermediate variable
REAL(r64)         :: ConvHeatGainToZoneAir ! Convective heat gain to zone air from window gap airflow (W)
REAL(r64)         :: TotAirflowGap         ! Total volumetric airflow through window gap (m3/s)
REAL(r64)         :: CpAirOutlet           ! Heat capacity of air from window gap (J/kg-K)
REAL(r64)         :: CpAirZone             ! Heat capacity of zone air (J/kg-K)
REAL(r64)         :: InletAirHumRat        ! Humidity ratio of air from window gap entering fan
!!unused REAL(r64)         :: RhoAir                ! Density of air from window gap entering fan (kg/m3)
!!unused REAL(r64)         :: MassFlow              ! Mass flow of air from window gap entering fan (kg/s)
REAL(r64)         :: ZoneTemp              ! Zone air temperature (C)
Integer           :: InsideFaceIndex       ! intermediate variable for index of inside face in thetas

iter = 0
ConvHeatFlowNatural = 0.0d0
ConvHeatFlowForced = 0.0d0
nglfacep = nglface
ShadeFlag = SurfaceWindow(SurfNum)%ShadingFlag
ZoneNum = Surface(SurfNum)%Zone
AbsRadShadeFace = 0.0d0
TGapNew=0.0d0

IF(ShadeFlag==IntShadeOn.OR.ShadeFlag==ExtShadeOn.OR.ShadeFlag==IntBlindOn.OR.ShadeFlag==ExtBlindOn.OR. &
   ShadeFlag==BGShadeOn.OR.ShadeFlag==BGBlindOn.OR.ShadeFlag==ExtScreenOn) THEN
  nglfacep = nglface + 2
  ShadeAbsFac1 = SurfaceWindow(SurfNum)%ShadeAbsFacFace(1)
  ShadeAbsFac2 = SurfaceWindow(SurfNum)%ShadeAbsFacFace(2)
  AbsRadShadeFace(1) = (SurfaceWindow(SurfNum)%ExtBeamAbsByShade + SurfaceWindow(SurfNum)%ExtDiffAbsByShade) * &
                          ShadeAbsFac1 + &
                       (SurfaceWindow(SurfNum)%IntBeamAbsByShade + SurfaceWindow(SurfNum)%IntSWAbsByShade) * &
                          ShadeAbsFac2
  AbsRadShadeFace(2) = (SurfaceWindow(SurfNum)%ExtBeamAbsByShade + SurfaceWindow(SurfNum)%ExtDiffAbsByShade) * &
                          ShadeAbsFac2 + &
                       (SurfaceWindow(SurfNum)%IntBeamAbsByShade + SurfaceWindow(SurfNum)%IntSWAbsByShade) * &
                          ShadeAbsFac1
  IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) &
    AbsRadShadeFace(2) = AbsRadShadeFace(2) + SurfaceWindow(SurfNum)%IntLWAbsByShade
  sconsh  = scon(ngllayer+1)
  TauShIR = tir(nglface+1)
  EpsShIR1 = emis(nglface+1)
  EpsShIR2 = emis(nglface+2)
  RhoShIR1 = MAX(0.d0,1.d0-TauShIR-EpsShIR1)
  RhoShIR2 = MAX(0.d0,1.d0-TauShIR-EpsShIR2)
  IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) THEN
    RhoGlIR2 = 1.d0-emis(2*ngllayer)
    ShGlReflFacIR = 1.d0-RhoGlIR2*RhoShIR1
  ELSE IF(ShadeFlag == ExtShadeOn .OR. ShadeFlag == ExtBlindOn .OR. ShadeFlag == ExtScreenOn) THEN
    RhoGlIR1 = 1.d0-emis(1)
    ShGlReflFacIR = 1.d0-RhoGlIR1*RhoShIR2
  END IF
END IF  ! End of check if shade or blind is on

! Initialize face temperatures.

CALL StartingWindowTemps(SurfNum,AbsRadShadeFace)

hcvPrev=0.0d0
VGapPrev=0.0d0

! Calculate radiative conductances

errtemp=errtemptol*2.0d0

DO WHILE (iter < MaxIterations .AND. errtemp > errtemptol)

  DO i = 1,nglfacep
    hr(i) = emis(i) * sigma * thetas(i)**3
      ! Following line is redundant since thetas is being relaxed;
      ! removed by FCW, 3/4/03
      !!fw if(iter >= 1) hr(i) = 0.5*(hrprev(i)+hr(i))
    hrprev(i) = hr(i)
  END DO

! call for new interior film coeff (since it is temperature dependent) if using Detailed inside coef model
  IF (((Surface(SurfNum)%IntConvCoeff == 0) .AND. (Zone(ZoneNum)%InsideConvectionAlgo == ASHRAETARP)) &
     .OR. (Surface(SurfNum)%IntConvCoeff == -2)) Then
    ! coef model is "detailed" and not prescribed by user
    !need to find inside face index, varies with shade/blind etc.
    IF (ShadeFlag==IntShadeOn .OR. ShadeFlag==IntBlindOn) Then
      InsideFaceIndex = nglfacep
    ELSE
      InsideFaceIndex = nglface
    ENDIF
    CALL CalcISO15099WindowIntConvCoeff(SurfNum,thetas(InsideFaceIndex)-KelvinConv,tin-KelvinConv)
    hcin = HconvIn(SurfNum)
  ENDIF

  Aface = 0.0d0
  Bface = 0.0d0

! If interior or exterior shade or blind is present, get heat transfer
! coefficient from glass and shade/blind to gap between glass and shade/blind,
! effective gap air temperature, velocity of air in gap and gap outlet temperature.

  IF(ShadeFlag==IntShadeOn.OR.ShadeFlag==ExtShadeOn.OR.ShadeFlag==IntBlindOn.OR.ShadeFlag==ExtBlindOn.OR. &
     ShadeFlag==ExtScreenOn) THEN
    CALL ExtOrIntShadeNaturalFlow(SurfNum,iter,VGap,TGapNew,TGapOutlet,hcv,ConvHeatFlowNatural)
    IF(iter >= 1) THEN
      hcv =  0.5d0*(hcvPrev + hcv)
      VGap = 0.5d0*(VGapPrev + VGap)
    END IF
    hcvPrev = hcv
    VGapPrev = VGap
  END IF

  TAirFlowGapOutlet=0.0d0
! If between-glass shade or blind is not present and this is an airflow window
! (i.e., with forced airflow in the gap for double glass or in the inner gap for triple glass)
! get glass-to-air forced convection heat transfer coefficient, average gap air temperature, and
! convective heat flow from gap.

  IF(ShadeFlag /= BGShadeOn .AND. ShadeFlag /= BGBlindOn .AND. SurfaceWindow(SurfNum)%AirflowThisTS > 0.0d0) THEN
    CALL BetweenGlassForcedFlow(SurfNum,iter,VAirflowGap,TAirFlowGapNew,TAirFlowGapOutlet,hcvAirflowGap,ConvHeatFlowForced)
  ENDIF

! If between-glass shade or blind is present, get convective heat transfer
! coefficients from glass and shade/blind to the two gaps on either side of the shade/blind.
! Also get average gas temperature in the two gaps, and, for airflow window, the sum of the
! convective heat flows from the gaps.

  IF(ShadeFlag == BGShadeOn .OR. ShadeFlag == BGBlindOn) THEN
    IF(SurfaceWindow(SurfNum)%AirflowThisTS == 0.0d0) THEN  ! Natural convection in gaps
      CALL BetweenGlassShadeNaturalFlow(SurfNum,iter,VGap,TGapNewBG,hcvBG)
    ELSE                                            ! Forced convection in gaps
      CALL BetweenGlassShadeForcedFlow(SurfNum,iter,VGap,TGapNewBG,TAirFlowGapOutlet,hcvBG,ConvHeatFlowForced)
    END IF
  END IF

  iter = iter + 1
  SurfaceWindow(SurfNum)%WindowCalcIterationsRep = iter

  ! Calculations based on number of glass layers
  SELECT CASE(ngllayer)

    CASE (1)
      Bface(1) = outir*emis(1) + hcout*tout + AbsRadGlassFace(1)
      Bface(2) = rmir*emis(2) + hcin*tin + AbsRadGlassFace(2)

      Aface(1,1) = hr(1) + scon(1) + hcout
      Aface(1,2) = -scon(1)
      Aface(2,1) = -scon(1)
      Aface(2,2) = hr(2) + scon(1) + hcin

      IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) THEN
        Bface(2) = rmir*emis(2)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(2)
        Bface(3) = rmir*TauShIR*RhoGlIR2*EpsShIR1/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(1)
        Bface(4) = rmir*EpsShIR2 + hcin*tin + AbsRadShadeFace(2)

        Aface(2,2) = hr(2)*(1-RhoShIR1)/ShGlReflFacIR + scon(1) + hcv
        Aface(2,3) = -emis(2)*hr(3)/ShGlReflFacIR
        Aface(3,2) = -hr(2)*EpsShIR1/ShGlReflFacIR
        Aface(3,3) = hr(3)*(1-RhoGlIR2*(EpsShIR1+RhoShIR1))/ShGlReflFacIR + sconsh + hcv
        Aface(3,4) = -sconsh
        Aface(4,3) = -sconsh
        Aface(4,4) = hr(4) + sconsh + hcin
      END IF

      IF(ShadeFlag == ExtShadeOn .OR. ShadeFlag == ExtBlindOn .OR. ShadeFlag == ExtScreenOn) THEN
        Bface(1) = outir*emis(1)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(1)
        Bface(3) = outir*EpsShIR1 + hcout*tout + AbsRadShadeFace(1)
        Bface(4) = outir*TauShIR*RhoGlIR1*EpsShIR2/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(2)

        Aface(1,1) = hr(1)*(1-RhoShIR2)/ShGlReflFacIR + scon(1) + hcv
        Aface(1,4) = -emis(1)*hr(4)/ShGlReflFacIR
        Aface(3,3) = hr(3) + sconsh + hcout
        Aface(3,4) = -sconsh
        Aface(4,1) = -hr(1)*EpsShIR2/ShGlReflFacIR
        Aface(4,3) = -sconsh
        Aface(4,4) = hr(4)*(1-RhoGlIR1*(EpsShIR2+RhoShIR2))/ShGlReflFacIR + sconsh + hcv
      END IF

    CASE (2)
      call WindowGasConductance(thetas(2),thetas(3),1,con,pr,gr)
      call NusseltNumber(SurfNum,thetas(2),thetas(3),1,gr,pr,nu)
      hgap(1) = con/gap(1)*nu
      IF(SurfaceWindow(SurfNum)%EdgeGlCorrFac > 1.0d0) THEN ! Edge of glass correction
        hrgap(1) = 0.5d0*ABS(A23)*(thetas(2)+thetas(3))**3
        hgap(1) = hgap(1) * SurfaceWindow(SurfNum)%EdgeGlCorrFac + &
                  hrgap(1) * (SurfaceWindow(SurfNum)%EdgeGlCorrFac - 1.0d0)
      END IF

      Bface(1) = outir*emis(1) + hcout*tout + AbsRadGlassFace(1)
      Bface(2) = AbsRadGlassFace(2)
      Bface(3) = AbsRadGlassFace(3)
      Bface(4) = rmir*emis(4) + hcin*tin + AbsRadGlassFace(4)

      Aface(1,1) = hr(1) + scon(1) + hcout
      Aface(1,2) = -scon(1)

      Aface(2,1) = -scon(1)
      Aface(2,2) = scon(1) + hgap(1) - A23P*hr(2)
      Aface(2,3) = -hgap(1) - A32P*hr(3)

      Aface(3,2) = -hgap(1) + A23P*hr(2)
      Aface(3,3) = hgap(1) + scon(2) + A32P*hr(3)
      Aface(3,4) = -scon(2)

      Aface(4,3) = -scon(2)
      Aface(4,4) = hr(4) + scon(2) + hcin

      IF(ShadeFlag /= BGShadeOn .AND. ShadeFlag /= BGBlindOn .AND. SurfaceWindow(SurfNum)%AirflowThisTS > 0.0d0) THEN
        Bface(2) = AbsRadGlassFace(2) + hcvAirflowGap*TAirflowGapNew
        Bface(3) = AbsRadGlassFace(3) + hcvAirflowGap*TAirflowGapNew
        Aface(2,2) = scon(1) + hcvAirflowGap - A23P*hr(2)
        Aface(2,3) = -A32P*hr(3)
        Aface(3,2) = A23P*hr(2)
        Aface(3,3) = hcvAirflowGap + scon(2) + A32P*hr(3)
      END IF

      IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) THEN
        Bface(4) = rmir*emis(4)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(4)
        Bface(5) = rmir*TauShIR*RhoGlIR2*EpsShIR1/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(1)
        Bface(6) = rmir*EpsShIR2 + hcin*tin + AbsRadShadeFace(2)

        Aface(4,4) = hr(4)*(1-RhoShIR1)/ShGlReflFacIR + scon(2) + hcv
        Aface(4,5) = -emis(4)*hr(5)/ShGlReflFacIR
        Aface(5,4) = -hr(4)*EpsShIR1/ShGlReflFacIR
        Aface(5,5) = hr(5)*(1-RhoGlIR2*(EpsShIR1+RhoShIR1))/ShGlReflFacIR + sconsh + hcv
        Aface(5,6) = -sconsh
        Aface(6,5) = -sconsh
        Aface(6,6) = hr(6) + sconsh + hcin
      END IF

      IF(ShadeFlag == ExtShadeOn .OR. ShadeFlag == ExtBlindOn .OR. ShadeFlag == ExtScreenOn) THEN
        Bface(1) = outir*emis(1)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(1)
        Bface(5) = outir*EpsShIR1 + hcout*tout + AbsRadShadeFace(1)
        Bface(6) = outir*TauShIR*RhoGlIR1*EpsShIR2/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(2)

        Aface(1,1) = hr(1)*(1-RhoShIR2)/ShGlReflFacIR + scon(1) + hcv
        Aface(1,6) = -emis(1)*hr(6)/ShGlReflFacIR
        Aface(5,5) = hr(5) + sconsh + hcout
        Aface(5,6) = -sconsh
        Aface(6,1) = -hr(1)*EpsShIR2/ShGlReflFacIR
        Aface(6,5) = -sconsh
        Aface(6,6) = hr(6)*(1-RhoGlIR1*(EpsShIR2+RhoShIR2))/ShGlReflFacIR + sconsh + hcv
      END IF

      IF(ShadeFlag == BGShadeOn .OR. ShadeFlag == BGBlindOn) THEN
        DO i = 1,6
          RhoIR(i) = MAX(0.0d0,1.d0-tir(i)-emis(i))
        END DO
        FacRhoIR25  = 1.d0-RhoIr(2)*RhoIR(5)
        FacRhoIR63  = 1.d0-RhoIr(6)*RhoIR(3)
        RhoIRfp     = RhoIR(5) + (tir(5)**2)*RhoIR(3)/FacRhoIR63
        RhoIRbp     = RhoIR(6) + (tir(5)**2)*RhoIR(2)/FacRhoIR25
        FacRhoIR2fp = 1.d0-RhoIRfp*RhoIR(2)
        FacRhoIR3bp = 1.d0-RhoIRbp*RhoIR(3)
        FacRhoIR2fpRhoIR63 = FacRhoIR2fp * FacRhoIR63
        FacRhoIR3bpRhoIR25 = FacRhoIR3bp * FacRhoIR25
        Aface(2,2)  = scon(1) + hcvBG(1) + hr(2)*(1-RhoIRfp*(emis(2)+RhoIR(2)))/FacRhoIR2fp
        Aface(2,3)  = -emis(2)*hr(3)*tir(5)/FacRhoIR2fpRhoIR63
        Aface(2,5)  = -emis(2)*hr(5)/FacRhoIR2fp
        Aface(2,6)  = -emis(2)*hr(6)*RhoIR(3)*tir(5)/FacRhoIR2fpRhoIR63
          Bface(2)  = hcvBG(1)*TGapNewBG(1) + AbsRadGlassFace(2)
        Aface(3,2)  = -emis(3)*hr(2)*tir(5)/FacRhoIR3bpRhoIR25
        Aface(3,3)  = scon(2) + hcvBG(2) + hr(3)*(1-RhoIRbp*(emis(3)+RhoIR(3)))/FacRhoIR3bp
        Aface(3,5)  = -emis(3)*hr(5)*RhoIR(2)*tir(5)/FacRhoIR3bpRhoIR25
        Aface(3,6)  = -emis(3)*hr(6)/FacRhoIR3bp
          Bface(3)  = hcvBG(2)*TGapNewBG(2) + AbsRadGlassFace(3)
        Aface(5,2)  = -emis(5)*hr(2)/FacRhoIR2fp
        Aface(5,3)  = -hr(3)*tir(5)*RhoIR(2)*emis(5)/FacRhoIR2fpRhoIR63
        Aface(5,5)  = sconsh + hcvBG(1) + hr(5)*(1-RhoIR(2)*emis(5)/FacRhoIr2fp)
        Aface(5,6)  = -sconsh - hr(6)*RhoIR(2)*tir(5)*RhoIR(3)*emis(5)/FacRhoIR2fpRhoIR63
          Bface(5)  = hcvBG(1)*TGapNewBG(1) + AbsRadShadeFace(1)
        Aface(6,2)  = -hr(2)*tir(5)*RhoIR(3)*emis(6)/FacRhoIR3bpRhoIR25
        Aface(6,3)  = -emis(6)*hr(3)/FacRhoIR3bp
        Aface(6,5)  = -sconsh - hr(5)*RhoIR(3)*tir(5)*RhoIR(2)*emis(6)/FacRhoIR3bpRhoIR25
        Aface(6,6)  = sconsh + hcvBG(2) + hr(6)*(1-RhoIR(3)*emis(6)/FacRhoIR3bp)
          Bface(6)  = hcvBG(2)*TGapNewBG(2) + AbsRadShadeFace(2)
      END IF

    CASE (3)
      call WindowGasConductance(thetas(2),thetas(3),1,con,pr,gr)
      call NusseltNumber(SurfNum,thetas(2),thetas(3),1,gr,pr,nu)
      hgap(1) = con/gap(1)*nu
      IF(SurfaceWindow(SurfNum)%EdgeGlCorrFac > 1.0d0) THEN ! Edge of glass correction
        hrgap(1) = 0.5d0*ABS(A23)*(thetas(2)+thetas(3))**3
        hgap(1) = hgap(1) * SurfaceWindow(SurfNum)%EdgeGlCorrFac + &
                  hrgap(1) * (SurfaceWindow(SurfNum)%EdgeGlCorrFac - 1.0d0)
      END IF

      call WindowGasConductance(thetas(4),thetas(5),2,con,pr,gr)
      call NusseltNumber(SurfNum,thetas(4),thetas(5),2,gr,pr,nu)
      hgap(2) = con/gap(2)*nu
      IF(SurfaceWindow(SurfNum)%EdgeGlCorrFac > 1.0d0) THEN ! Edge of glass correction
        hrgap(2) = 0.5d0*ABS(A45)*(thetas(4)+thetas(5))**3
        hgap(2) = hgap(2) * SurfaceWindow(SurfNum)%EdgeGlCorrFac + &
                  hrgap(2) * (SurfaceWindow(SurfNum)%EdgeGlCorrFac - 1.0d0)
      END IF

      Bface(1) = outir*emis(1) + hcout*tout + AbsRadGlassFace(1)
      Bface(2) = AbsRadGlassFace(2)
      Bface(3) = AbsRadGlassFace(3)
      Bface(4) = AbsRadGlassFace(4)
      Bface(5) = AbsRadGlassFace(5)
      Bface(6) = rmir*emis(6) + hcin*tin + AbsRadGlassFace(6)

      Aface(1,1) = hr(1) + scon(1) + hcout
      Aface(1,2) = -scon(1)

      Aface(2,1) = -scon(1)
      Aface(2,2) = scon(1) + hgap(1) - A23P*hr(2)
      Aface(2,3) = -hgap(1) - A32P*hr(3)

      Aface(3,2) = -hgap(1) + A23P*hr(2)
      Aface(3,3) = hgap(1) + scon(2) + A32P*hr(3)
      Aface(3,4) = -scon(2)

      Aface(4,3) = -scon(2)
      Aface(4,4) = scon(2) + hgap(2) - A45P*hr(4)
      Aface(4,5) = -hgap(2) - A54P*hr(5)

      Aface(5,4) = -hgap(2) + A45P*hr(4)
      Aface(5,5) = hgap(2) + scon(3) + A54P*hr(5)
      Aface(5,6) = -scon(3)

      Aface(6,5) = -scon(3)
      Aface(6,6) = hr(6) + scon(3) + hcin

      IF(ShadeFlag /= BGShadeOn .AND. ShadeFlag /= BGBlindOn .AND. SurfaceWindow(SurfNum)%AirflowThisTS > 0.0d0) THEN
        Bface(4) = AbsRadGlassFace(4) + hcvAirflowGap*TAirflowGapNew
        Bface(5) = AbsRadGlassFace(5) + hcvAirflowGap*TAirflowGapNew
        Aface(4,4) = scon(2) + hcvAirflowGap - A45P*hr(4)
        Aface(4,5) = -A54P*hr(5)
        Aface(5,4) = A45P*hr(4)
        Aface(5,5) = hcvAirflowGap + scon(3) + A54P*hr(5)
      END IF

      IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) THEN
        Bface(6) = rmir*emis(6)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(6)
        Bface(7) = rmir*TauShIR*RhoGlIR2*EpsShIR1/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(1)
        Bface(8) = rmir*EpsShIR2 + hcin*tin + AbsRadShadeFace(2)

        Aface(6,6) = hr(6)*(1-RhoShIR1)/ShGlReflFacIR + scon(3) + hcv
        Aface(6,7) = -emis(6)*hr(7)/ShGlReflFacIR
        Aface(7,6) = -hr(6)*EpsShIR1/ShGlReflFacIR
        Aface(7,7) = hr(7)*(1-RhoGlIR2*(EpsShIR1+RhoShIR1))/ShGlReflFacIR + sconsh + hcv
        Aface(7,8) = -sconsh
        Aface(8,7) = -sconsh
        Aface(8,8) = hr(8) + sconsh + hcin
      END IF

      IF(ShadeFlag == ExtShadeOn .OR. ShadeFlag == ExtBlindOn .OR. ShadeFlag == ExtScreenOn) THEN
        Bface(1) = outir*emis(1)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(1)
        Bface(7) = outir*EpsShIR1 + hcout*tout + AbsRadShadeFace(1)
        Bface(8) = outir*TauShIR*RhoGlIR1*EpsShIR2/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(2)

        Aface(1,1) = hr(1)*(1-RhoShIR2)/ShGlReflFacIR + scon(1) + hcv
        Aface(1,8) = -emis(1)*hr(8)/ShGlReflFacIR
        Aface(7,7) = hr(7) + sconsh + hcout
        Aface(7,8) = -sconsh
        Aface(8,1) = -hr(1)*EpsShIR2/ShGlReflFacIR
        Aface(8,7) = -sconsh
        Aface(8,8) = hr(8)*(1-RhoGlIR1*(EpsShIR2+RhoShIR2))/ShGlReflFacIR + sconsh + hcv
      END IF

      IF(ShadeFlag == BGShadeOn .OR. ShadeFlag == BGBlindOn) THEN
        DO i = 1,8
          RhoIR(i) = MAX(0.0d0,1.d0-tir(i)-emis(i))
        END DO
        FacRhoIR47  = 1-RhoIr(4)*RhoIR(7)
        FacRhoIR85  = 1-RhoIr(8)*RhoIR(5)
        RhoIRfp     = RhoIR(7) + (tir(7)**2)*RhoIR(5)/FacRhoIR85
        RhoIRbp     = RhoIR(8) + (tir(7)**2)*RhoIR(4)/FacRhoIR47
        FacRhoIR4fp = 1-RhoIRfp*RhoIR(4)
        FacRhoIR5bp = 1-RhoIRbp*RhoIR(5)
        FacRhoIR4fpRhoIR85 = FacRhoIR4fp * FacRhoIR85
        FacRhoIR5bpRhoIR47 = FacRhoIR5bp * FacRhoIR47
        Aface(4,4)  = scon(2) + hcvBG(1) + hr(4)*(1-RhoIRfp*(emis(4)+RhoIR(4)))/FacRhoIR4fp
        Aface(4,5)  = -emis(4)*hr(5)*tir(7)/FacRhoIR4fpRhoIR85
        Aface(4,7)  = -emis(4)*hr(7)/FacRhoIR4fp
        Aface(4,8)  = -emis(4)*hr(8)*RhoIR(5)*tir(7)/FacRhoIR4fpRhoIR85
          Bface(4)  = hcvBG(1)*TGapNewBG(1) + AbsRadGlassFace(4)
        Aface(5,4)  = -emis(5)*hr(4)*tir(7)/FacRhoIR5bpRhoIR47
        Aface(5,5)  = scon(3) + hcvBG(2) + hr(5)*(1-RhoIRbp*(emis(5)+RhoIR(5)))/FacRhoIR5bp
        Aface(5,7)  = -emis(5)*hr(7)*RhoIR(4)*tir(7)/FacRhoIR5bpRhoIR47
        Aface(5,8)  = -emis(5)*hr(8)/FacRhoIR5bp
          Bface(5)  = hcvBG(2)*TGapNewBG(2) + AbsRadGlassFace(5)
        Aface(7,4)  = -emis(7)*hr(4)/FacRhoIR4fp
        Aface(7,5)  = -hr(5)*tir(7)*RhoIR(4)*emis(7)/FacRhoIR4fpRhoIR85
        Aface(7,7)  = sconsh + hcvBG(1) + hr(7)*(1-RhoIR(4)*emis(7)/FacRhoIr4fp)
        Aface(7,8)  = -sconsh - hr(8)*RhoIR(4)*tir(7)*RhoIR(5)*emis(7)/FacRhoIR4fpRhoIR85
          Bface(7)  = hcvBG(1)*TGapNewBG(1) + AbsRadShadeFace(1)
        Aface(8,4)  = -hr(4)*tir(7)*RhoIR(5)*emis(8)/FacRhoIR5bpRhoIR47
        Aface(8,5)  = -emis(8)*hr(5)/FacRhoIR5bp
        Aface(8,7)  = -sconsh - hr(7)*RhoIR(5)*tir(7)*RhoIR(4)*emis(8)/FacRhoIR5bpRhoIR47
        Aface(8,8)  = sconsh + hcvBG(2) + hr(8)*(1-RhoIR(5)*emis(8)/FacRhoIR5bp)
          Bface(8)  = hcvBG(2)*TGapNewBG(2) + AbsRadShadeFace(2)
      END IF

    CASE (4)
      call WindowGasConductance(thetas(2),thetas(3),1,con,pr,gr)
      call NusseltNumber(SurfNum,thetas(2),thetas(3),1,gr,pr,nu)
      hgap(1) = con/gap(1)*nu
      IF(SurfaceWindow(SurfNum)%EdgeGlCorrFac > 1.0d0) THEN ! Edge of glass correction
        hrgap(1) = 0.5d0*ABS(A23)*(thetas(2)+thetas(3))**3
        hgap(1) = hgap(1) * SurfaceWindow(SurfNum)%EdgeGlCorrFac + &
                  hrgap(1) * (SurfaceWindow(SurfNum)%EdgeGlCorrFac - 1.0d0)
      END IF

      call WindowGasConductance(thetas(4),thetas(5),2,con,pr,gr)
      call NusseltNumber(SurfNum,thetas(4),thetas(5),2,gr,pr,nu)
      hgap(2) = con/gap(2)*nu
      IF(SurfaceWindow(SurfNum)%EdgeGlCorrFac > 1.0d0) THEN ! Edge of glass correction
        hrgap(2) = 0.5d0*ABS(A45)*(thetas(4)+thetas(5))**3
        hgap(2) = hgap(2) * SurfaceWindow(SurfNum)%EdgeGlCorrFac + &
                  hrgap(2) * (SurfaceWindow(SurfNum)%EdgeGlCorrFac - 1.0d0)
      END IF

      call WindowGasConductance(thetas(6),thetas(7),3,con,pr,gr)
      call NusseltNumber(SurfNum,thetas(6),thetas(7),3,gr,pr,nu)
      hgap(3) = con/gap(3)*nu
      IF(SurfaceWindow(SurfNum)%EdgeGlCorrFac > 1.0d0) THEN ! Edge of glass correction
        hrgap(3) = 0.5d0*ABS(A67)*(thetas(6)+thetas(7))**3
        hgap(3) = hgap(3) * SurfaceWindow(SurfNum)%EdgeGlCorrFac + &
                  hrgap(3) * (SurfaceWindow(SurfNum)%EdgeGlCorrFac - 1.0d0)
      END IF
      Bface(1) = outir*emis(1) + hcout*tout + AbsRadGlassFace(1)
      Bface(2) = AbsRadGlassFace(2)
      Bface(3) = AbsRadGlassFace(3)
      Bface(4) = AbsRadGlassFace(4)
      Bface(5) = AbsRadGlassFace(5)
      Bface(6) = AbsRadGlassFace(6)
      Bface(7) = AbsRadGlassFace(7)
      Bface(8) = rmir*emis(8) + hcin*tin + AbsRadGlassFace(8)

      Aface(1,1) = hr(1) + scon(1) + hcout
      Aface(1,2) = -scon(1)

      Aface(2,1) = -scon(1)
      Aface(2,2) = scon(1) + hgap(1) - A23P*hr(2)
      Aface(2,3) = -hgap(1) - A32P*hr(3)

      Aface(3,2) = -hgap(1) + A23P*hr(2)
      Aface(3,3) = hgap(1) + scon(2) + A32P*hr(3)
      Aface(3,4) = -scon(2)

      Aface(4,3) = -scon(2)
      Aface(4,4) = scon(2) + hgap(2) - A45P*hr(4)
      Aface(4,5) = -hgap(2) - A54P*hr(5)

      Aface(5,4) = -hgap(2) + A45P*hr(4)
      Aface(5,5) = hgap(2) + scon(3) + A54P*hr(5)
      Aface(5,6) = -scon(3)

      Aface(6,5) = -scon(3)
      Aface(6,6) = scon(3) + hgap(3) - A67P*hr(6)
      Aface(6,7) = -hgap(3) - A76P*hr(7)

      Aface(7,6) = -hgap(3) + A67P*hr(6)
      Aface(7,7) = hgap(3) + scon(4) + A76P*hr(7)
      Aface(7,8) = -scon(4)

      Aface(8,7) = -scon(4)
      Aface(8,8) = hr(8) + scon(4) + hcin

      IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) THEN
        Bface(8) = rmir*emis(8)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(8)
        Bface(9) = rmir*TauShIR*RhoGlIR2*EpsShIR1/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(1)
        Bface(10) = rmir*EpsShIR2 + hcin*tin + AbsRadShadeFace(2)

        Aface(8,8) = hr(8)*(1-RhoShIR1)/ShGlReflFacIR + scon(4) + hcv
        Aface(8,9) = -emis(8)*hr(9)/ShGlReflFacIR
        Aface(9,8) = -hr(8)*EpsShIR1/ShGlReflFacIR
        Aface(9,9) = hr(9)*(1-RhoGlIR2*(EpsShIR1+RhoShIR1))/ShGlReflFacIR + sconsh + hcv
        Aface(9,10) = -sconsh
        Aface(10,9) = -sconsh
        Aface(10,10) = hr(10) + sconsh + hcin
      END IF

      IF(ShadeFlag == ExtShadeOn .OR. ShadeFlag == ExtBlindOn .OR. ShadeFlag == ExtScreenOn) THEN
        Bface(1) = outir*emis(1)*TauShIR/ShGlReflFacIR + hcv*TGapNew + AbsRadGlassFace(1)
        Bface(9) = outir*EpsShIR1 + hcout*tout + AbsRadShadeFace(1)
        Bface(10) = outir*TauShIR*RhoGlIR1*EpsShIR2/ShGlReflFacIR + hcv*TGapNew + AbsRadShadeFace(2)

        Aface(1,1) = hr(1)*(1-RhoShIR2)/ShGlReflFacIR + scon(1) + hcv
        Aface(1,10) = -emis(1)*hr(10)/ShGlReflFacIR
        Aface(9,9) = hr(9) + sconsh + hcout
        Aface(9,10) = -sconsh
        Aface(10,1) = -hr(1)*EpsShIR2/ShGlReflFacIR
        Aface(10,9) = -sconsh
        Aface(10,10) = hr(10)*(1-RhoGlIR1*(EpsShIR2+RhoShIR2))/ShGlReflFacIR + sconsh + hcv
      END IF

    CASE DEFAULT
      CALL ShowFatalError('SolveForWindowTemperatures: Invalid number of Glass Layers='//  &
                                    TRIM(TrimSigDigits(ngllayer))//', up to 4 allowed.')
  END SELECT

  call LUdecomposition(Aface,nglfacep,indx,d)    ! Note that these routines change Aface;
  call LUsolution(Aface,nglfacep,indx,Bface)     ! face temperatures are returned in Bface

  DO i = 1,nglfacep
    thetasPrev(i) = thetas(i)
    IF(iter < MaxIterations/4) THEN
      thetas(i) = 0.5d0*thetas(i) + 0.5d0*Bface(i)
    ELSE
      thetas(i) = 0.75d0*thetas(i) + 0.25d0*Bface(i)
    END IF
  END DO

  errtemp = 0.0d0
  DO i = 1,nglfacep
    errtemp = errtemp + ABS(thetas(i)-thetasPrev(i))
  END DO
  errtemp = errtemp/nglfacep

END DO

       ! We have reached iteration limit or we have converged. If we have reached the
       ! iteration limit the following test relaxes the convergence tolerance.
       ! If we have converged (errtemp <= errtemptol) the following test has not effect.

IF(errtemp < 10*errtemptol) THEN

       ! Window heat balance solution has converged.

       ! For interior shade, add convective gain from glass/shade gap air flow to zone convective gain;
       ! For all cases, get total window heat gain for reporting. See CalcWinFrameAndDividerTemps for
       ! contribution of frame and divider.
  IncidentSolar = Surface(SurfNum)%Area * QRadSWOutIncident(SurfNum)
  IF(ShadeFlag == IntShadeOn .OR. ShadeFlag == IntBlindOn) THEN
       ! Interior shade or blind
    SurfaceWindow(SurfNum)%ConvHeatFlowNatural = ConvHeatFlowNatural
       ! Window heat gain from glazing and shade/blind to zone. Consists of transmitted solar, convection
       !   from air exiting gap, convection from zone-side of shade/blind, net IR to zone from shade and net IR to
       !   zone from the glass adjacent to the shade/blind (zero if shade/blind IR transmittance is zero).
       ! Following assumes glazed area = window area (i.e., dividers ignored) in calculating
       !   IR to zone from glass when interior shade/blind is present.
    ShadeArea = Surface(SurfNum)%Area + SurfaceWindow(SurfNum)%DividerArea
    CondHeatGainShade  = ShadeArea * sconsh * (thetas(nglfacep-1) - thetas(nglfacep))
    NetIRHeatGainShade = ShadeArea * &
    EpsShIR2*(sigma*thetas(nglfacep)**4 - rmir) + &
    EpsShIR1*(sigma*thetas(nglfacep-1)**4 - rmir)*RhoGlIR2*TauShIR/ShGlReflFacIR
    NetIRHeatGainGlass = ShadeArea * &
       (emis(2*ngllayer)*TauShIR/ShGlReflFacIR) * (sigma*thetas(2*ngllayer)**4 - rmir)
    ConvHeatGainFrZoneSideOfShade = ShadeArea * hcin*(thetas(nglfacep) - tin)
    WinHeatGain(SurfNum) = WinTransSolar(SurfNum) + ConvHeatFlowNatural + ConvHeatGainFrZoneSideOfShade + &
                         NetIRHeatGainGlass + NetIRHeatGainShade
    ! store components for reporting
    WinGainConvGlazShadGapToZoneRep(SurfNum) = ConvHeatFlowNatural
    WinGainConvShadeToZoneRep(SurfNum)       = ConvHeatGainFrZoneSideOfShade
    WinGainIRGlazToZoneRep(SurfNum)          = NetIRHeatGainGlass
    WinGainIRShadeToZoneRep(SurfNum)         =  NetIRHeatGainShade
  ELSE
       ! Interior shade or blind not present; innermost layer is glass
    CondHeatGainGlass = Surface(SurfNum)%Area * scon(ngllayer) * (thetas(2*ngllayer-1)-thetas(2*ngllayer))
    NetIRHeatGainGlass = Surface(SurfNum)%Area * emis(2*ngllayer)*(sigma*thetas(2*ngllayer)**4 - rmir)
    ConvHeatGainFrZoneSideOfGlass = Surface(SurfNum)%Area * hcin*(thetas(2*ngllayer) - tin)
    WinHeatGain(SurfNum) = WinTransSolar(SurfNum) + ConvHeatGainFrZoneSideOfGlass + NetIRHeatGainGlass
    ! store components for reporting
    WinGainConvGlazToZoneRep(SurfNum) = ConvHeatGainFrZoneSideOfGlass
    WinGainIRGlazToZoneRep(SurfNum) = NetIRHeatGainGlass
  END IF

      ! Add convective heat gain from airflow window
      ! Note: effect of fan heat on gap outlet temperature is neglected since fan power (based
      ! on pressure drop through the gap) is extremely small

  WinGapConvHtFlowRep(SurfNum) = 0.0d0
  WinGapConvHtFlowRepEnergy(SurfNum) = 0.0d0
  TotAirflowGap = SurfaceWindow(SurfNum)%AirFlowThisTS * Surface(SurfNum)%Width
  TAirflowGapOutletC = TAirflowGapOutlet-TKelvin
  SurfaceWindow(SurfNum)%TAirflowGapOutlet = TAirflowGapOutletC
  IF(SurfaceWindow(SurfNum)%AirFlowThisTS > 0.0d0) THEN
    WinGapConvHtFlowRep(SurfNum) = ConvHeatFlowForced
    WinGapConvHtFlowRepEnergy(SurfNum) = WinGapConvHtFlowRep(SurfNum) * TimeStepZone * SecInHour
    ! Add heat from gap airflow to zone air if destination is inside air; save the heat gain to return
    ! air in case it needs to be sent to the zone (due to no return air determined in HVAC simulation)
    IF(SurfaceWindow(SurfNum)%AirFlowDestination == AirFlowWindow_Destination_IndoorAir  .or.   &
       SurfaceWindow(SurfNum)%AirFlowDestination == AirFlowWindow_Destination_ReturnAir) THEN
      IF (SurfaceWindow(SurfNum)%AirflowSource == AirFlowWindow_Source_IndoorAir) THEN
        InletAirHumRat = ZoneAirHumRat(ZoneNum)
      ELSE  ! AirflowSource = outside air
        InletAirHumRat = OutHumRat
      END IF
      ZoneTemp    = MAT(ZoneNum)  ! this should be Tin (account for different reference temps)
      CpAirOutlet = PsyCpAirFnWTdb(InletAirHumRat,TAirflowGapOutletC)
      CpAirZone   = PsyCpAirFnWTdb(ZoneAirHumRat(ZoneNum),ZoneTemp)
      ConvHeatGainToZoneAir = TotAirflowGap * (CpAirOutlet*(TAirflowGapOutletC) - CpAirZone*ZoneTemp)
      IF (SurfaceWindow(SurfNum)%AirFlowDestination == AirFlowWindow_Destination_IndoorAir) THEN
        SurfaceWindow(SurfNum)%ConvHeatGainToZoneAir = ConvHeatGainToZoneAir
        WinHeatGain(SurfNum) = WinHeatGain(SurfNum) + ConvHeatGainToZoneAir
      ELSE
        SurfaceWindow(SurfNum)%RetHeatGainToZoneAir = ConvHeatGainToZoneAir
      END IF
    END IF
    ! For AirflowDestination = ReturnAir in a controlled (i.e., conditioned) zone with return air, see CalcZoneLeavingConditions
    ! for calculation of modification of return-air temperature due to airflow from window gaps into return air.
  END IF

      ! Correct WinHeatGain for interior diffuse shortwave (solar and shortwave from lights) transmitted
      ! back out window
  ConstrNum = Surface(SurfNum)%Construction
  ConstrNumSh = Surface(SurfNum)%ShadedConstruction
  IF(SurfaceWindow(SurfNum)%StormWinFlag==1) THEN
    ConstrNum = Surface(SurfNum)%StormWinConstruction
    ConstrNumSh = Surface(SurfNum)%StormWinShadedConstruction
  END IF

  TransDiff = Construct(ConstrNum)%TransDiff   ! Default value for TransDiff here
  IF(ShadeFlag <= 0) THEN
    TransDiff = Construct(ConstrNum)%TransDiff
  ELSE IF(ShadeFlag==IntShadeOn .OR. ShadeFlag==ExtShadeOn .OR. ShadeFlag==BGShadeOn .OR. ShadeFlag==ExtScreenOn) THEN
    TransDiff = Construct(ConstrNumSh)%TransDiff
  ELSE IF(ShadeFlag==IntBlindOn .OR. ShadeFlag==ExtBlindOn .OR.ShadeFlag==BGBlindOn) THEN
    TransDiff = InterpSlatAng(SurfaceWindow(SurfNum)%SlatAngThisTS,SurfaceWindow(SurfNum)%MovableSlats, &
                               Construct(ConstrNumSh)%BlTransDiff)
  ELSE IF(ShadeFlag == SwitchableGlazing) THEN
    TransDiff = InterpSW(SurfaceWindow(SurfNum)%SwitchingFactor,Construct(ConstrNum)%TransDiff, &
                               Construct(ConstrNumSh)%TransDiff)
  END IF
  WinHeatGain(SurfNum) = WinHeatGain(SurfNum) - QS(Surface(SurfNum)%Zone) * Surface(SurfNum)%Area * TransDiff
    ! shouldn't this be + outward flowing fraction of absorbed SW? -- do not know whose comment this is?  LKL (9/2012)
  WinLossSWZoneToOutWinRep(SurfNum) = QS(Surface(SurfNum)%Zone) * Surface(SurfNum)%Area * TransDiff

  IF(ShadeFlag==IntShadeOn.OR.ShadeFlag==ExtShadeOn.OR.ShadeFlag==IntBlindOn.OR.ShadeFlag==ExtBlindOn.OR. &
     ShadeFlag==BGShadeOn.OR.ShadeFlag==BGBlindOn.OR.ShadeFlag==ExtScreenOn) THEN
    WinShadingAbsorbedSolar(SurfNum) = (SurfaceWindow(SurfNum)%ExtBeamAbsByShade + &
                                       SurfaceWindow(SurfNum)%ExtDiffAbsByShade) * &
                                       (Surface(SurfNum)%Area+SurfaceWindow(SurfNum)%DividerArea)
    WinShadingAbsorbedSolarEnergy(SurfNum) = WinShadingAbsorbedSolar(SurfNum) * TimeStepZone * SecInHour
  END IF
  IF(SunIsUp) THEN
    WinSysSolTransmittance(SurfNum) = WinTransSolar(SurfNum) / &
      (QRadSWOutIncident(SurfNum)*(Surface(SurfNum)%Area+SurfaceWindow(SurfNum)%DividerArea)+0.0001d0)
    WinSysSolAbsorptance(SurfNum)   = (QRadSWwinAbsTot(SurfNum)+WinShadingAbsorbedSolar(SurfNum)) / &
      (QRadSWOutIncident(SurfNum)*(Surface(SurfNum)%Area+SurfaceWindow(SurfNum)%DividerArea)+0.0001d0)
    WinSysSolReflectance(SurfNum)   = 1.0d0 - WinSysSolTransmittance(SurfNum) - WinSysSolAbsorptance(SurfNum)
  ELSE
    WinSysSolTransmittance(SurfNum) = 0.0d0
    WinSysSolAbsorptance(SurfNum)   = 0.0d0
    WinSysSolReflectance(SurfNum)   = 0.0d0
  END IF

     ! Save hcv for use in divider calc with interior or exterior shade (see CalcWinFrameAndDividerTemps)
  IF(ShadeFlag==IntShadeOn .OR. ShadeFlag==ExtShadeOn .OR. ShadeFlag==IntBlindOn .OR. ShadeFlag==ExtBlindOn .OR. &
     ShadeFlag==ExtScreenOn) SurfaceWindow(SurfNum)%ConvCoeffWithShade = hcv
ELSE
     ! No convergence after MaxIterations even with relaxed error tolerance
  CALL ShowSevereError('Convergence error in SolveForWindowTemperatures for window '&
    //TRIM(Surface(SurfNum)%Name))
  Call ShowContinueErrorTimestamp(' ')

  If (DisplayExtraWarnings) Then
    !report out temperatures
    DO i = 1,nglfacep
      CAll ShowContinueError('Glazing face index = '//Trim(RoundSigDigits(i,1))// &
                         ' ; new temperature ='//Trim(RoundSigDigits(thetas(i)-KelvinConv,4))//'C '// &
                         ' ; previous temperature = '//Trim(RoundSigDigits(thetasPrev(i)-KelvinConv, 4))//'C' )
    END DO

  ENDIF

  CALL ShowFatalError('Program halted because of convergence error in SolveForWindowTemperatures for window '&
    //TRIM(Surface(SurfNum)%Name))

END IF

RETURN
END SUBROUTINE SolveForWindowTemperatures