AirflowNetworkSolver.f90 Source File

This File Depends On

sourcefile~~airflownetworksolver.f90~~EfferentGraph sourcefile~airflownetworksolver.f90 AirflowNetworkSolver.f90 sourcefile~dataloopnode.f90 DataLoopNode.f90 sourcefile~dataloopnode.f90->sourcefile~airflownetworksolver.f90 sourcefile~general.f90 General.f90 sourcefile~general.f90->sourcefile~airflownetworksolver.f90 sourcefile~psychroutines.f90 PsychRoutines.f90 sourcefile~general.f90->sourcefile~psychroutines.f90 sourcefile~dataenvironment.f90 DataEnvironment.f90 sourcefile~general.f90->sourcefile~dataenvironment.f90 sourcefile~datahvacglobals.f90 DataHVACGlobals.f90 sourcefile~datahvacglobals.f90->sourcefile~airflownetworksolver.f90 sourcefile~datahvacglobals.f90->sourcefile~general.f90 sourcefile~datasurfaces.f90 DataSurfaces.f90 sourcefile~datasurfaces.f90->sourcefile~airflownetworksolver.f90 sourcefile~datasurfaces.f90->sourcefile~general.f90 sourcefile~psychroutines.f90->sourcefile~airflownetworksolver.f90 sourcefile~dataprecisionglobals.f90 DataPrecisionGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~airflownetworksolver.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataloopnode.f90 sourcefile~dataprecisionglobals.f90->sourcefile~general.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datahvacglobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasurfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~psychroutines.f90 sourcefile~dataglobals.f90 DataGlobals.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataglobals.f90 sourcefile~dataairflownetwork.f90 DataAirflowNetwork.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataairflownetwork.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataenvironment.f90 sourcefile~dataairloop.f90 DataAirLoop.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataairloop.f90 sourcefile~datainterfaces.f90 DataInterfaces.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datainterfaces.f90 sourcefile~dataipshortcuts.f90 DataIPShortCuts.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataipshortcuts.f90 sourcefile~dataruntimelanguage.f90 DataRuntimeLanguage.f90 sourcefile~dataprecisionglobals.f90->sourcefile~dataruntimelanguage.f90 sourcefile~inputprocessor.f90 InputProcessor.f90 sourcefile~dataprecisionglobals.f90->sourcefile~inputprocessor.f90 sourcefile~datasizing.f90 DataSizing.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasizing.f90 sourcefile~datasystemvariables.f90 DataSystemVariables.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datasystemvariables.f90 sourcefile~datavectortypes.f90 DataVectorTypes.f90 sourcefile~dataprecisionglobals.f90->sourcefile~datavectortypes.f90 sourcefile~databsdfwindow.f90 DataBSDFWindow.f90 sourcefile~dataprecisionglobals.f90->sourcefile~databsdfwindow.f90 sourcefile~dataglobals.f90->sourcefile~airflownetworksolver.f90 sourcefile~dataglobals.f90->sourcefile~dataloopnode.f90 sourcefile~dataglobals.f90->sourcefile~general.f90 sourcefile~dataglobals.f90->sourcefile~datahvacglobals.f90 sourcefile~dataglobals.f90->sourcefile~datasurfaces.f90 sourcefile~dataglobals.f90->sourcefile~psychroutines.f90 sourcefile~dataglobals.f90->sourcefile~dataairflownetwork.f90 sourcefile~dataglobals.f90->sourcefile~dataenvironment.f90 sourcefile~dataglobals.f90->sourcefile~dataairloop.f90 sourcefile~dataglobals.f90->sourcefile~dataipshortcuts.f90 sourcefile~dataglobals.f90->sourcefile~dataruntimelanguage.f90 sourcefile~dataglobals.f90->sourcefile~inputprocessor.f90 sourcefile~sortandstringutilities.f90 SortAndStringUtilities.f90 sourcefile~dataglobals.f90->sourcefile~sortandstringutilities.f90 sourcefile~dataoutputs.f90 DataOutputs.f90 sourcefile~dataglobals.f90->sourcefile~dataoutputs.f90 sourcefile~dataglobals.f90->sourcefile~datasizing.f90 sourcefile~dataglobals.f90->sourcefile~databsdfwindow.f90 sourcefile~dataairflownetwork.f90->sourcefile~airflownetworksolver.f90 sourcefile~dataenvironment.f90->sourcefile~airflownetworksolver.f90 sourcefile~dataenvironment.f90->sourcefile~psychroutines.f90 sourcefile~dataairloop.f90->sourcefile~airflownetworksolver.f90 sourcefile~datainterfaces.f90->sourcefile~airflownetworksolver.f90 sourcefile~datainterfaces.f90->sourcefile~general.f90 sourcefile~datainterfaces.f90->sourcefile~psychroutines.f90 sourcefile~datainterfaces.f90->sourcefile~dataenvironment.f90 sourcefile~datainterfaces.f90->sourcefile~dataruntimelanguage.f90 sourcefile~datainterfaces.f90->sourcefile~inputprocessor.f90 sourcefile~dataipshortcuts.f90->sourcefile~general.f90 sourcefile~dataipshortcuts.f90->sourcefile~inputprocessor.f90 sourcefile~dataruntimelanguage.f90->sourcefile~general.f90 sourcefile~datastringglobals.f90 DataStringGlobals.f90 sourcefile~datastringglobals.f90->sourcefile~general.f90 sourcefile~datastringglobals.f90->sourcefile~inputprocessor.f90 sourcefile~datastringglobals.f90->sourcefile~datasystemvariables.f90 sourcefile~inputprocessor.f90->sourcefile~general.f90 sourcefile~sortandstringutilities.f90->sourcefile~inputprocessor.f90 sourcefile~dataoutputs.f90->sourcefile~inputprocessor.f90 sourcefile~datasizing.f90->sourcefile~inputprocessor.f90 sourcefile~datasystemvariables.f90->sourcefile~inputprocessor.f90 sourcefile~datavectortypes.f90->sourcefile~datasurfaces.f90 sourcefile~datavectortypes.f90->sourcefile~databsdfwindow.f90 sourcefile~databsdfwindow.f90->sourcefile~datasurfaces.f90
Help

Files Dependent On This One

sourcefile~~airflownetworksolver.f90~~AfferentGraph sourcefile~airflownetworksolver.f90 AirflowNetworkSolver.f90 sourcefile~airflownetworkbalancemanager.f90 AirflowNetworkBalanceManager.f90 sourcefile~airflownetworksolver.f90->sourcefile~airflownetworkbalancemanager.f90 sourcefile~systemavailabilitymanager.f90 SystemAvailabilityManager.f90 sourcefile~airflownetworkbalancemanager.f90->sourcefile~systemavailabilitymanager.f90 sourcefile~hvacmanager.f90 HVACManager.f90 sourcefile~airflownetworkbalancemanager.f90->sourcefile~hvacmanager.f90 sourcefile~systemavailabilitymanager.f90->sourcefile~hvacmanager.f90 sourcefile~zoneequipmentmanager.f90 Zoneequipmentmanager.f90 sourcefile~systemavailabilitymanager.f90->sourcefile~zoneequipmentmanager.f90 sourcefile~simairservingzones.f90 SimAirServingZones.f90 sourcefile~systemavailabilitymanager.f90->sourcefile~simairservingzones.f90 sourcefile~plantmanager.f90 PlantManager.f90 sourcefile~systemavailabilitymanager.f90->sourcefile~plantmanager.f90 sourcefile~heatbalanceairmanager.f90 HeatBalanceAirManager.f90 sourcefile~systemavailabilitymanager.f90->sourcefile~heatbalanceairmanager.f90 sourcefile~hvacmanager.f90->sourcefile~heatbalanceairmanager.f90 sourcefile~simulationmanager.f90 SimulationManager.f90 sourcefile~hvacmanager.f90->sourcefile~simulationmanager.f90 sourcefile~zoneequipmentmanager.f90->sourcefile~hvacmanager.f90 sourcefile~sizingmanager.f90 SizingManager.f90 sourcefile~zoneequipmentmanager.f90->sourcefile~sizingmanager.f90 sourcefile~simairservingzones.f90->sourcefile~hvacmanager.f90 sourcefile~simairservingzones.f90->sourcefile~sizingmanager.f90 sourcefile~plantmanager.f90->sourcefile~hvacmanager.f90 sourcefile~plantmanager.f90->sourcefile~simulationmanager.f90 sourcefile~utilityroutines.f90 UtilityRoutines.f90 sourcefile~plantmanager.f90->sourcefile~utilityroutines.f90 sourcefile~heatbalanceairmanager.f90->sourcefile~simulationmanager.f90 sourcefile~heatbalancesurfacemanager.f90 HeatBalanceSurfaceManager.f90 sourcefile~heatbalanceairmanager.f90->sourcefile~heatbalancesurfacemanager.f90 sourcefile~sizingmanager.f90->sourcefile~simulationmanager.f90 sourcefile~energyplus.f90 EnergyPlus.f90 sourcefile~simulationmanager.f90->sourcefile~energyplus.f90 sourcefile~simulationmanager.f90->sourcefile~utilityroutines.f90 sourcefile~heatbalancesurfacemanager.f90->sourcefile~simulationmanager.f90 sourcefile~heatbalancemanager.f90 HeatBalanceManager.f90 sourcefile~heatbalancesurfacemanager.f90->sourcefile~heatbalancemanager.f90 sourcefile~heatbalancemanager.f90->sourcefile~sizingmanager.f90 sourcefile~heatbalancemanager.f90->sourcefile~simulationmanager.f90
Help


Source Code

! define this variable to get new code, commenting should yield original
#define SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS

MODULE AirflowNetworkSolver


          ! MODULE INFORMATION:
          !       AUTHOR         Lixing Gu, Don Shirey, and Muthusamy V. Swami
          !       DATE WRITTEN   Jul. 2005
          !       MODIFIED       na
          !       RE-ENGINEERED  na


          ! PURPOSE OF THIS MODULE:
          ! This module is used to simulate airflows and pressures. The module is modified to
          ! meet requirements of EnergyPLus based on AIRNET, developed by
          ! National Institute of Standards and Technology (NIST).

          ! METHODOLOGY EMPLOYED:
          ! An airflow network approach is used. It consists of nodes connected by airflow elements.
          ! The Newton's method is applied to solve a sparse matrix. When a new solver is available, this
          ! module will be replaced or updated.

          ! REFERENCES:
          ! Walton, G. N., 1989, "AIRNET - A Computer Program for Building Airflow Network Modeling,"
          ! NISTIR 89-4072, National Institute of Standards and Technology, Gaithersburg, Maryland

          ! OTHER NOTES: none


          ! USE STATEMENTS:

    USE DataPrecisionGlobals
    USE DataGlobals, ONLY: MaxNameLength, PI, DegToRadians, KelvinConv,rTinyValue
    USE DataInterfaces, ONLY: ShowFatalError, ShowWarningError, ShowSevereError, ShowMessage, ShowContinueError,  &
                           SetupOutputVariable, ShowContinueErrorTimeStamp, ShowRecurringWarningErrorAtEnd
    USE DataEnvironment, ONLY: StdBaroPress, OutBaroPress, OutDryBulbTemp, OutHumRat,  &
                               Latitude
    USE DataSurfaces, ONLY: Surface
    USE Psychrometrics, ONLY: PsyRhoAirFnPbTdbW,PsyCpAirFnWTdb,PsyHFnTdbW
    USE DataAirflowNetwork


IMPLICIT NONE   ! Enforce explicit typing of all variables

INTEGER, PUBLIC  :: NetworkNumOfLinks = 0
INTEGER, PUBLIC  :: NetworkNumOfNodes = 0

PRIVATE ! Everything private unless explicitly made public

INTEGER, PARAMETER :: NrInt=20 ! Number of intervals for a large opening

! Common block AFEDAT
REAL(r64), ALLOCATABLE, DIMENSION(:) :: AFECTL
REAL(r64), ALLOCATABLE, DIMENSION(:) :: AFLOW2
REAL(r64), ALLOCATABLE, DIMENSION(:) :: AFLOW
REAL(r64), ALLOCATABLE, DIMENSION(:) :: PS
REAL(r64), ALLOCATABLE, DIMENSION(:) :: PW


! Common block CONTRL
REAL(r64), PRIVATE :: PB=0.0d0
INTEGER, PRIVATE :: LIST=0


! Common block ZONL
REAL(r64), ALLOCATABLE, DIMENSION(:) :: RHOZ
REAL(r64), ALLOCATABLE, DIMENSION(:) :: SQRTDZ
REAL(r64), ALLOCATABLE, DIMENSION(:) :: VISCZ
REAL(r64), ALLOCATABLE, DIMENSION(:) :: SUMAF
REAL(r64), ALLOCATABLE, DIMENSION(:) :: TZ              ! Temperature [C]
REAL(r64), ALLOCATABLE, DIMENSION(:) :: WZ              ! Humidity ratio [kg/kg]
REAL(r64), ALLOCATABLE, DIMENSION(:) :: PZ  ! Pressure [Pa]

! Other array variables
INTEGER, ALLOCATABLE, DIMENSION(:) :: ID
INTEGER, ALLOCATABLE, DIMENSION(:) :: IK
REAL(r64), ALLOCATABLE, DIMENSION(:) :: AD
REAL(r64), ALLOCATABLE, DIMENSION(:) :: AU

#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS
INTEGER, ALLOCATABLE, DIMENSION(:) :: newIK   !noel
REAL(r64), ALLOCATABLE, DIMENSION(:) :: newAU !noel
#endif

!REAL(r64), ALLOCATABLE, DIMENSION(:) :: AL
REAL(r64), ALLOCATABLE, DIMENSION(:) :: SUMF
INTEGER :: Unit11=0
INTEGER :: Unit21=0

! Large opening variables
REAL(r64), ALLOCATABLE, DIMENSION(:) :: DpProf    ! Differential pressure profile for Large Openings [Pa]
REAL(r64), ALLOCATABLE, DIMENSION(:) :: RhoProfF  ! Density profile in FROM zone [kg/m3]
REAL(r64), ALLOCATABLE, DIMENSION(:) :: RhoProfT  ! Density profile in TO zone [kg/m3]
REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: DpL     ! Array of stack pressures in link

PUBLIC SETSKY
PUBLIC AIRMOV
Public AllocateAirflowNetworkData
Public InitAirflowNetworkData


CONTAINS


SUBROUTINE AllocateAirflowNetworkData


          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Aug. 2003
          !       MODIFIED       na
          !       RE-ENGINEERED  na


          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine allocates dynamic arrays for AirflowNetworkSolver.


          ! METHODOLOGY EMPLOYED:
          ! na


          ! REFERENCES:
          ! na


          ! USE STATEMENTS:
          ! na


  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 i,j,N
INTEGER, EXTERNAL :: GetNewUnitNumber  ! External  function to "get" a unit number

! Assume a network to simulate multizone airflow is a subset of the network to simulate air distribution system.
! Network array size is allocated based on the network of air distribution system.
! If multizone airflow is simulated only, the array size is allocated based on the multizone network.
          ! FLOW:
   NetworkNumOfLinks = AirflowNetworkNumOfLinks
   NetworkNumOfNodes = AirflowNetworkNumOfNodes

   ALLOCATE(AFECTL(NetworkNumOfLinks))
   ALLOCATE(AFLOW2(NetworkNumOfLinks))
   ALLOCATE(AFLOW(NetworkNumOfLinks))
   ALLOCATE(PW(NetworkNumOfLinks))
   ALLOCATE(PS(NetworkNumOfLinks))

   ALLOCATE(TZ(NetworkNumOfNodes))
   ALLOCATE(WZ(NetworkNumOfNodes))
   ALLOCATE(PZ(NetworkNumOfNodes))
   ALLOCATE(RHOZ(NetworkNumOfNodes))
   ALLOCATE(SQRTDZ(NetworkNumOfNodes))
   ALLOCATE(VISCZ(NetworkNumOfNodes))
   ALLOCATE(SUMAF(NetworkNumOfNodes))


   ALLOCATE(ID(NetworkNumOfNodes))
   ALLOCATE(IK(NetworkNumOfNodes+1))
#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS
   ALLOCATE(newIK(NetworkNumOfNodes+1))
#endif
   ALLOCATE(AD(NetworkNumOfNodes))
   ALLOCATE(SUMF(NetworkNumOfNodes))

   N = 0
   Do i=1,AirflowNetworkNumOfLinks
      j = AirflowNetworkCompData(AirflowNetworkLinkageData(i)%CompNum)%CompTypeNum
      if (j .EQ. CompTypeNum_DOP) then
         n = n+1
      end if
   End do

   ALLOCATE(DpProf(N*(NrInt+2)))
   ALLOCATE(RhoProfF(N*(NrInt+2)))
   ALLOCATE(RhoProfT(N*(NrInt+2)))
   ALLOCATE(DpL(2,AirflowNetworkNumOfLinks))

   PB = 101325.0d0
!   LIST = 5
   LIST = 0
   If (LIST .GE. 1) then
      Unit21 = GetNewUnitNumber()
      open(Unit21,file='eplusADS.out')
   end if


   DO N=1,NetworkNumOfNodes
      ID(N) = N
   END DO
   DO I=1,NetworkNumOfLinks
      AFECTL(I) = 1.0d0
      AFLOW(I) = 0.0d0
      AFLOW2(I) = 0.0d0
   END DO


   DO I=1,NetworkNumOfNodes
      TZ(I) = AirflowNetworkNodeSimu(I)%TZ
      WZ(I) = AirflowNetworkNodeSimu(I)%WZ
      PZ(I) = AirflowNetworkNodeSimu(I)%PZ
   END DO

   ! Assign linkage values
   DO I=1,NetworkNumOfLinks
      PW(I) = 0.0d0
   END DO
   ! Write an ouput file used for AIRNET input
   if (LIST .GE. 5) then
     Unit11 = GetNewUnitNumber()
     open(Unit11,file='eplusADS.inp')
     DO i=1,NetworkNumOfNodes
        write(Unit11,901) i,AirflowNetworkNodeData(I)%NodeTypeNum,AirflowNetworkNodeData(I)%NodeHeight,TZ(I),PZ(I)
     end do
     write(Unit11,900) 0
     DO i=1,AirflowNetworkNumOfComps
        j = AirflowNetworkCompData(i)%TypeNum
        select case (AirflowNetworkCompData(i)%CompTypeNum)
           CASE (CompTypeNum_PLR) !'PLR'  Power law component
!              WRITE(Unit11,902) AirflowNetworkCompData(i)%CompNum,1,DisSysCompLeakData(j)%FlowCoef, &
!                  DisSysCompLeakData(j)%FlowCoef,DisSysCompLeakData(j)%FlowCoef,DisSysCompLeakData(j)%FlowExpo
           CASE (CompTypeNum_SCR) !'SCR'  Surface crack component
              WRITE(Unit11,902) AirflowNetworkCompData(i)%CompNum,1,MultizoneSurfaceCrackData(j)%FlowCoef, &
                  MultizoneSurfaceCrackData(j)%FlowCoef,MultizoneSurfaceCrackData(j)%FlowCoef,MultizoneSurfaceCrackData(j)%FlowExpo
           CASE (CompTypeNum_DWC) !'DWC' Duct component
!              WRITE(Unit11,902) AirflowNetworkCompData(i)%CompNum,2,DisSysCompDuctData(j)%L,DisSysCompDuctData(j)%D, &
!                               DisSysCompDuctData(j)%A,DisSysCompDuctData(j)%Rough
!              WRITE(Unit11,903) DisSysCompDuctData(i)%TurDynCoef,DisSysCompDuctData(j)%LamFriCoef, &
!                               DisSysCompDuctData(j)%LamFriCoef,DisSysCompDuctData(j)%InitLamCoef
!           CASE (CompTypeNum_CVF) ! 'CVF' Constant volume fan component
!              WRITE(Unit11,904) AirflowNetworkCompData(i)%CompNum,4,DisSysCompCVFData(j)%FlowRate
           CASE (CompTypeNum_EXF) ! 'EXF' Zone exhaust fan
              WRITE(Unit11,904) AirflowNetworkCompData(i)%CompNum,4,MultizoneCompExhaustFanData(j)%FlowRate
           CASE Default
        END SELECT
     end do
     write(Unit11,900) 0
     DO i=1,NetworkNumOfLinks
        write(Unit11,910) i,AirflowNetworkLinkageData(I)%NodeNums(1),AirflowNetworkLinkageData(I)%NodeHeights(1), &
                         AirflowNetworkLinkageData(I)%NodeNums(2),AirflowNetworkLinkageData(I)%NodeHeights(2), &
                         AirflowNetworkLinkageData(I)%CompNum,0,0
     end do
     write(Unit11,900) 0
   end if

   CALL SETSKY

   !SETSKY figures out the IK stuff -- which is why E+ doesn't allocate AU until here
#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS
!   ! only printing to screen, can be commented
!   print*, "SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS is defined"
!
!   write(*,'(2(a,i8))') "AllocateAirflowNetworkData: after SETSKY, allocating AU.  NetworkNumOfNodes=", &
!        NetworkNumOfNodes, " IK(NetworkNumOfNodes+1)= NNZE=", IK(NetworkNumOfNodes+1)
!   print*, " NetworkNumOfLinks=", NetworkNumOfLinks
!
   ! allocate same size as others -- this will be maximum  !noel
   ALLOCATE(newAU(IK(NetworkNumOfNodes+1)))
#endif

   ! noel, GNU says the AU is indexed above its upper bound
   !ALLOCATE(AU(IK(NetworkNumOfNodes+1)-1))
   ALLOCATE(AU(IK(NetworkNumOfNodes+1)))

900 Format(1X,i2)
901 Format(1X,2I4,4F9.4)
902 Format(1X,2I4,4F9.4)
903 Format(9X,4F9.4)
904 Format(1X,2I4,1F9.4)
910 Format(1X,I4,2(I4,F9.4),I4,2F4.1)

RETURN
END SUBROUTINE AllocateAirflowNetworkData

SUBROUTINE InitAirflowNetworkData

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Aug. 2003
          !       MODIFIED       na
          !       RE-ENGINEERED  na


          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine initializes variables for AirflowNetworkSolver.


          ! METHODOLOGY EMPLOYED:
          ! na


          ! REFERENCES:
          ! na


          ! USE STATEMENTS:
          ! na

  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 I,N

      ! FLOW:
   DO N=1,NetworkNumOfNodes
      ID(N) = N
   END DO
   DO I=1,NetworkNumOfLinks
      AFECTL(I) = 1.0d0
      AFLOW(I) = 0.0d0
      AFLOW2(I) = 0.0d0
   END DO

   DO I=1,NetworkNumOfNodes
      TZ(I) = AirflowNetworkNodeSimu(I)%TZ
      WZ(I) = AirflowNetworkNodeSimu(I)%WZ
      PZ(I) = AirflowNetworkNodeSimu(I)%PZ
   END DO

RETURN

END SUBROUTINE InitAirflowNetworkData


      SUBROUTINE SETSKY
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   1998
          !       MODIFIED       Feb. 2006 (L. Gu) to meet requirements of AirflowNetwork
          !       RE-ENGINEERED  na


          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine sets up the "IK" array describing the sparse matrix [A] in skyline
          !     form by using the location matrix.

          ! METHODOLOGY EMPLOYED:
          ! na


          ! REFERENCES:
          ! AIRNET


          ! USE STATEMENTS:
          ! na

  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:
          ! IK(K) - pointer to the top of column/row "K".
      INTEGER  I, J, K, L, M, N1, N2

      ! FLOW:
      ! Initialize "IK".
      DO I=1,NetworkNumOfNodes+1
         IK(I) = 0
      END DO
      ! Determine column heights.
      DO M=1,NetworkNumOfLinks
        J = AirflowNetworkLinkageData(M)%NodeNums(2)
        IF(J.EQ.0) CYCLE
        L = ID(J)
        I = AirflowNetworkLinkageData(M)%NodeNums(1)
        K = ID(I)
        N1 = ABS(L-K)
        N2 = MAX(K,L)
        IK(N2) = MAX(IK(N2),N1)
      END DO
      ! Convert heights to column addresses.
      J = IK(1)
      IK(1) = 1
      DO K=1,NetworkNumOfNodes
        I = IK(K+1)
        IK(K+1) = IK(K)+J
        J = I
      END DO

      RETURN
      END SUBROUTINE SETSKY


      SUBROUTINE AIRMOV
!
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine is a driver for AIRNET to calculate nodal pressures and linkage airflows

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:

  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  I, M, N, ITER

      ! FLOW:
      ! Compute zone air properties.
      DO N=1,NetworkNumOfNodes
        RHOZ(N) = PsyRhoAirFnPbTdbW(StdBaroPress+PZ(N),TZ(N),WZ(N))
        if (AirflowNetworkNodeData(N)%ExtNodeNum > 0) then
           RHOZ(N) = PsyRhoAirFnPbTdbW(StdBaroPress+PZ(N),OutDryBulbTemp,OutHumRat)
           TZ(N) = OutDryBulbTemp
           WZ(N) = OutHumRat
        end if
        SQRTDZ(N) = SQRT(RHOZ(N))
        VISCZ(N) = 1.71432d-5+4.828d-8*TZ(N)
        IF(LIST.GE.2) WRITE(Unit21,903) 'D,V:',N,RHOZ(N),VISCZ(N)
      END DO
      IF(AirflowNetworkSimu%InitFlag.EQ.0) THEN
        DO N=1,NetworkNumOfNodes
          IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.0) PZ(N) = 0.0d0
        END DO
      END IF
      ! Compute stack pressures.
      DO I=1,NetworkNumOfLinks
        N = AirflowNetworkLinkageData(I)%NodeNums(1)
        M = AirflowNetworkLinkageData(I)%NodeNums(2)
        IF(AFLOW(I).GT.0.0d0) THEN
          PS(I) = 9.80d0*( RHOZ(N)*(AirflowNetworkNodeData(N)%NodeHeight-AirflowNetworkNodeData(M)%NodeHeight) &
                    +AirflowNetworkLinkageData(I)%NodeHeights(2)*(RHOZ(M)-RHOZ(N)) )
        ELSE IF(AFLOW(I).LT.0.0d0) THEN
          PS(I) = 9.80d0*( RHOZ(M)*(AirflowNetworkNodeData(N)%NodeHeight-AirflowNetworkNodeData(M)%NodeHeight) &
                     +AirflowNetworkLinkageData(I)%NodeHeights(1)*(RHOZ(M)-RHOZ(N)) )
        ELSE
          PS(I) = 4.90d0*( (RHOZ(N)+RHOZ(M))*(AirflowNetworkNodeData(N)%NodeHeight-AirflowNetworkNodeData(M)%NodeHeight) &
            + (AirflowNetworkLinkageData(I)%NodeHeights(1)+AirflowNetworkLinkageData(I)%NodeHeights(2))*(RHOZ(M)-RHOZ(N)))
        END IF
      END DO

      ! Calculate pressure field in a large opening
      CALL PSTACK

      CALL SOLVZP(IK,AD,AU,ITER)

      ! Report element flows and zone pressures.
      DO N=1,NetworkNumOfNodes
        SUMAF(N) = 0.0d0
      END DO
      IF(LIST.GE.1) WRITE(Unit21,900)
      DO I=1,NetworkNumOfLinks
        N = AirflowNetworkLinkageData(I)%NodeNums(1)
        M = AirflowNetworkLinkageData(I)%NodeNums(2)
        IF(LIST.GE.1) THEN
          WRITE(Unit21,901) 'Flow: ',I,N,M,AirflowNetworkLinkSimu(I)%DP,AFLOW(I),AFLOW2(I)
        END IF
        SUMAF(N) = SUMAF(N)-AFLOW(I)-AFLOW2(I)
        SUMAF(M) = SUMAF(M)+AFLOW(I)+AFLOW2(I)
      END DO
      DO N=1,NetworkNumOfNodes
        IF(LIST.GE.1) WRITE(Unit21,903) 'Room: ',N,PZ(N),SUMAF(N),TZ(N)
      END DO

   DO I=1,NetworkNumOfLinks
      if (AFLOW2(i) .ne. 0.0d0) then
         Continue
      end if
      if (AFLOW(i) > 0.0d0) then
         AirflowNetworkLinkSimu(I)%FLOW = AFLOW(I)
         AirflowNetworkLinkSimu(I)%FLOW2 = 0.0d0
      else
         AirflowNetworkLinkSimu(I)%FLOW = 0.0d0
         AirflowNetworkLinkSimu(I)%FLOW2 = -AFLOW(I)
      end if
      If (AirflowNetworkCompData(AirflowNetworkLinkageData(I)%CompNum)%CompTypeNum == CompTypeNum_HOP) Then
         if (AFLOW2(I) .ne. 0.0d0) then
            AirflowNetworkLinkSimu(I)%FLOW = AFLOW(I)+AFLOW2(I)
            AirflowNetworkLinkSimu(I)%FLOW2 = AFLOW2(I)
         end if
      end if
      if (AirflowNetworkLinkageData(I)%DetOpenNum > 0) then
         if (AFLOW2(I) .ne. 0.0d0) then
            AirflowNetworkLinkSimu(I)%FLOW = AFLOW(I)+AFLOW2(I)
            AirflowNetworkLinkSimu(I)%FLOW2 = AFLOW2(I)
         end if
      end if
      If (AirflowNetworkCompData(AirflowNetworkLinkageData(I)%CompNum)%CompTypeNum == CompTypeNum_SOP .AND. &
          AFLOW2(I) .ne. 0.0d0) then
        if (AFLOW(i) >= 0.0d0) then
         AirflowNetworkLinkSimu(I)%FLOW = AFLOW(I)
         AirflowNetworkLinkSimu(I)%FLOW2 = abs(AFLOW2(I))
        else
         AirflowNetworkLinkSimu(I)%FLOW = abs(AFLOW2(I))
         AirflowNetworkLinkSimu(I)%FLOW2 = -AFLOW(I)
        End If
      end if
   END DO

   DO I=1,NetworkNumOfNodes
      AirflowNetworkNodeSimu(I)%PZ = PZ(I)
   END DO

  900 FORMAT(/,11X,'i    n    m       DP',12x,'F1',12X,'F2')
  901 FORMAT(1X,A6,3I5,3F14.6)
  902 FORMAT(/,11X,'n       P',12x,'sumF')
  903 FORMAT(1X,A6,I5,3F14.6)
  907 FORMAT(/,' CPU seconds for ',A,F12.3)

      RETURN
      END SUBROUTINE AIRMOV


      SUBROUTINE SOLVZP(IK,AD,AU,ITER)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves zone pressures by modified Newton-Raphson iteration

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          USE General, ONLY: RoundSigDigits

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(INOUT)  :: IK(NetworkNumOfNodes+1) ! pointer to the top of column/row "K"
          INTEGER, INTENT(INOUT)  :: ITER  ! number of iterations
          !noel GNU says AU is being indexed beyound bounds
          !REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)-1) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AD(NetworkNumOfNodes)  ! the main diagonal of [A] before and after factoring


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

          !     NNZE   - number of nonzero entries in the "AU" array.
          !     LFLAG   - if = 1, use laminar relationship (initialization).
          !     I       - element number.
          !     N       - number of node/zone 1.
          !     M       - number of node/zone 2.
          !     F       - flows through the element (kg/s).
          !     DF      - partial derivatives:  DF/DP.
          !     NF      - number of flows, 1 or 2.
          !     SUMF    - sum of flows into node/zone.
          !     CCF     - current pressure correction (Pa).
          !     PCF     - previous pressure correction (Pa).
          !     CEF     - convergence enhancement factor.
!
          INTEGER  N, NNZE, NSYM, LFLAG, CONVG, ACCEL
          REAL(r64) PCF(NetworkNumOfNodes), CEF(NetworkNumOfNodes)
          REAL(r64) C, SSUMF, SSUMAF, ACC0, ACC1
          REAL(r64) CCF(NetworkNumOfNodes)

      ! FLOW:
      ACC1 = 0.0d0
      ACCEL = 0
      NSYM = 0
      NNZE = IK(NetworkNumOfNodes+1)-1
      IF(LIST.GE.2) WRITE(Unit21,*) 'Initialization',NetworkNumOfNodes,NetworkNumOfLinks,NNZE
      ITER = 0

      DO N=1,NetworkNumOfNodes
        PCF(N) = 0.0d0
        CEF(N) = 0.0d0
      END DO

      IF(AirflowNetworkSimu%InitFlag.NE.1) THEN
        ! Initialize node/zone pressure values by assuming only linear relationship between
        ! airflows and pressure drops.
        LFLAG = 1
        CALL FILJAC(NNZE,LFLAG)
        DO N=1,NetworkNumOfNodes
          IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.0) PZ(N) = SUMF(N)
        END DO
        ! Data dump.
        IF(LIST.GE.3) THEN
          CALL DUMPVD('AD:',AD,NetworkNumOfNodes,Unit21)
          CALL DUMPVD('AU:',AU,NNZE,Unit21)
          CALL DUMPVR('AF:',SUMF,NetworkNumOfNodes,Unit21)
        END IF
        ! Solve linear system for approximate PZ.
#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS
        CALL FACSKY(newAU,AD,newAU,newIK,NetworkNumOfNodes,NSYM) !noel
        CALL SLVSKY(newAU,AD,newAU,PZ,newIK,NetworkNumOfNodes,NSYM) !noel
#else
        CALL FACSKY(AU,AD,AU,IK,NetworkNumOfNodes,NSYM)
        CALL SLVSKY(AU,AD,AU,PZ,IK,NetworkNumOfNodes,NSYM)
#endif
        IF(LIST.GE.2) CALL DUMPVD('PZ:',PZ,NetworkNumOfNodes,Unit21)
      END IF
      ! Solve nonlinear airflow network equations by modified Newton's method.

      DO WHILE (ITER .LT. AirflowNetworkSimu%MaxIteration)
        LFLAG = 0
        ITER = ITER+1
        IF(LIST.GE.2) WRITE(Unit21,*) 'Begin iteration ',ITER
        ! Set up the Jacobian matrix.
        CALL FILJAC(NNZE,LFLAG)
        ! Data dump.
        IF(LIST.GE.3) THEN
          CALL DUMPVR('SUMF:',SUMF,NetworkNumOfNodes,Unit21)
          CALL DUMPVR('SUMAF:',SUMAF,NetworkNumOfNodes,Unit21)
        END IF
        ! Check convergence.
        CONVG = 1
        SSUMF = 0.0d0
        SSUMAF = 0.0d0
        DO N=1,NetworkNumOfNodes
          SSUMF = SSUMF + ABS(SUMF(N))
          SSUMAF = SSUMAF + SUMAF(N)
          IF(CONVG.EQ.1) THEN
            IF(ABS(SUMF(N)).LE.AirflowNetworkSimu%AbsTol) CYCLE
            IF(ABS(SUMF(N)/SUMAF(N)).GT.AirflowNetworkSimu%RelTol) CONVG = 0
          END IF
        END DO
        ACC0 = ACC1
        IF(SSUMAF.GT.0.0d0) ACC1 = SSUMF / SSUMAF
        IF(CONVG.EQ.1 .AND. ITER.GT.1) RETURN
        IF(ITER.GE.AirflowNetworkSimu%MaxIteration) EXIT
        ! Data dump.
        IF(LIST.GE.3) THEN
          CALL DUMPVD('AD:',AD,NetworkNumOfNodes,Unit21)
          CALL DUMPVD('AU:',AU,NNZE,Unit21)
        END IF
        ! Solve AA * CCF = SUMF.
        DO N=1,NetworkNumOfNodes
          CCF(N) = SUMF(N)
        END DO
#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS
        CALL FACSKY(newAU,AD,newAU,newIK,NetworkNumOfNodes,NSYM) !noel
        CALL SLVSKY(newAU,AD,newAU,CCF,newIK,NetworkNumOfNodes,NSYM) !noel
#else
        CALL FACSKY(AU,AD,AU,IK,NetworkNumOfNodes,NSYM)
        CALL SLVSKY(AU,AD,AU,CCF,IK,NetworkNumOfNodes,NSYM)
#endif
        ! Revise PZ (Steffensen iteration on the N-R correction factors to handle oscillating corrections).
        IF(ACCEL.EQ.1) THEN
          ACCEL = 0
        ELSE
          IF(ITER.GT.2 .AND. ACC1.GT.0.5d0*ACC0) ACCEL = 1
        END IF
        DO N=1,NetworkNumOfNodes
          IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.1) CYCLE
          CEF(N) = 1.0d0
          IF(ACCEL.EQ.1) THEN
            C = CCF(N)/PCF(N)
            IF(C.LT.AirflowNetworkSimu%ConvLimit) CEF(N) = 1.0/(1.0-C)
            C = CCF(N) * CEF(N)
          ELSE
!            IF (CCF(N) .EQ. 0.0) CCF(N)=TINY(CCF(N))  ! 1.0E-40
            IF (CCF(N) .EQ. 0.0d0) CCF(N)=rTinyValue     ! 1.0E-40 (Epsilon)
            PCF(N) = CCF(N)
            C = CCF(N)
          END IF
          IF(ABS(C).GT.AirflowNetworkSimu%MaxPressure) THEN
            CEF(N) = CEF(N)*AirflowNetworkSimu%MaxPressure/ABS(C)
            PZ(N) = PZ(N)-CCF(N)*CEF(N)
          ELSE
            PZ(N) = PZ(N)-C
          END IF
        END DO
        ! Data revision dump.
        IF(LIST.GE.2) THEN
          DO N=1,NetworkNumOfNodes
            IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.0) &
            WRITE(Unit21,901) ' Rev:',N,SUMF(N),CCF(N),CEF(N),PZ(N)
          END DO
        ENDIF
      END DO

      ! Error termination.
      CALL ShowSevereError('Too many iterations (SOLVZP) in Airflow Network simulation')
      AirflowNetworkSimu%ExtLargeOpeningErrCount = AirflowNetworkSimu%ExtLargeOpeningErrCount + 1
      if (AirflowNetworkSimu%ExtLargeOpeningErrCount < 2) then
        CALL ShowWarningError('AirflowNetwork: SOLVER, Changing values for initialization flag, Relative airflow convergence, ' &
         //'Absolute airflow convergence, Convergence acceleration limit or Maximum Iteration Number may solve the problem.')
        CALL ShowContinueErrorTimeStamp(' ')
        CALL ShowContinueError('..Iterations='//TRIM(RoundSigDigits(ITER))//', Max allowed='//  &
               TRIM(RoundSigDigits(AirflowNetworkSimu%MaxIteration)))
        CALL ShowFatalError('AirflowNetwork: SOLVER, The previous error causes termination.')
      else
        CALL ShowRecurringWarningErrorAtEnd('AirFlowNetwork: Too many iterations (SOLVZP) in AirflowNetwork ' &
                                            //' simulation continues.',AirflowNetworkSimu%ExtLargeOpeningErrIndex)
      end if

  901 FORMAT(A5,I3,2E14.6,0PF8.4,F24.14)
      RETURN
      END SUBROUTINE SOLVZP


      SUBROUTINE FILJAC(NNZE,LFLAG)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine creates matrices for solution of flows

          ! 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)  :: NNZE   ! number of nonzero entries in the "AU" array.
          INTEGER, INTENT(IN)  :: LFLAG  ! if = 1, use laminar relationship (initialization).


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

          ! I       - component number.
          ! N       - number of node/zone 1.
          ! M       - number of node/zone 2.
          ! F       - flows through the element (kg/s).
          ! DF      - partial derivatives:  DF/DP.
          ! NF      - number of flows, 1 or 2.
!
       INTEGER  I, J, M, N, FLAG, NF
#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS
       INTEGER LHK,JHK,JHK1, newsum,newh,Nzeros, ispan, thisIK !noel
       LOGICAL allZero ! noel
       LOGICAL, save :: firsttime=.true. ! noel
#endif
       REAL(r64)     X(4)
       REAL(r64) DP, F(2), DF(2)

      ! FLOW:
      DO N=1,NetworkNumOfNodes
        SUMF(N) = 0.0d0
        SUMAF(N) = 0.0d0
        IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.1) THEN
          AD(N) = 1.0d0
        ELSE
          AD(N) = 0.0d0
        END IF
      END DO
      DO N=1,NNZE
        AU(N) = 0.0d0
      END DO
!                              Set up the Jacobian matrix.
      DO I=1,NetworkNumOfLinks
        N = AirflowNetworkLinkageData(I)%NodeNums(1)
        M = AirflowNetworkLinkageData(I)%NodeNums(2)
        !!!! Check array of DP. DpL is used for multizone air flow calculation only
        !!!! and is not for forced air calculation
        if (I .GT. NumOfLinksMultiZone) then
          DP = PZ(N)-PZ(M)+PS(I)+PW(I)
        else
          DP = PZ(N)-PZ(M)+DpL(1,I)+PW(I)
        end if
        IF(LIST.GE.4) WRITE(Unit21,901) 'PS:',I,N,M,PS(I),PW(I),AirflowNetworkLinkSimu(I)%DP
        J = AirflowNetworkLinkageData(I)%CompNum
        Select Case (AirflowNetworkCompData(J)%CompTypeNum)
           Case (CompTypeNum_PLR) ! Distribution system crack component
              CALL AFEPLR(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_DWC) ! Distribution system duct component
              CALL AFEDWC(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_CVF) ! Distribution system constant volume fan component
              CALL AFECFR(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_FAN) ! Distribution system detailed fan component
              CALL AFEFAN(J,LFLAG,DP,I,N,M,F,DF,NF)
!           Case (CompTypeNum_CPF) ! not currently used in EnergyPlus code -- left for compatibility with AirNet
!              CALL AFECPF(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_DMP) ! Distribution system damper component
              CALL AFEDMP(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_ELR) ! Distribution system effective leakage ratio component
              CALL AFEELR(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_CPD) ! Distribution system constant pressure drop component
              CALL AFECPD(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_DOP) ! Detailed opening
              CALL AFEDOP(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_SOP) ! Simple opening
              CALL AFESOP(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_SCR) ! Surface crack component
              CALL AFESCR(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_SEL) ! Surface effective leakage ratio component
              CALL AFESEL(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_COI) ! Distribution system coil component
              CALL AFECOI(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_TMU) ! Distribution system terminal unit component
              CALL AFETMU(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_EXF) ! Exhaust fan component
              CALL AFEEXF(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_HEX) ! Distribution system heat exchanger component
              CALL AFEHEX(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case (CompTypeNum_HOP) ! Horizontal opening
              CALL AFEHOP(J,LFLAG,DP,I,N,M,F,DF,NF)
           Case Default
              CYCLE
        End Select
        AirflowNetworkLinkSimu(I)%DP = DP
        AFLOW(I) = F(1)
        AFLOW2(I) = 0.0d0
        if (AirflowNetworkCompData(J)%CompTypeNum .EQ. CompTypeNum_DOP) then
           AFLOW2(I) = F(2)
        end if
        if (AirflowNetworkCompData(J)%CompTypeNum .EQ. CompTypeNum_HOP) then
           AFLOW2(I) = F(2)
        end if
        IF(LIST.GE.3) WRITE(Unit21,901) ' NRi:',I,N,M,AirflowNetworkLinkSimu(I)%DP,F(1),DF(1)
        FLAG = 1
        IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.0) THEN
          FLAG = FLAG+1
          X(1) =  DF(1)
          X(2) = -DF(1)
          SUMF(N) = SUMF(N)+F(1)
          SUMAF(N) = SUMAF(N)+ABS(F(1))
        END IF
        IF(AirflowNetworkNodeData(M)%NodeTypeNum.EQ.0) THEN
          FLAG = FLAG+2
          X(4) =  DF(1)
          X(3) = -DF(1)
          SUMF(M) = SUMF(M)-F(1)
          SUMAF(M) = SUMAF(M)+ABS(F(1))
        END IF
        IF(FLAG.NE.1) CALL FILSKY(X,AirflowNetworkLinkageData(I)%NodeNums,IK,AU,AD,FLAG)
        IF(NF.EQ.1) CYCLE
        AFLOW2(I) = F(2)
        IF(LIST.GE.3) WRITE(Unit21,901) ' NRj:',I,N,M,AirflowNetworkLinkSimu(I)%DP,F(2),DF(2)
        FLAG = 1
        IF(AirflowNetworkNodeData(N)%NodeTypeNum.EQ.0) THEN
          FLAG = FLAG+1
          X(1) =  DF(2)
          X(2) = -DF(2)
          SUMF(N) = SUMF(N)+F(2)
          SUMAF(N) = SUMAF(N)+ABS(F(2))
        END IF
        IF(AirflowNetworkNodeData(M)%NodeTypeNum.EQ.0) THEN
          FLAG = FLAG+2
          X(4) =  DF(2)
          X(3) = -DF(2)
          SUMF(M) = SUMF(M)-F(2)
          SUMAF(M) = SUMAF(M)+ABS(F(2))
        END IF
        IF(FLAG.NE.1) CALL FILSKY(X,AirflowNetworkLinkageData(I)%NodeNums,IK,AU,AD,FLAG)
      END DO
  901 FORMAT(A5,3I3,4E16.7)

#ifdef SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS

     ! After the matrix values have been set, we can look at them and see if any columns are filled with zeros.
     ! If they are, let's remove them from the matrix -- but only for the purposes of doing the solve.
     ! They way I do this is building a separate IK array (newIK) that simply changes the column heights.
     ! So the affected SOLVEs would use this newIK and nothing else changes.
     DO n=1,NetworkNumOfNodes+1
        newIK(n)=IK(n)
        !print*, " NetworkNumOfNodes  n=", n, " IK(n)=", IK(n)
     enddo

     newsum=IK(2)-IK(1) ! always 0?

     JHK=1
     DO n=2,NetworkNumOfNodes
        JHK1 = IK(n+1)  ! starts at IK(3)-IK(2)
        LHK = JHK1-JHK
        IF(LHK.LE.0) then
          newIK(n+1) = newIK(n)
          CYCLE
        endif
        !write(*,'(4(a,i8))') "n=", n, " ik=", ik(n), " JHK=", JHK, " LHK=", LHK

        ! is the entire column zero?  noel
        allZero=.True.
        DO i=0,LHK-1
           if (AU(JHK+i) .ne. 0.0d0) then
              allZero=.False.
              exit
           endif
        enddo

        newh=LHK
        if (allZero .eqv. .True.) then
           !print*, "allzero n=", n
           newh=0
        else
           !DO i=0,LHK-1
           !   write(*, '(2(a,i8),a, f15.3)') "  n=", n, " i=", i, " AU(JHK+i)=", AU(JHK+i)
           !enddo
        endif
        newIK(n+1) = newIK(n) + newh
        newsum = newsum + newh

        !do i = LHK-1,0, -1
        !   write(*, '(2(a,i8),a, f15.3)') "  n=", n, " i=", i, " AU(JHK+i)=", AU(JHK+i)
        !enddo
        JHK = JHK1
     enddo

     ! this is just a print to screen, is not necessary
!     if (firsttime) then
!        write(*, '(2(a,i8))') " After SKYLINE_MATRIX_REMOVE_ZERO_COLUMNS: newsum=", newsum, " oldsum=", IK(NetworkNumOfNodes+1)
!        firsttime=.false.
!     endif


     ! Now fill newAU from AU, using newIK
     thisIK=1
     DO n=2,NetworkNumOfNodes
        thisIK = newIK(n)
        ispan = newIK(n+1) - thisIK

        IF(ispan.LE.0) CYCLE
        DO i=0,ispan-1
           newAU(thisIK+i) = AU(IK(n)+i)
        enddo

     enddo
#endif

      RETURN
      END SUBROUTINE FILJAC


      SUBROUTINE AFEPLR(J,LFLAG,PDROP,I,N,M,F,DF,NF)
!
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a power law resistance airflow component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     CDM, FL, FT
         REAL(r64)     RhozNorm, VisczNorm, expn, Ctl, coef
         INTEGER  CompNum

      ! FLOW:
      ! Crack standard condition: T=20C, p=101325 Pa and 0 g/kg
      CompNum = AirflowNetworkCompData(J)%TypeNum
      RhozNorm = PsyRhoAirFnPbTdbW(101325.0d0,20.0d0,0.0d0)
      VisczNorm = 1.71432d-5+4.828d-8*20.0d0
      expn = DisSysCompLeakData(CompNum)%FlowExpo
      coef = DisSysCompLeakData(CompNum)%FlowCoef

      IF(PDROP.GE.0.0d0) THEN
         coef = coef/SQRTDZ(N)
      ELSE
         coef = Coef/SQRTDZ(M)
      END IF

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          ctl = (RhozNorm/RHOZ(N))**(expn-1.0d0)*(VisczNorm/VISCZ(N))**(2.0d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(N)/VISCZ(N)*ctl
        ELSE
          ctl = (RhozNorm/RHOZ(M))**(expn-1.0d0)*(VisczNorm/VISCZ(M))**(2.0d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(M)/VISCZ(M)*ctl
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction for laminar flow.
          ctl = (RhozNorm/RHOZ(N))**(expn-1.0d0)*(VisczNorm/VISCZ(N))**(2.0d0*expn-1.0d0)
          CDM = Coef*RHOZ(N)/VISCZ(N)*ctl
          FL = CDM*PDROP
          ! Flow in positive direction for turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = Coef*SQRTDZ(N)*SQRT(PDROP)
          ELSE
            FT = Coef*SQRTDZ(N)*(PDROP**expn)
          END IF
        ELSE
          ! Flow in negative direction for laminar flow
          CDM = Coef*RHOZ(M)/VISCZ(M)
          FL = CDM*PDROP
          ! Flow in negative direction for turbulent flow
          IF(expn .EQ. 0.5d0) THEN
            FT = -Coef*SQRTDZ(M)*SQRT(-PDROP)
          ELSE
            FT = -Coef*SQRTDZ(M)*(-PDROP)**expn
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(LIST.GE.4) WRITE(Unit21,901) ' plr: ',I,PDROP,FL,FT
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = FT*expn/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,I3,6X,4E16.7)
      RETURN
      END SUBROUTINE AFEPLR

      SUBROUTINE AFESCR(J,LFLAG,PDROP,I,N,M,F,DF,NF)
!
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a surface crack component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     CDM, FL, FT
         REAL(r64)     RhozNorm, VisczNorm, expn, Ctl, coef, Corr, VisAve, Tave, RhoCor
         INTEGER  CompNum

      ! FLOW:
      ! Crack standard condition from given inputs
      Corr = MultizoneSurfaceData(I)%Factor
      CompNum = AirflowNetworkCompData(J)%TypeNum
      RhozNorm = PsyRhoAirFnPbTdbW(MultizoneSurfaceCrackData(CompNum)%StandardP, &
                 MultizoneSurfaceCrackData(CompNum)%StandardT,MultizoneSurfaceCrackData(CompNum)%StandardW)
      VisczNorm = 1.71432d-5+4.828d-8*MultizoneSurfaceCrackData(CompNum)%StandardT

      expn = MultizoneSurfaceCrackData(CompNum)%FlowExpo
      VisAve = (VISCZ(N)+VISCZ(M))/2.0d0
      Tave = (TZ(N)+TZ(M))/2.0d0
      IF(PDROP.GE.0.0d0) THEN
         coef = MultizoneSurfaceCrackData(CompNum)%FlowCoef/SQRTDZ(N)*Corr
      ELSE
         coef = MultizoneSurfaceCrackData(CompNum)%FlowCoef/SQRTDZ(M)*Corr
      END IF

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          RhoCor = (TZ(N)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(N)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(N)/VISCZ(N)*ctl
        ELSE
          RhoCor = (TZ(M)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(M)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(M)/VISCZ(M)*ctl
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow.
          RhoCor = (TZ(N)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(N)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          CDM = Coef*RHOZ(N)/VISCZ(N)*ctl
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = Coef*SQRTDZ(N)*SQRT(PDROP)*ctl
          ELSE
            FT = Coef*SQRTDZ(N)*(PDROP**expn)*ctl
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow.
          RhoCor = (TZ(M)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(M)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          CDM = Coef*RHOZ(M)/VISCZ(M)*ctl
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = -Coef*SQRTDZ(M)*SQRT(-PDROP)*ctl
          ELSE
            FT = -Coef*SQRTDZ(M)*(-PDROP)**expn*ctl
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(LIST.GE.4) WRITE(Unit21,901) ' scr: ',I,PDROP,FL,FT
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = FT*expn/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,I3,6X,4E16.7)
      RETURN
      END SUBROUTINE AFESCR


      SUBROUTINE AFEDWC(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a duct/pipe component using Colebrook equation for the
          ! turbulent friction factor

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          !     RE      - Reynolds number.
          !     FL      - friction factor for laminar flow.
          !     FT      - friction factor for turbulent flow.

         REAL(r64)     A0, A1, A2, B, C, D, EPS, S2, CDM, FL, FT, FTT, RE
         INTEGER  CompNum
         REAL(r64)     ed, ld, g, AA1
!
      DATA C,EPS / 0.868589d0, 0.001d0 /

      ! FLOW:
      CompNum = AirflowNetworkCompData(J)%TypeNum
      ed = DisSysCompDuctData(CompNum)%Rough/DisSysCompDuctData(CompNum)%D
      ld = DisSysCompDuctData(CompNum)%L/DisSysCompDuctData(CompNum)%D
      g  = 1.14d0 - 0.868589d0*LOG(ed)
      AA1 = g

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = (2.d0*RHOZ(N)*DisSysCompDuctData(CompNum)%A*DisSysCompDuctData(CompNum)%D)/ &
                  (VISCZ(N)*DisSysCompDuctData(CompNum)%InitLamCoef*ld)
        ELSE
          DF(1) = (2.d0*RHOZ(M)*DisSysCompDuctData(CompNum)%A*DisSysCompDuctData(CompNum)%D)/ &
                  (VISCZ(M)*DisSysCompDuctData(CompNum)%InitLamCoef*ld)
        END IF
        F(1) = -DF(1)*PDROP
        IF(LIST.GE.4) WRITE(Unit21,901) ' dwi:',I,DisSysCompDuctData(CompNum)%InitLamCoef,F(1),DF(1)
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow coefficient !=0
          IF(DisSysCompDuctData(CompNum)%LamFriCoef.GE.0.001d0) THEN
            A2 = DisSysCompDuctData(CompNum)%LamFriCoef/(2.0d0*RHOZ(N)*DisSysCompDuctData(CompNum)%A* &
                 DisSysCompDuctData(CompNum)%A)
            A1 = (VISCZ(N)*DisSysCompDuctData(CompNum)%LamDynCoef*ld)/ &
                 (2.0d0*RHOZ(N)*DisSysCompDuctData(CompNum)%A*DisSysCompDuctData(CompNum)%D)
            A0 = -PDROP
            CDM = SQRT(A1*A1-4.0d0*A2*A0)
            FL = (CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(N)*DisSysCompDuctData(CompNum)%A*DisSysCompDuctData(CompNum)%D)/ &
                  (VISCZ(N)*DisSysCompDuctData(CompNum)%LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = FL*DisSysCompDuctData(CompNum)%D/(VISCZ(N)*DisSysCompDuctData(CompNum)%A)
          IF(LIST.GE.4) WRITE(Unit21,901) ' dwl:',I,PDROP,FL,CDM,RE
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(2.d0*RHOZ(N)*PDROP)*DisSysCompDuctData(CompNum)%A
            FTT = S2 / SQRT(ld/g**2+DisSysCompDuctData(CompNum)%TurDynCoef)
            IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,S2,FTT,g
            DO
              FT = FTT
              B = (9.3d0*VISCZ(N)*DisSysCompDuctData(CompNum)%A)/(FT*DisSysCompDuctData(CompNum)%Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+DisSysCompDuctData(CompNum)%TurDynCoef)
              IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,B,FTT,g
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = FTT
          ELSE
            FT = FL
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow coefficient !=0
          IF(DisSysCompDuctData(CompNum)%LamFriCoef.GE.0.001d0) THEN
            A2 = DisSysCompDuctData(CompNum)%LamFriCoef/(2.d0*RHOZ(M)*DisSysCompDuctData(CompNum)%A* &
                 DisSysCompDuctData(CompNum)%A)
            A1 = (VISCZ(M)*DisSysCompDuctData(CompNum)%LamDynCoef*ld)/ &
                 (2.d0*RHOZ(M)*DisSysCompDuctData(CompNum)%A*DisSysCompDuctData(CompNum)%D)
            A0 = PDROP
            CDM = SQRT(A1*A1-4.0d0*A2*A0)
            FL = -(CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(M)*DisSysCompDuctData(CompNum)%A*DisSysCompDuctData(CompNum)%D)/ &
                  (VISCZ(M)*DisSysCompDuctData(CompNum)%LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = -FL*DisSysCompDuctData(CompNum)%D/(VISCZ(M)*DisSysCompDuctData(CompNum)%A)
          IF(LIST.GE.4) WRITE(Unit21,901) ' dwl:',I,PDROP,FL,CDM,RE
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(-2.d0*RHOZ(M)*PDROP)*DisSysCompDuctData(CompNum)%A
            FTT = S2 / SQRT(ld/g**2+DisSysCompDuctData(CompNum)%TurDynCoef)
            IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,S2,FTT,g
            DO
              FT = FTT
              B = (9.3d0*VISCZ(M)*DisSysCompDuctData(CompNum)%A)/(FT*DisSysCompDuctData(CompNum)%Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+ DisSysCompDuctData(CompNum)%TurDynCoef)
              IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,B,FTT,g
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = -FTT
          ELSE
            FT = FL
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = 0.5d0*FT/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,I3,6X,4E16.7)

      RETURN
      END SUBROUTINE AFEDWC

      SUBROUTINE AFESOP(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a Doorway airflow component using standard interface.
          ! A doorway may have two-way airflows. Heights measured relative to the bottom of the door.

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          !     DPMID   - pressure drop at mid-height of doorway.
          !     DRHO    - difference in air densities between rooms.
          !     Y       - height of neutral plane rel. to bottom of door (m).
          !     F0      - flow factor at the bottom of the door.
          !     FH      - flow factor at the top of the door.
          !     DF0     - derivative factor at the bottom of the door.
          !     DFH     - derivative factor at the top of the door.
!
         REAL(r64) DPMID
         REAL(r64)     C, DF0, DFH, DRHO, GDRHO, F0, FH, Y, SQRT2
         INTEGER  CompNum
         REAL(r64) FlowCoef,FlowExpo,MinRhoDiff,DischCoeff,Width,Height,OpenFactor
!
      DATA SQRT2 / 1.414214d0 /

      ! FLOW:
      CompNum = AirflowNetworkCompData(J)%TypeNum
      MinRhoDiff = MultizoneCompSimpleOpeningData(CompNum)%MinRhoDiff
      DischCoeff = MultizoneCompSimpleOpeningData(CompNum)%DischCoeff
      Width = MultizoneSurfaceData(i)%Width
      Height = MultizoneSurfaceData(i)%Height
      FlowCoef = MultizoneCompSimpleOpeningData(CompNum)%FlowCoef*2.0d0*(Width+Height)
      FlowExpo = MultizoneCompSimpleOpeningData(CompNum)%FlowExpo
      OpenFactor = MultizoneSurfaceData(I)%OpenFactor
      If (OpenFactor > 0.0d0) then
        Width = Width*OpenFactor
        If (Surface(MultizoneSurfaceData(I)%SurfNum)%Tilt .LT. 90.d0) then
          Height = Height*Surface(MultizoneSurfaceData(I)%SurfNum)%SinTilt
        End If
      End If

      IF(PDROP.GE.0.0d0) THEN
         FlowCoef = FlowCoef/SQRTDZ(N)
      ELSE
         FlowCoef = FlowCoef/SQRTDZ(M)
      END IF

      ! Add window multiplier with window close
      If (MultizoneSurfaceData(I)%Multiplier > 1.0d0) FlowCoef=FlowCoef*MultizoneSurfaceData(I)%Multiplier
      ! Add window multiplier with window open
      If (OpenFactor .gt. 0.0d0) then
        If (MultizoneSurfaceData(I)%Multiplier > 1.0d0) Width=Width*MultizoneSurfaceData(I)%Multiplier
      End If

      NF = 1
      DRHO = RHOZ(N)-RHOZ(M)
      GDRHO = 9.8d0*DRHO
      IF(LIST.GE.4) WRITE(Unit21,903) ' DOR:',I,N,M,PDROP,ABS(DRHO),MinRhoDiff
      IF(OpenFactor .eq. 0.0d0) THEN
        CALL GenericCrack(FlowCoef,FlowExpo,LFLAG,PDROP,N,M,F,DF,NF)
        Return
      end if
      IF(ABS(DRHO).LT.MinRhoDiff .OR. LFLAG.EQ.1) THEN
        DPMID = PDROP-0.5d0*Height*GDRHO
        ! Initialization or identical temps: treat as one-way flow.
        CALL GenericCrack(FlowCoef,FlowExpo,LFLAG,DPMID,N,M,F,DF,NF)
        IF(LIST.GE.4) WRITE(Unit21,900) ' Drs:',DPMID,F(1),DF(1)
      ELSE
        ! Possible two-way flow:
        Y = PDROP/GDRHO
        IF(LIST.GE.4) WRITE(Unit21,900) ' DrY:',PDROP,GDRHO,Y
        ! F0 = lower flow, FH = upper flow.
        C = SQRT2*Width*DischCoeff
        DF0 = C*SQRT(ABS(PDROP))/ABS(GDRHO)
!        F0 = 0.666667d0*C*SQRT(ABS(GDRHO*Y))*ABS(Y)
        F0 = (2.0d0/3.0d0)*C*SQRT(ABS(GDRHO*Y))*ABS(Y)
        DFH = C*SQRT(ABS((Height-Y)/GDRHO))
!        FH = 0.666667d0*DFH*ABS(GDRHO*(Height-Y))
        FH = (2.0d0/3.0d0)*DFH*ABS(GDRHO*(Height-Y))
        IF(LIST.GE.4) WRITE(Unit21,900) ' DrF:',F0,DF0,FH,DFH
        IF(Y.LE.0.0d0) THEN
          ! One-way flow (negative).
          IF(DRHO.GE.0.0d0) THEN
            F(1) = -SQRTDZ(M)*ABS(FH-F0)
            DF(1) = SQRTDZ(M)*ABS(DFH-DF0)
          ELSE
            F(1) =  SQRTDZ(N)*ABS(FH-F0)
            DF(1) = SQRTDZ(N)*ABS(DFH-DF0)
          END IF
          IF(LIST.GE.4) WRITE(Unit21,900) ' Dr1:',C,F(1),DF(1)
        ELSE IF(Y.GE.Height) THEN
          ! One-way flow (positive).
          IF(DRHO.GE.0.0d0) THEN
            F(1) =  SQRTDZ(N)*ABS(FH-F0)
            DF(1) = SQRTDZ(N)*ABS(DFH-DF0)
          ELSE
            F(1) = -SQRTDZ(M)*ABS(FH-F0)
            DF(1) = SQRTDZ(M)*ABS(DFH-DF0)
          END IF
          IF(LIST.GE.4) WRITE(Unit21,900) ' Dr2:',C,F(1),DF(1)
        ELSE
          ! Two-way flow.
          NF = 2
          IF(DRHO.GE.0.0d0) THEN
            F(1) = -SQRTDZ(M)*FH
            DF(1) = SQRTDZ(M)*DFH
            F(2) =  SQRTDZ(N)*F0
            DF(2) = SQRTDZ(N)*DF0
          ELSE
            F(1) =  SQRTDZ(N)*FH
            DF(1) = SQRTDZ(N)*DFH
            F(2) = -SQRTDZ(M)*F0
            DF(2) = SQRTDZ(M)*DF0
          END IF
          IF(LIST.GE.4) WRITE(Unit21,900) ' Dr3:',C,F(1),DF(1)
          IF(LIST.GE.4) WRITE(Unit21,900) ' Dr4:',C,F(2),DF(2)
        ENDIF
      END IF
  900 FORMAT(A5,9X,4E16.7)
  903 FORMAT(A5,3I3,4E16.7)

      RETURN
      END SUBROUTINE AFESOP


      SUBROUTINE AFECFR(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a constant flow rate airflow component -- using standard interface.

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          USE DataLoopNode, ONLY: Node
          USE DataAirLoop,  ONLY: LoopSystemOnMassFlowrate,LoopSystemOffMassFlowrate,LoopFanOperationMode,LoopCompCycRatio
          USE DataHVACGlobals, ONLY: FanType_SimpleOnOff, FanType_SimpleConstVolume, FanType_SimpleVAV

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          INTEGER, PARAMETER :: CycFanCycComp    = 1 ! fan cycles with compressor operation


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          INTEGER CompNum, k,k1
          REAL(r64) :: SumTermFlow  ! Sum of all Terminal flows [kg/s]
          REAL(r64) :: SumFracSuppLeak  ! Sum of all supply leaks as a fraction of supply fan flow rate
          Integer Node1, Node2

      ! FLOW:
      CompNum = AirflowNetworkCompData(J)%TypeNum

      NF = 1
      If (DisSysCompCVFData(CompNum)%FanTypeNum .eq. FanType_SimpleOnOff) then
        If (LoopFanOperationMode .EQ. CycFanCycComp .and. LoopSystemOnMassFlowrate .GT. 0.0d0) then
          F(1) = LoopSystemOnMassFlowrate
        else
          F(1) = Node(DisSysCompCVFData(CompNum)%InletNode)%MassFlowrate*DisSysCompCVFData(CompNum)%Ctrl
          If (MultiSpeedHPIndicator .EQ. 2) Then
            F(1) = LoopSystemOnMassFlowrate*LoopCompCycRatio+LoopSystemOffMassFlowrate*(1.0d0-LoopCompCycRatio)
          End If
        End If
      Else If (DisSysCompCVFData(CompNum)%FanTypeNum .eq. FanType_SimpleConstVolume) Then
        If (DisSysCompCVFData(CompNum)%FlowRate > 0) Then
          F(1) = DisSysCompCVFData(CompNum)%FlowRate*DisSysCompCVFData(CompNum)%Ctrl
        Else
          F(1) = Node(DisSysCompCVFData(CompNum)%InletNode)%MassFlowrate*DisSysCompCVFData(CompNum)%Ctrl
        End If
        If (MultiSpeedHPIndicator .EQ. 2) Then
          F(1) = LoopSystemOnMassFlowrate
        End If
      Else If (DisSysCompCVFData(CompNum)%FanTypeNum .eq. FanType_SimpleVAV) Then
        ! Check VAV termals with a damper
        SumTermFlow = 0.d0
        SumFracSuppLeak = 0.d0
        Do k=1,NetworkNumOfLinks
          If (AirflowNetworkLinkageData(k)%VAVTermDamper) Then
            k1 = AirflowNetworkNodeData(AirflowNetworkLinkageData(k)%NodeNums(1))%EPlusNodeNum
            If (Node(k1)%MassFlowRate .GT. 0.0d0) Then
              SumTermFlow = SumTermFlow + Node(k1)%MassFlowRate
            End If
          End If
          If (AirflowNetworkCompData(AirflowNetworkLinkageData(k)%CompNum)%CompTypeNum == CompTypeNum_ELR) then
            ! Calculate supply leak sensible losses
            Node1 = AirflowNetworkLinkageData(k)%NodeNums(1)
            Node2 = AirflowNetworkLinkageData(k)%NodeNums(2)
            if ((AirflowNetworkNodeData(Node2)%EPlusZoneNum > 0) .AND. (AirflowNetworkNodeData(Node1)%EPlusNodeNum == 0)) Then
              SumFracSuppLeak = SumFracSuppLeak + &
                DisSysCompELRData(AirflowNetworkCompData(AirflowNetworkLinkageData(k)%CompNum)%TypeNum)%ELR
            End If
          End If
        End Do
        F(1) = SumTermFlow/(1.d0-SumFracSuppLeak)
        VAVTerminalRatio = 0.d0
        If (F(1) .GT. DisSysCompCVFData(CompNum)%MaxAirMassFlowRate) Then
          VAVTerminalRatio = DisSysCompCVFData(CompNum)%MaxAirMassFlowRate/F(1)
          F(1) = DisSysCompCVFData(CompNum)%MaxAirMassFlowRate
        End If
      End If
      DF(1) = 0.0d0

      RETURN
      END SUBROUTINE AFECFR


      SUBROUTINE AFEFAN(JA,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a detailed fan component -- using standard interface.

          ! 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)  :: JA    ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          !     PRISE   - pressure rise (negative of pressure drop) (Pa).
!
         INTEGER  J, K, L
         REAL(r64)     DPDF, PRISE, BX, BY, CX, CY, CCY, DX, DY, TOL
         INTEGER  CompNum, NumCur
         REAL(r64)     FlowCoef, FlowExpo

      DATA TOL / 0.00001d0 /

      ! FLOW:
      CompNum = AirflowNetworkCompData(JA)%TypeNum
      NumCur = DisSysCompDetFanData(CompNum)%n
      FlowCoef = DisSysCompDetFanData(CompNum)%FlowCoef
      FlowExpo = DisSysCompDetFanData(CompNum)%FlowExpo

      NF = 1
      IF(AFECTL(I).LE.0.0d0) THEN
        ! Speed = 0; treat fan as resistance.
        CALL GenericCrack(FlowCoef,FlowExpo,LFLAG,PDROP,N,M,F,DF,NF)
        RETURN
      END IF
      ! Pressure rise at reference fan speed.
      IF(AFECTL(I).GE.DisSysCompDetFanData(CompNum)%TranRat) THEN
        PRISE = -PDROP*(DisSysCompDetFanData(CompNum)%RhoAir/RHOZ(N))/AFECTL(I)**2
      ELSE
        PRISE = -PDROP*(DisSysCompDetFanData(CompNum)%RhoAir/RHOZ(N))/ &
                (DisSysCompDetFanData(CompNum)%TranRat*AFECTL(I))
      END IF
      IF(LIST.GE.4) WRITE(Unit21,901) ' fan:',I,PDROP,PRISE,AFECTL(I),DisSysCompDetFanData(CompNum)%TranRat
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear approximation.
        F(1) = -DisSysCompDetFanData(CompNum)%Qfree*AFECTL(I)*(1.0d0-PRISE/DisSysCompDetFanData(CompNum)%Pshut)
        DPDF = -DisSysCompDetFanData(CompNum)%Pshut/DisSysCompDetFanData(CompNum)%Qfree
        IF(LIST.GE.4) WRITE(Unit21,901) ' fni:',JA,DisSysCompDetFanData(CompNum)%Qfree, &
                                               DisSysCompDetFanData(CompNum)%Pshut
      ELSE
        ! Solution of the fan performance curve.
        ! Determine curve fit range.
        J = 1
        K = 5*(J-1)+1
        BX = DisSysCompDetFanData(CompNum)%Coeff(K)
        BY = DisSysCompDetFanData(CompNum)%Coeff(K+1)+BX*(DisSysCompDetFanData(CompNum)%Coeff(K+2)+ &
          BX*(DisSysCompDetFanData(CompNum)%Coeff(K+3)+BX*DisSysCompDetFanData(CompNum)%Coeff(K+4)))-PRISE
        IF(BY.LT.0.0d0) CALL ShowFatalError('Out of range, too low in an AirflowNetwork detailed Fan')

        DO
          DX = DisSysCompDetFanData(CompNum)%Coeff(K+5)
          DY = DisSysCompDetFanData(CompNum)%Coeff(K+1)+DX*(DisSysCompDetFanData(CompNum)%Coeff(K+2)+ &
               DX*(DisSysCompDetFanData(CompNum)%Coeff(K+3)+DX*DisSysCompDetFanData(CompNum)%Coeff(K+5)))-PRISE
          IF(LIST.GE.4) WRITE(Unit21,901) ' fp0:',J,BX,BY,DX,DY
          IF(BY*DY.LE.0.0d0) EXIT
          J = J+1
          IF(J.GT.NumCur) &
            CALL ShowFatalError('Out of range, too high (FAN) in ADS simulation')
          K = K+5
          BX = DX
          BY = DY
        END DO
        ! Determine reference mass flow rate by false position method.
        L = 0
        CY = 0.0d0
   40     CONTINUE
          L = L+1
          IF(L.GT.100) CALL ShowFatalError('Too many iterations (FAN) in AirflowNtework simulation')
          CCY = CY
          CX = BX-BY*((DX-BX)/(DY-BY))
          CY = DisSysCompDetFanData(CompNum)%Coeff(K+1)+CX*(DisSysCompDetFanData(CompNum)%Coeff(K+2)+ &
             CX*(DisSysCompDetFanData(CompNum)%Coeff(K+3)+CX*DisSysCompDetFanData(CompNum)%Coeff(K+4)))-PRISE
          IF(BY*CY .EQ. 0.0d0) GOTO 90
          IF(BY*CY .GT. 0.0d0) GOTO 60
   50     DX = CX
          DY = CY
          IF(CY*CCY.GT.0.0d0) BY = 0.5d0*BY
          GO TO 70
   60     BX = CX
          BY = CY
          IF(CY*CCY.GT.0.0d0) DY = 0.5d0*DY
   70     CONTINUE
          IF(LIST.GE.4) WRITE(Unit21,901) ' fpi:',J,BX,CX,DX,BY,DY
          IF(DX-BX.LT.TOL*CX) GO TO 80
          IF(DX-BX.LT.TOL) GO TO 80
          GO TO 40
   80   CX = 0.5d0*(BX+DX)
   90   F(1) = CX
        DPDF = DisSysCompDetFanData(CompNum)%Coeff(K+2)+CX*(2.0d0*DisSysCompDetFanData(CompNum)%Coeff(K+3)+ &
               CX*3.0d0*DisSysCompDetFanData(CompNum)%Coeff(K+4))
      END IF
      ! Convert to flow at given speed.
      F(1) = F(1)*(RHOZ(N)/DisSysCompDetFanData(CompNum)%RhoAir)*AFECTL(I)
      ! Set derivative w/r pressure drop (-).
      IF(AFECTL(I).GE.DisSysCompDetFanData(CompNum)%TranRat) THEN
        DF(1) = -AFECTL(I)/DPDF
      ELSE
        DF(1) = -1.0d0/DPDF
      END IF
  901 FORMAT(A5,I3,5E14.6)

  999 RETURN
      END SUBROUTINE AFEFAN

      ! The above subroutine is not used. Leave it for the time being and revise later.

      SUBROUTINE AFECPF(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a constant power simple fan component -- using standard interface.

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          ! na

      ! FLOW:
      NF = 1
      IF(LFLAG.EQ.1) THEN
        F(1) = AFECTL(I)
        DF(1) = F(1)
      ELSE
        F(1) = -AFECTL(I)/PDROP
        DF(1) = -F(1)/PDROP
      END IF

      RETURN
      END SUBROUTINE AFECPF


! Leave it for the time being and revise later. Or drop this component ???????????

      SUBROUTINE AFEDMP(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a Controlled power law resistance airflow component (damper)

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     C
         INTEGER  CompNum

      ! FLOW:
      ! Get component number
      CompNum = AirflowNetworkCompData(J)%TypeNum

      NF = 1
      C = AFECTL(I)
      IF(C.LT.DisSysCompDamperData(CompNum)%FlowMin) C = DisSysCompDamperData(CompNum)%FlowMin
      IF(C.GT.DisSysCompDamperData(CompNum)%FlowMax) C = DisSysCompDamperData(CompNum)%FlowMax
      C = DisSysCompDamperData(CompNum)%A0+C*(DisSysCompDamperData(CompNum)%A1+C* &
          (DisSysCompDamperData(CompNum)%A2+C*DisSysCompDamperData(CompNum)%A3))
      IF(LIST.GE.4) WRITE(Unit21,901) ' Dmp:',I,AFECTL(I),DisSysCompDamperData(CompNum)%FlowMin, &
                                      DisSysCompDamperData(CompNum)%FlowMax,C
      IF(LFLAG.EQ.1 .OR. ABS(PDROP).LE.DisSysCompDamperData(CompNum)%LTP) THEN
!                              Laminar flow.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = C*DisSysCompDamperData(CompNum)%LamFlow*RHOZ(N)/VISCZ(N)
        ELSE
          DF(1) = C*DisSysCompDamperData(CompNum)%LamFlow*RHOZ(M)/VISCZ(M)
        END IF
        F(1) = DF(1)*PDROP
      ELSE
!                              Turbulent flow.
        IF(PDROP.GE.0.0d0) THEN
          F(1) = C*DisSysCompDamperData(CompNum)%TurFlow*SQRTDZ(N)*PDROP**DisSysCompDamperData(CompNum)%FlowExpo
        ELSE
          F(1) = -C*DisSysCompDamperData(CompNum)%TurFlow*SQRTDZ(M)*(-PDROP)**DisSysCompDamperData(CompNum)%FlowExpo
        ENDIF
        DF(1) = F(1)*DisSysCompDamperData(CompNum)%FlowExpo/PDROP
      END IF
  901 FORMAT(A5,I3,6X,4E16.7)

      RETURN
      END SUBROUTINE AFEDMP


      SUBROUTINE AFESEL(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a Surface effective leakage area component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     CDM, FL, FT
         REAL(r64) FlowCoef,FlowExpo
         INTEGER CompNum

      ! FLOW:
      ! Get component properties
      CompNum  = AirflowNetworkCompData(J)%TypeNum
      FlowExpo = MultizoneSurfaceELAData(CompNum)%FlowExpo
      FlowCoef = MultizoneSurfaceELAData(CompNum)%ELA*MultizoneSurfaceELAData(CompNum)%DischCoeff*(2)**0.5d0* &
                 MultizoneSurfaceELAData(CompNum)%RefDeltaP**(0.5d0-FlowExpo)

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = FlowCoef*RHOZ(N)/VISCZ(N)
        ELSE
          DF(1) = FlowCoef*RHOZ(M)/VISCZ(M)
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow.
          CDM = FlowCoef*RHOZ(N)/VISCZ(N)
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(MultizoneSurfaceELAData(CompNum)%FlowExpo.EQ.0.5d0) THEN
            FT = FlowCoef*SQRTDZ(N)*SQRT(PDROP)
          ELSE
            FT = FlowCoef*SQRTDZ(N)*(PDROP**FlowExpo)
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow.
          CDM = FlowCoef*RHOZ(M)/VISCZ(M)
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(FlowExpo .EQ.0.5d0) THEN
            FT = -FlowCoef*SQRTDZ(M)*SQRT(-PDROP)
          ELSE
            FT = -FlowCoef*SQRTDZ(M)*(-PDROP)**FlowExpo
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(LIST.GE.4) WRITE(Unit21,901) ' plr: ',I,PDROP,FL,FT
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = FT*FlowExpo/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,I3,6X,4E16.7)

      RETURN
      END SUBROUTINE AFESEL


      SUBROUTINE AFEELR(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a Effective leakage ratio component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     CDM, FL, FT
         REAL(r64) FlowCoef
         INTEGER CompNum

      ! FLOW:
      ! Get component properties
      CompNum  = AirflowNetworkCompData(J)%TypeNum
      FlowCoef = DisSysCompELRData(CompNum)%ELR*DisSysCompELRData(CompNum)%FlowRate/RHOZ(N)* &
                 DisSysCompELRData(CompNum)%RefPres**(-DisSysCompELRData(CompNum)%FlowExpo)

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = FlowCoef*RHOZ(N)/VISCZ(N)
        ELSE
          DF(1) = FlowCoef*RHOZ(M)/VISCZ(M)
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow.
          CDM = FlowCoef*RHOZ(N)/VISCZ(N)
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(DisSysCompELRData(CompNum)%FlowExpo.EQ.0.5d0) THEN
            FT = FlowCoef*SQRTDZ(N)*SQRT(PDROP)
          ELSE
            FT = FlowCoef*SQRTDZ(N)*(PDROP**DisSysCompELRData(CompNum)%FlowExpo)
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow.
          CDM = FlowCoef*RHOZ(M)/VISCZ(M)
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(DisSysCompELRData(CompNum)%FlowExpo.EQ.0.5d0) THEN
            FT = -FlowCoef*SQRTDZ(M)*SQRT(-PDROP)
          ELSE
            FT = -FlowCoef*SQRTDZ(M)*(-PDROP)**DisSysCompELRData(CompNum)%FlowExpo
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(LIST.GE.4) WRITE(Unit21,901) ' plr: ',I,PDROP,FL,FT
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = FT*DisSysCompELRData(CompNum)%FlowExpo/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,I3,6X,4E16.7)

      RETURN
      END SUBROUTINE AFEELR


      SUBROUTINE AFECPD(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a Constant pressure drop component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(INOUT) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     Co
         INTEGER  CompNum,K

      ! FLOW:
      ! Get component properties
      ! A  = Cross section area [m2]
      ! DP = Pressure difference across the element [Pa]
      CompNum  = AirflowNetworkCompData(J)%TypeNum

      NF = 1
      if (PDROP .EQ. 0.0d0) then
         F(1) = SQRT(2.0d0*RHOZ(N))*DisSysCompCPDData(CompNum)%A*DisSysCompCPDData(CompNum)%DP**0.5d0
         DF(1) = 0.5d0*F(1)/DisSysCompCPDData(CompNum)%DP
      else
         DO k=1,NetworkNumOfLinks
            if (AirflowNetworkLinkageData(k)%NodeNums(2) .EQ. N) then
               F(1) = AFLOW(k)
               Exit
            end if
         end do
         PDROP = DisSysCompCPDData(CompNum)%DP
         PZ(M) = PZ(N) -PDROP
         Co=F(1)/DisSysCompCPDData(CompNum)%DP
         DF(1) = 10.d10
      end if

      RETURN
      END SUBROUTINE AFECPD

      SUBROUTINE AFECOI(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a coil component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     A0, A1, A2, B, C, D, EPS, S2, CDM, FL, FT, FTT, RE
         INTEGER  CompNum
         REAL(r64)     ed, ld, g, AA1, rough,area,InitLamCoef,LamDynCoef,LamFriCoef,TurDynCoef

      DATA C,EPS / 0.868589d0, 0.001d0 /
      DATA rough/0.0001d0/
      DATA InitLamCoef,LamDynCoef,LamFriCoef,TurDynCoef/128.0d0,64.0d0,0.0001d0,0.0001d0/

      ! FLOW:
      ! Get component properties
      CompNum = AirflowNetworkCompData(J)%TypeNum
      ed = Rough/DisSysCompDuctData(CompNum)%D
      area = (DisSysCompCoilData(CompNum)%D)**2*PI
      ld = DisSysCompCoilData(CompNum)%L/DisSysCompCoilData(CompNum)%D
      g  = 1.14d0 - 0.868589d0*LOG(ed)
      AA1 = g

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = (2.d0*RHOZ(N)*area*DisSysCompCoilData(CompNum)%D)/(VISCZ(N)*InitLamCoef*ld)
        ELSE
          DF(1) = (2.d0*RHOZ(M)*Area*DisSysCompCoilData(CompNum)%D)/(VISCZ(M)*InitLamCoef*ld)
        END IF
        F(1) = -DF(1)*PDROP
        IF(LIST.GE.4) WRITE(Unit21,901) ' dwi:',I,InitLamCoef,F(1),DF(1)
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow coefficient !=0
          IF(LamFriCoef.GE.0.001d0) THEN
            A2 = LamFriCoef/(2.0d0*RHOZ(N)*Area*Area)
            A1 = (VISCZ(N)*LamDynCoef*ld)/(2.d0*RHOZ(N)*Area*DisSysCompCoilData(CompNum)%D)
            A0 = -PDROP
            CDM = SQRT(A1*A1-4.d0*A2*A0)
            FL = (CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(N)*Area*DisSysCompCoilData(CompNum)%D)/(VISCZ(N)*LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = FL*DisSysCompCoilData(CompNum)%D/(VISCZ(N)*Area)
          IF(LIST.GE.4) WRITE(Unit21,901) ' dwl:',I,PDROP,FL,CDM,RE
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(2.d0*RHOZ(N)*PDROP)*Area
            FTT = S2 / SQRT(ld/g**2+TurDynCoef)
            IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,S2,FTT,g
            DO
              FT = FTT
              B = (9.3d0*VISCZ(N)*Area)/(FT*Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+TurDynCoef)
              IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,B,FTT,g
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = FTT
          ELSE
            FT = FL
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow coefficient !=0
          IF(LamFriCoef.GE.0.001d0) THEN
            A2 = LamFriCoef/(2.d0*RHOZ(M)*Area*Area)
            A1 = (VISCZ(M)*LamDynCoef*ld)/(2.d0*RHOZ(M)*Area*DisSysCompCoilData(CompNum)%D)
            A0 = PDROP
            CDM = SQRT(A1*A1-4.d0*A2*A0)
            FL = -(CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(M)*Area*DisSysCompCoilData(CompNum)%D)/(VISCZ(M)*LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = -FL*DisSysCompCoilData(CompNum)%D/(VISCZ(M)*Area)
          IF(LIST.GE.4) WRITE(Unit21,901) ' dwl:',I,PDROP,FL,CDM,RE
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(-2.d0*RHOZ(M)*PDROP)*Area
            FTT = S2 / SQRT(ld/g**2+TurDynCoef)
            IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,S2,FTT,g
            DO
              FT = FTT
              B = (9.3d0*VISCZ(M)*Area)/(FT*Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+TurDynCoef)
              IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,B,FTT,g
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = -FTT
          ELSE
            FT = FL
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = 0.5d0*FT/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,I3,6X,4E16.7)

      RETURN
      END SUBROUTINE AFECOI

      SUBROUTINE AFETMU(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a terminal unit component

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          USE DataLoopNode, ONLY: Node

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     A0, A1, A2, B, C, D, EPS, S2, CDM, FL, FT, FTT, RE
         INTEGER  CompNum
         REAL(r64)     ed, ld, g, AA1, rough,Area,InitLamCoef,LamDynCoef,LamFriCoef,TurDynCoef

      DATA C,EPS / 0.868589d0, 0.001d0 /
      DATA rough/0.0001d0/
      DATA InitLamCoef,LamDynCoef,LamFriCoef,TurDynCoef/128.0d0,64.0d0,0.0001d0,0.0001d0/

      ! FLOW:
      ! Get component properties
      CompNum = AirflowNetworkCompData(J)%TypeNum
      ed = Rough/DisSysCompTermUnitData(CompNum)%D
      Area = (DisSysCompTermUnitData(CompNum)%D)**2*PI
      ld = DisSysCompTermUnitData(CompNum)%L/DisSysCompTermUnitData(CompNum)%D
      g  = 1.14d0 - 0.868589d0*LOG(ed)
      AA1 = g

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = (2.d0*RHOZ(N)*area*DisSysCompTermUnitData(CompNum)%D)/(VISCZ(N)*InitLamCoef*ld)
        ELSE
          DF(1) = (2.d0*RHOZ(M)*Area*DisSysCompTermUnitData(CompNum)%D)/(VISCZ(M)*InitLamCoef*ld)
        END IF
        F(1) = -DF(1)*PDROP
        IF(LIST.GE.4) WRITE(Unit21,901) ' dwi:',I,InitLamCoef,F(1),DF(1)
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow coefficient !=0
          IF(LamFriCoef.GE.0.001d0) THEN
            A2 = LamFriCoef/(2.d0*RHOZ(N)*Area*Area)
            A1 = (VISCZ(N)*LamDynCoef*ld)/(2.d0*RHOZ(N)*Area*DisSysCompTermUnitData(CompNum)%D)
            A0 = -PDROP
            CDM = SQRT(A1*A1-4.d0*A2*A0)
            FL = (CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(N)*Area*DisSysCompTermUnitData(CompNum)%D)/(VISCZ(N)*LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = FL*DisSysCompTermUnitData(CompNum)%D/(VISCZ(N)*Area)
          IF(LIST.GE.4) WRITE(Unit21,901) ' dwl:',I,PDROP,FL,CDM,RE
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(2.d0*RHOZ(N)*PDROP)*Area
            FTT = S2 / SQRT(ld/g**2+TurDynCoef)
            IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,S2,FTT,g
            DO
              FT = FTT
              B = (9.3d0*VISCZ(N)*Area)/(FT*Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+TurDynCoef)
              IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,B,FTT,g
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = FTT
          ELSE
            FT = FL
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow coefficient !=0
          IF(LamFriCoef.GE.0.001d0) THEN
            A2 = LamFriCoef/(2.d0*RHOZ(M)*Area*Area)
            A1 = (VISCZ(M)*LamDynCoef*ld)/(2.0d0*RHOZ(M)*Area*DisSysCompTermUnitData(CompNum)%D)
            A0 = PDROP
            CDM = SQRT(A1*A1-4.d0*A2*A0)
            FL = -(CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(M)*Area*DisSysCompTermUnitData(CompNum)%D)/(VISCZ(M)*LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = -FL*DisSysCompTermUnitData(CompNum)%D/(VISCZ(M)*Area)
          IF(LIST.GE.4) WRITE(Unit21,901) ' dwl:',I,PDROP,FL,CDM,RE
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(-2.d0*RHOZ(M)*PDROP)*Area
            FTT = S2 / SQRT(ld/g**2+TurDynCoef)
            IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,S2,FTT,g
            DO
              FT = FTT
              B = (9.3d0*VISCZ(M)*Area)/(FT*Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+TurDynCoef)
              IF(LIST.GE.4) WRITE(Unit21,901) ' dwt:',I,B,FTT,g
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = -FTT
          ELSE
            FT = FL
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = 0.5d0*FT/PDROP
        END IF
      END IF
      ! If damper, setup the airflows from nodal values calculated from teminal
      If (AirflowNetworkLinkageData(I)%VAVTermDamper) Then
        F(1) = Node(DisSysCompTermUnitData(CompNum)%DamperInletNode)%MassFlowrate
        If (VAVTerminalRatio .gt. 0.d0) Then
          F(1) = F(1)*VAVTerminalRatio
        End If
        DF(1) = 0.0d0
      End IF

!
  901 FORMAT(A5,I3,6X,4E16.7)

      RETURN
      END SUBROUTINE AFETMU

      SUBROUTINE AFEEXF(J,LFLAG,PDROP,I,N,M,F,DF,NF)
!
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 12/17/06
          !                      Revised for zone exhaust fan
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a surface crack component

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          USE DataLoopNode,    ONLY: Node
          USE DataHVACGlobals, ONLY: VerySmallMassFlow

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     CDM, FL, FT
         REAL(r64)     RhozNorm, VisczNorm, expn, Ctl, coef, Corr, VisAve, Tave, RhoCor
         INTEGER  CompNum,InletNode

    ! FLOW:
    CompNum = AirflowNetworkCompData(J)%TypeNum
    InletNode = MultizoneCompExhaustFanData(CompNum)%InletNode
    If (Node(InletNode)%MassFlowRate > VerySmallMassFlow) then
      ! Treat the component as an exhaust fan
      NF = 1
      F(1) = Node(InletNode)%MassFlowRate
      DF(1) = 0.0d0
      RETURN
    Else
      ! Treat the component as a surface crack
      ! Crack standard condition from given inputs
      Corr = MultizoneSurfaceData(I)%Factor
      RhozNorm = PsyRhoAirFnPbTdbW(MultizoneCompExhaustFanData(CompNum)%StandardP, &
                 MultizoneCompExhaustFanData(CompNum)%StandardT,MultizoneCompExhaustFanData(CompNum)%StandardW)
      VisczNorm = 1.71432d-5+4.828d-8*MultizoneCompExhaustFanData(CompNum)%StandardT

      expn = MultizoneCompExhaustFanData(CompNum)%FlowExpo
      VisAve = (VISCZ(N)+VISCZ(M))/2.0d0
      Tave = (TZ(N)+TZ(M))/2.0d0
      IF(PDROP.GE.0.0d0) THEN
         coef = MultizoneCompExhaustFanData(CompNum)%FlowCoef/SQRTDZ(N)*Corr
      ELSE
         coef = MultizoneCompExhaustFanData(CompNum)%FlowCoef/SQRTDZ(M)*Corr
      END IF

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          RhoCor = (TZ(N)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(N)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(N)/VISCZ(N)*ctl
        ELSE
          RhoCor = (TZ(M)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(M)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(M)/VISCZ(M)*ctl
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow.
          RhoCor = (TZ(N)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(N)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          CDM = Coef*RHOZ(N)/VISCZ(N)*ctl
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = Coef*SQRTDZ(N)*SQRT(PDROP)*ctl
          ELSE
            FT = Coef*SQRTDZ(N)*(PDROP**expn)*ctl
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow.
          RhoCor = (TZ(M)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(M)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          CDM = Coef*RHOZ(M)/VISCZ(M)*ctl
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = -Coef*SQRTDZ(M)*SQRT(-PDROP)*ctl
          ELSE
            FT = -Coef*SQRTDZ(M)*(-PDROP)**expn*ctl
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(LIST.GE.4) WRITE(Unit21,901) ' scr: ',I,PDROP,FL,FT
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = FT*expn/PDROP
        END IF
      END IF
    End If
!
  901 FORMAT(A5,I3,6X,4E16.7)
      RETURN
      END SUBROUTINE AFEEXF

      SUBROUTINE AFEHEX(J,LFLAG,PDROP,I,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 1/18/09
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a heat exchanger component

          ! 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)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     A0, A1, A2, B, C, D, EPS, S2, CDM, FL, FT, FTT, RE
         INTEGER  CompNum
         REAL(r64)     ed, ld, g, AA1, rough,area,InitLamCoef,LamDynCoef,LamFriCoef,TurDynCoef

      DATA C,EPS / 0.868589d0, 0.001d0 /
      DATA rough/0.0001d0/
      DATA InitLamCoef,LamDynCoef,LamFriCoef,TurDynCoef/128.0d0,64.0d0,0.0001d0,0.0001d0/

      ! FLOW:
      ! Get component properties
      CompNum = AirflowNetworkCompData(J)%TypeNum
      ed = Rough/DisSysCompHXData(CompNum)%D
      area = (DisSysCompHXData(CompNum)%D)**2*PI
      ld = DisSysCompHXData(CompNum)%L/DisSysCompHXData(CompNum)%D
      g  = 1.14d0 - 0.868589d0*LOG(ed)
      AA1 = g

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          DF(1) = (2.d0*RHOZ(N)*area*DisSysCompHXData(CompNum)%D)/(VISCZ(N)*InitLamCoef*ld)
        ELSE
          DF(1) = (2.d0*RHOZ(M)*Area*DisSysCompHXData(CompNum)%D)/(VISCZ(M)*InitLamCoef*ld)
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow coefficient !=0
          IF(LamFriCoef.GE.0.001d0) THEN
            A2 = LamFriCoef/(2.0d0*RHOZ(N)*Area*Area)
            A1 = (VISCZ(N)*LamDynCoef*ld)/(2.d0*RHOZ(N)*Area*DisSysCompHXData(CompNum)%D)
            A0 = -PDROP
            CDM = SQRT(A1*A1-4.d0*A2*A0)
            FL = (CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(N)*Area*DisSysCompHXData(CompNum)%D)/(VISCZ(N)*LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = FL*DisSysCompHXData(CompNum)%D/(VISCZ(N)*Area)
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.0d0) THEN
            S2 = SQRT(2.d0*RHOZ(N)*PDROP)*Area
            FTT = S2 / SQRT(ld/g**2+TurDynCoef)
            DO
              FT = FTT
              B = (9.3d0*VISCZ(N)*Area)/(FT*Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+TurDynCoef)
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = FTT
          ELSE
            FT = FL
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow coefficient !=0
          IF(LamFriCoef.GE.0.001d0) THEN
            A2 = LamFriCoef/(2.d0*RHOZ(M)*Area*Area)
            A1 = (VISCZ(M)*LamDynCoef*ld)/(2.d0*RHOZ(M)*Area*DisSysCompHXData(CompNum)%D)
            A0 = PDROP
            CDM = SQRT(A1*A1-4.d0*A2*A0)
            FL = -(CDM-A1)/(2.d0*A2)
            CDM = 1.0d0/CDM
          ELSE
            CDM = (2.d0*RHOZ(M)*Area*DisSysCompHXData(CompNum)%D)/(VISCZ(M)*LamDynCoef*ld)
            FL = CDM*PDROP
          END IF
          RE = -FL*DisSysCompHXData(CompNum)%D/(VISCZ(M)*Area)
          ! Turbulent flow; test when Re>10.
          IF(RE.GE.10.d0) THEN
            S2 = SQRT(-2.d0*RHOZ(M)*PDROP)*Area
            FTT = S2 / SQRT(ld/g**2+TurDynCoef)
            DO
              FT = FTT
              B = (9.3d0*VISCZ(M)*Area)/(FT*Rough)
              D = 1.0d0 + g*B
              g = g - (g-AA1+C*LOG(D))/(1.0d0+C*B/D)
              FTT = S2 / SQRT(ld/g**2+TurDynCoef)
              IF(ABS(FTT-FT)/FTT .LT. EPS) EXIT
            ENDDO
            FT = -FTT
          ELSE
            FT = FL
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = 0.5d0*FT/PDROP
        END IF
      END IF

      RETURN
      END SUBROUTINE AFEHEX

      SUBROUTINE AFEHOP(J,LFLAG,PDROP,I,N,M,F,DF,NF)
!
          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Apr. 2009
          !       MODIFIED       na
          !       MODIFIED       na
          !
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a horizontal opening component. The subroutine was
          ! developed based on the subroutine AFEPLR of AIRNET.

          ! METHODOLOGY EMPLOYED:
          ! Combine forced and buyancy airflows together with a cap

          ! REFERENCES:
          ! Cooper, L., 1989, "Calculation of the Flow Through a Horizontal Ceiling/Floor Vent,"
          ! NISTIR 89-4052, National Institute of Standards and Technology, Gaithersburg, MD

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN)     :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: I     ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT)    :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT)    :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          REAL(r64), PARAMETER :: PI = 3.14159265358979323846d0


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     RhozAver, expn, coef
         INTEGER  CompNum
         REAL(r64)     Width            ! Opening width
         REAL(r64)     Height           ! Opening height
         REAL(r64)     Fact             ! Opening factor
         REAL(r64)     Slope            ! Opening slope
         REAL(r64)     DischCoeff       ! Discharge coefficient
         REAL(r64)     fma12            ! massflow in direction "from-to" [kg/s]
         REAL(r64)     fma21            ! massflow in direction "to-from" [kg/s]
         REAL(r64)     dp1fma12         ! derivative d fma12 / d Dp [kg/s/Pa]
         REAL(r64)     dp1fma21         ! derivative d fma21 / d Dp [kg/s/Pa]
         REAL(r64)     PurgedP          ! Purge pressure [Pa]
         REAL(r64)     BuoFlow          ! Buoyancy flow rate [Pa]
         REAL(r64)     BuoFlowMax       ! Maximum buoyancy flow rate [Pa]
         REAL(r64)     dPBuoFlow        ! Derivative of buoyancy flow rate [kg/s/Pa]
         REAL(r64)     DH               ! Hydraulic diameter [m]
         REAL(r64)     Cshape           ! Shape factor [dimensionless]
         REAL(r64)     OpenArea         ! Opening area [m2]


      ! FLOW:
      ! Get information on the horizontal opening
      CompNum = AirflowNetworkCompData(J)%TypeNum
      RhozAver = (RHOZ(N)+RHOZ(M))/2.0d0
      Width = MultizoneSurfaceData(I)%Width
      Height = MultizoneSurfaceData(I)%Height
      fact = MultizoneSurfaceData(I)%OpenFactor
      expn = MultizoneCompHorOpeningData(CompNum)%FlowExpo
      coef = MultizoneCompHorOpeningData(CompNum)%FlowCoef
      slope = MultizoneCompHorOpeningData(CompNum)%slope
      DischCoeff = MultizoneCompHorOpeningData(CompNum)%DischCoeff
      Cshape = 0.942d0*Width/Height
      OpenArea = Width*Height*fact*sin(slope*PI/180.0d0)*(1.0d0+cos(slope*PI/180.0d0))
      DH = 4.0d0*(Width*Height)/2.0/(Width+Height)*fact

      ! Check which zone is higher

      IF(fact .eq. 0.0d0) THEN
        CALL GenericCrack(coef,expn,LFLAG,PDROP,N,M,F,DF,NF)
        Return
      end if


      fma12     = 0.0d0
      fma21     = 0.0d0
      dp1fma12  = 0.0d0
      dp1fma21  = 0.0d0
      BuoFlow   = 0.0d0
      dPBuoFlow = 0.0d0

      If (AirflowNetworkLinkageData(I)%NodeHeights(1) > AirflowNetworkLinkageData(I)%NodeHeights(2)) Then
      ! Node N is upper zone
        If (RHOZ(N) > RHOZ(M)) Then
          BuoFlowMax = RhozAver*0.055d0*SQRT(9.81*abs(RHOZ(N)-RHOZ(M))*DH**5/RhozAver)
          PurgedP = Cshape*Cshape*9.81d0*ABS(RHOZ(N)-RHOZ(M))*DH**5/(2.0d0*(OpenArea)**2)
          If (abs(PDROP) .LE. PurgedP) Then
            BuoFlow = BuoFlowMax*(1.0d0-abs(PDROP)/PurgedP)
            dPBuoFlow = BuoFlowMax/PurgedP
          End If
        End If
      Else
      ! Node M is upper zone
        If (RHOZ(N) < RHOZ(M)) Then
          BuoFlowMax = RhozAver*0.055d0*SQRT(9.81d0*abs(RHOZ(N)-RHOZ(M))*DH**5/RhozAver)
          PurgedP = Cshape*Cshape*9.81d0*ABS(RHOZ(N)-RHOZ(M))*DH**5/(2.0d0*(OpenArea)**2)
          If (abs(PDROP) .LE. PurgedP) Then
            BuoFlow = BuoFlowMax*(1.0d0-abs(PDROP)/PurgedP)
            dPBuoFlow = BuoFlowMax/PurgedP
          End If
        End If
      End If

      IF(PDROP.EQ.0.0d0) THEN
         fma12 = BuoFlow
         fma21 = BuoFlow
         dp1fma12 = 0.0d0
         dp1fma21 = 0.0d0
      ELSE IF (PDROP.GT.0.0d0) THEN
         fma12 = RHOZ(N)*OpenArea*fact*DischCoeff*SQRT(2.0d0*PDROP/RhozAver)+BuoFlow
         dp1fma12 = RHOZ(N)*OpenArea*DischCoeff/SQRT(2.0d0*PDROP*RhozAver)+dPBuoFlow
         If (BuoFlow .gt. 0.0d0) Then
           fma21 = BuoFlow
           dp1fma21 = dPBuoFlow
         End If
      ELSE ! PDROP.LT.0.0
         fma21 = RHOZ(M)*OpenArea*fact*DischCoeff*SQRT(2.0d0*ABS(PDROP)/RhozAver)+BuoFlow
         dp1fma21 = -RHOZ(M)*OpenArea*DischCoeff/SQRT(2.0d0*ABS(PDROP)*RhozAver)+dPBuoFlow
         If (BuoFlow .gt. 0.0d0) Then
           fma12 = BuoFlow
           dp1fma12 = dPBuoFlow
         End If
      END IF

      F(1) = fma12-fma21
      DF(1) = dp1fma12-dp1fma21
      F(2) = 0.0d0
      if (fma12 .NE. 0.0d0 .and. fma21 .NE. 0.0d0) then
        F(2) = fma21
      End if
      DF(2) = 0.0d0

      RETURN
      END SUBROUTINE AFEHOP

      SUBROUTINE GenericCrack(Coef,Expn,LFLAG,PDROP,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  This subroutine is revised from AFEPLR developed by George Walton, NIST

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves airflow for a power law component

          ! METHODOLOGY EMPLOYED:
          ! Using Q=C(dP)^n

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          REAL(r64), INTENT(INOUT)  :: Coef     ! Flow coefficient
          REAL(r64), INTENT(IN)  :: Expn     ! Flow exponent
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         REAL(r64)     CDM, FL, FT
         REAL(r64)     RhozNorm, VisczNorm, Ctl, VisAve, Tave, RhoCor

      ! FLOW:
      ! Calculate normal density and viscocity at Crack standard condition: T=20C, p=101325 Pa and 0 g/kg
      RhozNorm = PsyRhoAirFnPbTdbW(101325.0d0,20.0d0,0.0d0)
      VisczNorm = 1.71432d-5+4.828d-8*20.0d0
      VisAve = (VISCZ(N)+VISCZ(M))/2.0d0
      Tave = (TZ(N)+TZ(M))/2.0d0
      IF(PDROP.GE.0.0d0) THEN
         coef = Coef/SQRTDZ(N)
      ELSE
         coef = Coef/SQRTDZ(M)
      END IF

      NF = 1
      IF(LFLAG.EQ.1) THEN
        ! Initialization by linear relation.
        IF(PDROP.GE.0.0d0) THEN
          RhoCor = (TZ(N)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(N)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(N)/VISCZ(N)*ctl
        ELSE
          RhoCor = (TZ(M)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(M)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          DF(1) = Coef*RHOZ(M)/VISCZ(M)*ctl
        END IF
        F(1) = -DF(1)*PDROP
      ELSE
        ! Standard calculation.
        IF(PDROP.GE.0.0d0) THEN
          ! Flow in positive direction.
          ! Laminar flow.
          RhoCor = (TZ(N)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(N)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          CDM = Coef*RHOZ(N)/VISCZ(N)*ctl
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = Coef*SQRTDZ(N)*SQRT(PDROP)*ctl
          ELSE
            FT = Coef*SQRTDZ(N)*(PDROP**expn)*ctl
          END IF
        ELSE
          ! Flow in negative direction.
          ! Laminar flow.
          RhoCor = (TZ(M)+KelvinConv)/(Tave+KelvinConv)
          ctl = (RhozNorm/RHOZ(M)/RhoCor)**(expn-1.0d0)*(VisczNorm/VisAve)**(2.d0*expn-1.0d0)
          CDM = Coef*RHOZ(M)/VISCZ(M)*ctl
          FL = CDM*PDROP
          ! Turbulent flow.
          IF(expn .EQ. 0.5d0) THEN
            FT = -Coef*SQRTDZ(M)*SQRT(-PDROP)*ctl
          ELSE
            FT = -Coef*SQRTDZ(M)*(-PDROP)**expn*ctl
          END IF
        END IF
        ! Select laminar or turbulent flow.
        IF(LIST.GE.4) WRITE(Unit21,901) ' generic crack: ',PDROP,FL,FT
        IF(ABS(FL).LE.ABS(FT)) THEN
          F(1) = FL
          DF(1) = CDM
        ELSE
          F(1) = FT
          DF(1) = FT*expn/PDROP
        END IF
      END IF
!
  901 FORMAT(A5,6X,4E16.7)
      RETURN
      END SUBROUTINE GenericCrack


      SUBROUTINE FACSKY(AU,AD,AL,IK,NEQ,NSYM)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  This subroutine is revised from FACSKY developed by George Walton, NIST

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine performs L-U factorization of a skyline ordered matrix, [A]
          ! The algorithm has been restructured for clarity.
          ! Note dependence on compiler for optimizing the inner do loops.

          ! METHODOLOGY EMPLOYED:
          !     L-U factorization of a skyline ordered matrix, [A], used for
          !     solution of simultaneous linear algebraic equations [A] * X = B.
          !     No pivoting!  No scaling!  No warnings!!!
          !     Related routines:  SLVSKY, SETSKY, FILSKY.

          ! REFERENCES:
          !     Algorithm is described in "The Finite Element Method Displayed",
          !     by G. Dhatt and G. Touzot, John Wiley & Sons, New York, 1984.

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: IK(NetworkNumOfNodes+1) ! pointer to the top of column/row "K"
          INTEGER, INTENT(IN)  :: NEQ   ! number of equations
          INTEGER, INTENT(IN)  :: NSYM  ! symmetry:  0 = symmetric matrix, 1 = non-symmetric
          ! noel, GNU says the AU is indexed above its upper bound
          !REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)-1) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AD(NetworkNumOfNodes)  ! the main diagonal of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AL(IK(NetworkNumOfNodes+1)-1) ! the lower triangle of [A] before and after factoring


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         INTEGER  JHK, JHK1, LHK, LHK1, IMIN, IMIN1
         INTEGER  JHJ, JHJ1, IC, I, J, K
         REAL(r64)  T1, T2, SDOT, SUMD

      ! FLOW:
      AD(1) = 1.0d0/AD(1)
      JHK = 1
      DO K=2,NEQ
        SUMD = 0.0d0
        JHK1 = IK(K+1)
        LHK = JHK1-JHK
        IF(LHK.GT.0) THEN
          LHK1 = LHK-1
          IMIN = K-LHK1
          IMIN1 = IMIN-1
          IF(NSYM.EQ.1) AL(JHK) = AL(JHK)*AD(IMIN1)
          IF(LHK1.NE.0) THEN
            JHJ = IK(IMIN)
            IF(NSYM.EQ.0) THEN
              DO J=1,LHK1
                JHJ1 = IK(IMIN+J)
                IC = MIN(J,JHJ1-JHJ)
                IF(IC.GT.0) THEN
                  SDOT = 0.0d0
                  DO I=0,IC-1
                    SDOT = SDOT+AU(JHJ1-IC+I)*AU(JHK+J-IC+I)
                  END DO
                  AU(JHK+J) = AU(JHK+J)-SDOT
                END IF
                JHJ = JHJ1
              END DO
            ELSE
              DO J=1,LHK1
                JHJ1 = IK(IMIN+J)
                IC = MIN(J,JHJ1-JHJ)
                SDOT = 0.0d0
                IF(IC.GT.0) THEN
                  DO I=0,IC-1
                    SDOT = SDOT+AL(JHJ1-IC+I)*AU(JHK+J-IC+I)
                  END DO
                  AU(JHK+J) = AU(JHK+J)-SDOT
                  SDOT = 0.0d0
                  DO I=0,IC-1
                    SDOT = SDOT+AU(JHJ1-IC+I)*AL(JHK+J-IC+I)
                  END DO
                END IF
                AL(JHK+J) = (AL(JHK+J)-SDOT)*AD(IMIN1+J)
                JHJ = JHJ1
              END DO
            END IF
!
          END IF
          IF(NSYM.EQ.0) THEN
            DO I=0,LHK1
              T1 = AU(JHK+I)
              T2 = T1*AD(IMIN1+I)
              AU(JHK+I) = T2
              SUMD = SUMD+T1*T2
            END DO
          ELSE
            DO I=0,LHK1
              SUMD = SUMD+AU(JHK+I)*AL(JHK+I)
            END DO
          END IF
        END IF
        If (AD(K)-SUMD .EQ. 0.d0) Then
          CALL ShowSevereError('AirflowNetworkSolver: L-U factorization in Subroutine FACSKY.')
          CALL ShowContinueError('The denominator used in L-U factorizationis equal to 0.0 at node = ' &
                                //TRIM(AirflowNetworkNodeData(K)%Name)//'.')
          CALL ShowContinueError('One possible cause is that this node may not be connected directly, or indirectly via airflow ' &
                                //'network connections ')
          CALL ShowContinueError('(e.g., AirflowNetwork:Multizone:SurfaceCrack, AirflowNetwork:Multizone:Component:' &
                                 // 'SimpleOpening, etc.), to an external')
          CALL ShowContinueError('node (AirflowNetwork:MultiZone:Surface).')
          CALL ShowContinueError('Please send your input file and weather file to EnergyPlus support/development team' &
                                 //' for further investigation.')
          CALL ShowFatalError('Preceding condition causes termination.')
        End If
        AD(K) = 1.0d0/(AD(K)-SUMD)
        JHK = JHK1
      END DO
!
      RETURN
      END SUBROUTINE FACSKY


      SUBROUTINE SLVSKY(AU,AD,AL,B,IK,NEQ,NSYM)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  This subroutine is revised from CLVSKY developed by George Walton, NIST

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine solves simultaneous linear algebraic equations [A] * X = B
          ! using L-U factored skyline form of [A] from "FACSKY"

          ! 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)  :: IK(NetworkNumOfNodes+1) ! pointer to the top of column/row "K"
          INTEGER, INTENT(IN)  :: NEQ   ! number of equations
          INTEGER, INTENT(IN)  :: NSYM  ! symmetry:  0 = symmetric matrix, 1 = non-symmetric
          ! noel, GNU says the AU is indexed above its upper bound
          !REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)-1) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(IN) :: AU(IK(NetworkNumOfNodes+1)) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(IN) :: AD(NetworkNumOfNodes)  ! the main diagonal of [A] before and after factoring
          REAL(r64), INTENT(IN) :: AL(IK(NetworkNumOfNodes+1)-1) ! the lower triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: B(NetworkNumOfNodes) ! "B" vector (input); "X" vector (output).


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:

      INTEGER  I, JHK, JHK1, K, LHK
      REAL(r64)  SDOT, T1

      ! FLOW:
      JHK = 1
      DO K=2,NEQ
        JHK1 = IK(K+1)
        LHK = JHK1-JHK
        IF(LHK.LE.0) CYCLE
        SDOT = 0.0d0
        IF(NSYM.EQ.0) THEN
          DO I=0,LHK-1
            SDOT = SDOT+AU(JHK+I)*B(K-LHK+I)
          END DO
        ELSE
          DO I=0,LHK-1
            SDOT = SDOT+AL(JHK+I)*B(K-LHK+I)
          END DO
        ENDIF
        B(K) = B(K)-SDOT
        JHK = JHK1
      END DO
!
      IF(NSYM.EQ.0) THEN
        DO K=1,NEQ
          B(K) = B(K)*AD(K)
        END DO
      END IF
!
      K = NEQ+1
      JHK1 = IK(K)
      DO WHILE (K .NE. 1)
        K = K-1
        IF(NSYM.EQ.1) B(K) = B(K)*AD(K)
        IF(K.EQ.1) EXIT
!        IF(K.EQ.1) RETURN
        JHK = IK(K)
        T1 = B(K)
        DO I=0,JHK1-JHK-1
          B(K-JHK1+JHK+I) = B(K-JHK1+JHK+I)-AU(JHK+I)*T1
        END DO
        JHK1 = JHK
      END DO

      RETURN
      END SUBROUTINE SLVSKY


      SUBROUTINE FILSKY(X,LM,IK,AU,AD,FLAG)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  This subroutine is revised from FILSKY developed by George Walton, NIST

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine adds element array "X" to the sparse skyline matrix [A]

          ! 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)  :: LM(2) ! location matrix
          INTEGER, INTENT(IN)  :: IK(NetworkNumOfNodes+1) ! pointer to the top of column/row "K"
          INTEGER, INTENT(IN)  :: FLAG  ! mode of operation
          REAL(r64), INTENT(IN)     :: X(4)  ! element array (row-wise sequence)
          ! noel, GNU says the AU is indexed above its upper bound
          !REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)-1) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AU(IK(NetworkNumOfNodes+1)) ! the upper triangle of [A] before and after factoring
          REAL(r64), INTENT(INOUT) :: AD(NetworkNumOfNodes)  ! the main diagonal of [A] before and after factoring


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
         INTEGER  J, K, L

      ! FLOW:
      ! K = row number, L = column number.
      IF(FLAG.GT.1) THEN
        K = LM(1)
        L = LM(2)
        IF(FLAG.EQ.4) THEN
          AD(K) = AD(K)+X(1)
          IF(K.LT.L) THEN
            J = IK(L+1)-L+K
            AU(J) = AU(J)+X(2)
          ELSE
            J = IK(K+1)-K+L
            AU(J) = AU(J)+X(3)
          END IF
          AD(L) = AD(L)+X(4)
        ELSE IF(FLAG.EQ.3) THEN
          AD(L) = AD(L)+X(4)
        ELSE IF(FLAG.EQ.2) THEN
          AD(K) = AD(K)+X(1)
        END IF
      END IF
!
      RETURN
      END SUBROUTINE FILSKY


      SUBROUTINE DUMPVD(S,V,N,UOUT)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  This subroutine is revised from DUMPVD developed by George Walton, NIST

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine prints the contents of the REAL(r64) "V" vector

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          CHARACTER(*), INTENT(IN):: S  ! Description
          REAL(r64), INTENT(IN)     :: V(*) ! Output values
          INTEGER, INTENT(IN)  :: N  ! Array size
          INTEGER, INTENT(IN)  :: UOUT ! Output file unit


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          INTEGER I

      ! FLOW:
      ! Write values for debug
      WRITE(UOUT,901) S
      WRITE(UOUT,902) (V(I),I=1,N)
  901 FORMAT(1X,A)
  902 FORMAT(1X,5E15.07)
!
      RETURN
      END SUBROUTINE DUMPVD


      SUBROUTINE DUMPVR(S,V,N,UOUT)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         George Walton
          !       DATE WRITTEN   Extracted from AIRNET
          !       MODIFIED       Lixing Gu, 2/1/04
          !                      Revised the subroutine to meet E+ needs
          !       MODIFIED       Lixing Gu, 6/8/05
          !
          !       RE-ENGINEERED  This subroutine is revised from DUMPVR developed by George Walton, NIST

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine prints the contents of the REAL(r64) "V" vector

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          CHARACTER(*), INTENT(IN):: S  ! Description
          REAL(r64), INTENT(IN)     :: V(*) ! Output values
          INTEGER, INTENT(IN)  :: N  ! Array size
          INTEGER, INTENT(IN)  :: UOUT ! Output file unit


          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
          INTEGER I

      ! FLOW:
      WRITE(UOUT,901) S
      WRITE(UOUT,902) (V(I),I=1,N)
  901 FORMAT(1X,A)
  902 FORMAT(1X,5E15.07)
!
      RETURN
      END SUBROUTINE DUMPVR


      SUBROUTINE AFEDOP(J,LFLAG,PDROP,IL,N,M,F,DF,NF)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Oct. 2005
          !       MODIFIED       na
          !       RE-ENGINEERED  This subroutine is revised based on a vertical large opening subroutine from COMIS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine simulates airflow and pressure of a detailed large opening component.

          ! METHODOLOGY EMPLOYED:
          ! Purpose:  This routine calculates the massflow and its derivative
          !       through a large opening in both flow directions. As input
          !       the density profiles RhoProfF/T are required aswell as the
          !       effective pressure difference profile DpProfNew, which is the
          !       sum of the stack pressure difference profile DpProf and the
          !       difference of the actual pressures at reference height. The
          !       profiles are calculated in the routine PresProfile.
          !       The massflow and its derivative are calculated for each
          !       interval representing a step of the pressure difference
          !       profile. The total flow and derivative are obtained by
          !       summation over the whole opening.
          !       The calculation is split into different cases representing
          !       different situations of the opening:
          !       - closed opening (opening factor = 0): summation of top and
          !         bottom crack (crack length = lwmax) plus "integration" over
          !         a vertically distributed crack of length (2*lhmax+lextra).
          !       - type 1: normal rectangular opening: "integration" over NrInt
          !         openings with width actlw and height actlh/NrInt
          !       - type 2: horizontally pivoted window: flow direction assumed
          !         strictly perpendicular to the plane of the opening
          !         -> "integration" over normal rectangular openings at top
          !         and bottom of LO plus a rectangular opening in series with two
          !         triangular openings in the middle of the LO (most general
          !         situation). The geometry is defined by the input parameters
          !         actlw(=lwmax), actlh, axisheight, opening angle.
          !       Assuming the massflow perpendicular to the opening plane in all
          !       cases the ownheightfactor has no influence on the massflow.

          ! REFERENCES:
          ! Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
          ! Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: J     ! Component number
          INTEGER, INTENT(IN)  :: LFLAG ! Initialization flag.If = 1, use laminar relationship
          REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
          INTEGER, INTENT(IN)  :: IL    ! Linkage number
          INTEGER, INTENT(IN)  :: N     ! Node 1 number
          INTEGER, INTENT(IN)  :: M     ! Node 2 number
          INTEGER, INTENT(OUT) :: NF    ! Number of flows, either 1 or 2
          REAL(r64), INTENT(OUT) :: F(2)  ! Airflow through the component [kg/s]
          REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative:  DF/DP

          ! SUBROUTINE PARAMETER DEFINITIONS:
          REAL(r64), PARAMETER :: RealMax=0.1d+37
          REAL(r64), PARAMETER :: RealMin=1d-37


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
      INTEGER  CompNum
      REAL(r64) Width,Height

      REAL(r64)  fma12    ! massflow in direction "from-to" [kg/s]
      REAL(r64)  fma21    ! massflow in direction "to-from" [kg/s]
      REAL(r64)  dp1fma12 ! derivative d fma12 / d Dp [kg/s/Pa]
      REAL(r64)  dp1fma21 ! derivative d fma21 / d Dp [kg/s/Pa]
      REAL(r64)  DpProfNew(NrInt+2) ! Differential pressure profile for Large Openings, taking into account fixed
                                           ! pressures and actual zone pressures at reference height
      REAL(r64) fact ! Actual opening factor
      REAL(r64) DifLim ! Limit for the pressure difference where laminarization takes place [Pa]
      REAL(r64) Cfact,FvVeloc

      REAL(r64) ActLh,ActLw,Lextra,Axishght,ActCD,Cs,expn,Type
      REAL(r64) Interval,fmasum,dfmasum,Prefact,EvalHghts(NrInt+2)
      REAL(r64) h2,h4,alpha,rholink,c1,c2
      REAL(r64) DpZeroOffset,area,WFact,HFact
      INTEGER i,Loc,iNum

      ! FLOW:
      ! Get component properties
      DifLim = 1.0d-4
      CompNum = AirflowNetworkCompData(J)%TypeNum
      Width = MultizoneSurfaceData(IL)%Width
      Height = MultizoneSurfaceData(IL)%Height
      fact = MultizoneSurfaceData(IL)%OpenFactor
      Loc = (AirflowNetworkLinkageData(IL)%DetOpenNum-1)*(NrInt+2)
      iNum = MultizoneCompDetOpeningData(CompNum)%NumFac
      ActCD=0.0d0

      if (iNum .eq. 2) then
        if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac2) then
          WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%WidthFac2-MultizoneCompDetOpeningData(CompNum)%WidthFac1)
          HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%HeightFac2-MultizoneCompDetOpeningData(CompNum)%HeightFac1)
          CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%DischCoeff2-MultizoneCompDetOpeningData(CompNum)%DischCoeff1)
        Else
          CALL ShowFatalError('Open Factor is above the maximum input range for opening factors in ' &
                              //'AirflowNetwork:MultiZone:Component:DetailedOpening = ' &
                                   //TRIM(MultizoneCompDetOpeningData(CompNum)%Name))
        End if
      end if

      if (iNum .eq. 3) then
        if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac2) then
          WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%WidthFac2-MultizoneCompDetOpeningData(CompNum)%WidthFac1)
          HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%HeightFac2-MultizoneCompDetOpeningData(CompNum)%HeightFac1)
          CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%DischCoeff2-MultizoneCompDetOpeningData(CompNum)%DischCoeff1)
        else if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac3) then
          WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac2+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
           (MultizoneCompDetOpeningData(CompNum)%WidthFac3-MultizoneCompDetOpeningData(CompNum)%WidthFac2)
          HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac2+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
           (MultizoneCompDetOpeningData(CompNum)%HeightFac3-MultizoneCompDetOpeningData(CompNum)%HeightFac2)
          CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff2+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
           (MultizoneCompDetOpeningData(CompNum)%DischCoeff3-MultizoneCompDetOpeningData(CompNum)%DischCoeff2)
        Else
          CALL ShowFatalError('Open Factor is above the maximum input range for opening factors in ' &
                              //'AirflowNetwork:MultiZone:Component:DetailedOpening = ' &
                                   //TRIM(MultizoneCompDetOpeningData(CompNum)%Name))
        end if
      end if

      if (iNum .eq. 4) then
        if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac2) then
          WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%WidthFac2-MultizoneCompDetOpeningData(CompNum)%WidthFac1)
          HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%HeightFac2-MultizoneCompDetOpeningData(CompNum)%HeightFac1)
          CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff1+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
           (MultizoneCompDetOpeningData(CompNum)%DischCoeff2-MultizoneCompDetOpeningData(CompNum)%DischCoeff1)
        else if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac3) then
          WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac2+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
           (MultizoneCompDetOpeningData(CompNum)%WidthFac3-MultizoneCompDetOpeningData(CompNum)%WidthFac2)
          HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac2+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
           (MultizoneCompDetOpeningData(CompNum)%HeightFac3-MultizoneCompDetOpeningData(CompNum)%HeightFac2)
          CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff2+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
           (MultizoneCompDetOpeningData(CompNum)%DischCoeff3-MultizoneCompDetOpeningData(CompNum)%DischCoeff2)
        else if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac4) then
          WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac3+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac3)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac4-MultizoneCompDetOpeningData(CompNum)%OpenFac3)* &
           (MultizoneCompDetOpeningData(CompNum)%WidthFac4-MultizoneCompDetOpeningData(CompNum)%WidthFac3)
          HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac3+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac3)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac4-MultizoneCompDetOpeningData(CompNum)%OpenFac3)* &
           (MultizoneCompDetOpeningData(CompNum)%HeightFac4-MultizoneCompDetOpeningData(CompNum)%HeightFac3)
          CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff3+ &
           (fact-MultizoneCompDetOpeningData(CompNum)%OpenFac3)/ &
           (MultizoneCompDetOpeningData(CompNum)%OpenFac4-MultizoneCompDetOpeningData(CompNum)%OpenFac3)* &
           (MultizoneCompDetOpeningData(CompNum)%DischCoeff4-MultizoneCompDetOpeningData(CompNum)%DischCoeff3)
        Else
          CALL ShowFatalError('Open Factor is above the maximum input range for opening factors in ' &
                              //'AirflowNetwork:MultiZone:Component:DetailedOpening = ' &
                                   //TRIM(MultizoneCompDetOpeningData(CompNum)%Name))
        end if
      end if

      ! calculate DpProfNew
      Do i=1,NrInt+2
        DpProfNew(i) = PDROP+DpProf(Loc+i)-DpL(1,iL)
      End Do

      ! Get opening data based on the opening factor
      if (fact .eq. 0) then
        ActLw=MultizoneSurfaceData(IL)%Width
        ActLh=MultizoneSurfaceData(IL)%Height
        Cfact = 0.0d0
      else
        ActLw=MultizoneSurfaceData(IL)%Width*Wfact
        ActLh=MultizoneSurfaceData(IL)%Height*HFact
        ActCD=CFact
      end if

      Cs=MultizoneCompDetOpeningData(CompNum)%FlowCoef
      Expn=MultizoneCompDetOpeningData(CompNum)%FlowExpo
      Type=MultizoneCompDetOpeningData(CompNum)%LVOType
      IF (Type.EQ.1) THEN
        Lextra=MultizoneCompDetOpeningData(CompNum)%LVOValue
        AxisHght=0.0d0
      ELSE IF (Type.EQ.2) THEN
        Lextra=0.0d0
        AxisHght=MultizoneCompDetOpeningData(CompNum)%LVOValue
        ActLw=MultizoneSurfaceData(IL)%Width
        ActLh=MultizoneSurfaceData(IL)%Height
      ENDIF

      ! Add window multiplier with window close
      If (MultizoneSurfaceData(IL)%Multiplier > 1.0d0) Cs=Cs*MultizoneSurfaceData(IL)%Multiplier
      ! Add window multiplier with window open
      If (fact .gt. 0.0d0) then
        If (MultizoneSurfaceData(IL)%Multiplier > 1.0d0) ActLw=ActLw*MultizoneSurfaceData(IL)%Multiplier
      End If

      ! Add recurring warnings
      If (fact .gt. 0.0d0) then
        If (ActLw .eq. 0.0d0) Then
          MultizoneCompDetOpeningData(CompNum)%WidthErrCount = MultizoneCompDetOpeningData(CompNum)%WidthErrCount + 1
          if (MultizoneCompDetOpeningData(CompNum)%WidthErrCount< 2) then
            CALL ShowWarningError('The actual width of the AirflowNetwork:MultiZone:Component:DetailedOpening of '//&
              TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0.')
            CALL ShowContinueError('The actual width is set to 1.0E-6 m.')
            CALL ShowContinueErrorTimeStamp(' Occurrence info: ')
          else
            CALL ShowRecurringWarningErrorAtEnd('The actual width of the '// &
              'AirflowNetwork:MultiZone:Component:DetailedOpening of ' &
              //TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0 error continues.', &
              MultizoneCompDetOpeningData(CompNum)%WidthErrIndex, ActLw,ActLw)
          end if
          ActLw = 1.0d-6
        End If
        If (ActLh .eq. 0.0d0) Then
          MultizoneCompDetOpeningData(CompNum)%HeightErrCount = MultizoneCompDetOpeningData(CompNum)%HeightErrCount + 1
          if (MultizoneCompDetOpeningData(CompNum)%HeightErrCount< 2) then
            CALL ShowWarningError('The actual height of the AirflowNetwork:MultiZone:Component:DetailedOpening of '//&
              TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0.')
            CALL ShowContinueError('The actual height is set to 1.0E-6 m.')
            CALL ShowContinueErrorTimeStamp(' Occurrence info: ')
          else
            CALL ShowRecurringWarningErrorAtEnd('The actual width of the '// &
              'AirflowNetwork:MultiZone:Component:DetailedOpening of ' &
              //TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0 error continues.', &
              MultizoneCompDetOpeningData(CompNum)%HeightErrIndex, ActLh,ActLh)
          end if
          ActLh = 1.0d-6
        End If
      End If
      ! Initialization:
      NF = 1
      Interval=ActLh/NrInt
      fma12=0.0d0
      fma21=0.0d0
      dp1fma12=0.0d0
      dp1fma21=0.0d0

      ! Closed LO
      IF (Cfact.EQ.0) THEN
        DpZeroOffset=DifLim
        ! bottom crack
        IF (DpProfNew(1).GT.0) THEN
          IF (Abs(DpProfNew(1)).LE.DpZeroOffset) THEN
            dfmasum=Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
            fmasum=DpProfNew(1)*dfmasum
          ELSE
            fmasum=Cs*ActLw*(DpProfNew(1))**expn
            dfmasum=fmasum*expn/DpProfNew(1)
          ENDIF
          fma12=fma12+fmasum
          dp1fma12=dp1fma12+dfmasum
        ELSE
          IF (Abs(DpProfNew(1)).LE.DpZeroOffset) THEN
            dfmasum=-Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
            fmasum=DpProfNew(1)*dfmasum
          ELSE
            fmasum=Cs*ActLw*(-DpProfNew(1))**expn
            dfmasum=fmasum*expn/DpProfNew(1)
          ENDIF
          fma21=fma21+fmasum
          dp1fma21=dp1fma21+dfmasum
        ENDIF
        ! top crack
        IF (DpProfNew(NrInt+2).GT.0) THEN
          IF (Abs(DpProfNew(NrInt+2)).LE.DpZeroOffset) THEN
            dfmasum=Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
            fmasum=DpProfNew(NrInt+2)*dfmasum
          ELSE
            fmasum=Cs*ActLw*(DpProfNew(NrInt+2))**expn
            dfmasum=fmasum*expn/DpProfNew(NrInt+2)
          ENDIF
          fma12=fma12+fmasum
          dp1fma12=dp1fma12+dfmasum
        ELSE
          IF (Abs(DpProfNew(NrInt+2)).LE.DpZeroOffset) THEN
            dfmasum=-Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
            fmasum=DpProfNew(NrInt+2)*dfmasum
          ELSE
            fmasum=Cs*ActLw*(-DpProfNew(NrInt+2))**expn
            dfmasum=fmasum*expn/DpProfNew(NrInt+2)
          ENDIF
          fma21=fma21+fmasum
          dp1fma21=dp1fma21+dfmasum
        ENDIF
        ! side and extra cracks
        Prefact=Interval*(2+lextra/ActLh)*Cs
        DO i=2,NrInt+1
          IF (DpProfNew(i).GT.0) THEN
            IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
              dfmasum=Prefact*DpZeroOffset**expn/DpZeroOffset
              fmasum=DpProfNew(i)*dfmasum
            ELSE
              fmasum=Prefact*(DpProfNew(i))**expn
              dfmasum=fmasum*expn/DpProfNew(i)
            ENDIF
            fma12=fma12+fmasum
            dp1fma12=dp1fma12+dfmasum
          ELSE
            IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
              dfmasum=-Prefact*DpZeroOffset**expn/DpZeroOffset
              fmasum=DpProfNew(i)*dfmasum
            ELSE
              fmasum=Prefact*(-DpProfNew(i))**expn
              dfmasum=fmasum*expn/DpProfNew(i)
            ENDIF
            fma21=fma21+fmasum
            dp1fma21=dp1fma21+dfmasum
          ENDIF
        ENDDO
      ENDIF

      ! Open LO, type 1
      IF ((Cfact.NE.0).AND.(type.EQ.1))  THEN
        DpZeroOffset=DifLim*1d-3
        Prefact=ActLw*ActCd*Interval*SQRT(2.0d0)
        DO i=2,NrInt+1
          IF (DpProfNew(i).GT.0) THEN
            IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
              dfmasum=SQRT(RhoProfF(Loc+i)*DpZeroOffset)/DpZeroOffset
              fmasum=DpProfNew(i)*dfmasum
            ELSE
              fmasum=SQRT(RhoProfF(Loc+i)*DpProfNew(i))
              dfmasum=0.5d0*fmasum/DpProfNew(i)
            ENDIF
            fma12=fma12+fmasum
            dp1fma12=dp1fma12+dfmasum
          ELSE
            IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
              dfmasum=-SQRT(RhoProfT(Loc+i)*DpZeroOffset)/DpZeroOffset
              fmasum=DpProfNew(i)*dfmasum
            ELSE
              fmasum=SQRT(-RhoProfT(Loc+i)*DpProfNew(i))
              dfmasum=0.5d0*fmasum/DpProfNew(i)
            ENDIF
            fma21=fma21+fmasum
            dp1fma21=dp1fma21+dfmasum
          ENDIF
        ENDDO

        fma12=prefact*fma12
        fma21=prefact*fma21
        dp1fma12=prefact*dp1fma12
        dp1fma21=prefact*dp1fma21

      ENDIF

      ! Open LO, type 2
      IF ((Cfact.NE.0).AND.(type.EQ.2))  THEN
        ! Initialization
        DpZeroOffset=DifLim*1d-3
        ! New definition for opening factors for LVO type 2: opening angle = 90 degrees --> opening factor = 1.0
        ! should be PIOvr2 in below?
        alpha=fact*3.14159d0/2.d0
        h2=Axishght*(1.0d0-COS(alpha))
        h4=Axishght+(ActLh-Axishght)*COS(alpha)
        EvalHghts(1)=0.0d0
        EvalHghts(NrInt+2)=ActLh
        ! New definition for opening factors for LVO type 2: pening angle = 90 degrees --> opening factor = 1.0
        IF (fact.EQ.1.0d0) THEN
          h2=Axishght
          h4=Axishght
        ENDIF

        DO i=2,NrInt+1
          EvalHghts(i)=Interval*(i-1.5d0)
        ENDDO

        ! Calculation of massflow and its derivative
        DO i=2,NrInt+1
          IF (DpProfNew(i).GT.0) THEN
            rholink=RhoProfF(Loc+i)
          ELSE
            rholink=RhoProfT(Loc+i)
          ENDIF

          IF ((EvalHghts(i).LE.h2).OR.(EvalHghts(i).GE.h4)) THEN
            IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
              dfmasum=ActCd*ActLw*Interval*SQRT(2.0*rholink* &
              DpZeroOffset)/DpZeroOffset*DSIGN(1d0,DpProfNew(i))
              fmasum=DpProfNew(i)*dfmasum
            ELSE
              fmasum=ActCd*ActLw*Interval*SQRT(2.0*rholink*ABS(DpProfNew(i)))
              dfmasum=0.5d0*fmasum/DpProfNew(i)
            ENDIF
          ELSE
            ! triangular opening at the side of LO
            c1=ActCd*ActLw*Interval*SQRT(2.0*rholink)
            c2=2*ActCd*ABS((Axishght-Evalhghts(i)))*TAN(alpha)* &
              Interval*SQRT(2.0d0*rholink)
            IF ((c1.NE.0).AND.(c2.NE.0)) THEN
              IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
                dfmasum=SQRT(DpZeroOffset/(1/c1/c1+1/c2/c2))/ &
                        DpZeroOffset*DSIGN(1d0,DpProfNew(i))
                fmasum=DpProfNew(i)*dfmasum
              ELSE
                fmasum=SQRT(ABS(DpProfNew(i))/(1/c1/c1+1/c2/c2))
                dfmasum=0.5d0*fmasum/DpProfNew(i)
              ENDIF
            ELSE
              fmasum=0.0d0
              dfmasum=0.0d0
            ENDIF
          ENDIF

          IF (DpProfNew(i).GT.0) THEN
            fma12=fma12+fmasum
            dp1fma12=dp1fma12+dfmasum
          ELSE
            fma21=fma21+fmasum
            dp1fma21=dp1fma21+dfmasum
          ENDIF
        ENDDO

      ENDIF

      ! Calculate some velocity in the large opening
      area=ActLh*ActLw*actCD
      if (area.gt.(cs+RealMin)) then
        if (area.gt.RealMin) then
          FvVeloc=(fma21+fma12)/area
        else
          FvVeloc=0.0d0
        end if
      else
        ! here the average velocity over the full area, may blow half in half out.
        ! velocity= Fva/Nett area=Fma/Rho/(Cm/( (2**N)* sqrt(1.2) ) )
        if (cs.gt.0.0d0) then
          ! get the average Rho for this closed window
          DO i=2,NrInt+1
            rholink=0.0d0
            IF (DpProfNew(i).GT.0) THEN
              rholink=RhoProfF(Loc+i)
            ELSE
              rholink=RhoProfT(Loc+i)
            ENDIF
            rholink=rholink/nrInt
            rholink=1.2d0
          END DO
          FvVeloc=(Fma21+fma12)*(2.d0**ExpN)*SQRT(1.2d0)/(rhoLink*Cs)
        else
          FvVeloc=0.0d0
        end if
      end if

      ! Output mass flow rates and associated derivatives
      F(1) = fma12-fma21
      DF(1) = dp1fma12-dp1fma21
      F(2) = 0.0d0
      if (fma12 .NE. 0.0d0 .and. fma21 .NE. 0.0d0) then
        F(2) = fma21
      End if
      DF(2) = 0.0d0

      RETURN
    END SUBROUTINE AFEDOP

    SUBROUTINE PresProfile(il,Pprof,G,DpF,DpT,BetaF,BetaT,RhoStF, &
            RhoStT,From,To,ActLh,OwnHeightFactor)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Oct. 2005
          !       MODIFIED       na
          !       RE-ENGINEERED  This subroutine is revised based on PresProfile subroutine from COMIS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine calculates for a large opening profiles of stack pressure difference and
          ! densities in the zones linked by the a detailed opening cmponent.

          ! METHODOLOGY EMPLOYED:
          ! The profiles are obtained in the following
          ! way:    - the opening is divided into NrInt vertical intervals
          !         - the stack pressure difference and densities in From-
          !           and To-zone are calculated at the centre of each
          !           interval aswell as at the top and bottom of the LO
          !          - these values are stored in the (NrInt+2)-dimensional
          !             arrays DpProf, RhoProfF, RhoProfT.
          ! The calculation of stack pressure and density in the two zones
          ! is based on the arrays DpF/T, RhoStF/T, BetaF/T. These arrays
          ! are calculated in the COMIS routine Lclimb. They contain the
          ! values of stack pressure and density at the startheight of the
          ! opening and at startheights of all layers lying inside the
          ! opening, and the density gradients across the layers.
          ! The effective startheight zl(1/2) in the From/To zone and the
          ! effective length actLh of the LO take into account the
          ! startheightfactor, heightfactor and ownheightfactor. Thus for
          ! slanted windows the range of the profiles is the vertical
          ! projection of the actual opening.

          ! REFERENCES:
          ! Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
          ! Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          INTEGER, INTENT(IN)  :: IL       ! Linkage number
          INTEGER, INTENT(IN)  :: Pprof    ! Opening number
          REAL(r64), INTENT(IN)     :: G        ! gravitation field strength [N/kg]
          REAL(r64), INTENT(IN) :: DpF(2) ! Stack pressures at start heights of Layers
                                                 ! in the FROM zone (starting at linkheight) [Pa]
          REAL(r64), INTENT(IN) :: DpT(2) ! Stack pressures at start heights of Layers
          REAL(r64), INTENT(IN)     :: BetaF(2) ! Density gradients in the FROM zone (starting at linkheight) [Kg/m3/m]
          REAL(r64), INTENT(IN)     :: BetaT(2) ! Density gradients in the TO zone (starting at linkheight) [Kg/m3/m]
          REAL(r64), INTENT(IN)     :: RhoStF(2) ! Density at the start heights of Layers in the FROM zone
                                            ! (starting at linkheight) [Kg/m3]
          REAL(r64), INTENT(IN)     :: RhoStT(2) ! Density at the start heights of Layers in the TO zone
                                            ! (starting at linkheight) [Kg/m3]
          REAL(r64), INTENT(IN)     :: ActLh    ! Actual height of opening [m]
          REAL(r64), INTENT(IN)     :: OwnHeightFactor ! Cosine of deviation angle of the opening plane from the vertical direction
          INTEGER, INTENT(IN)  :: From     ! Number of FROM zone
          INTEGER, INTENT(IN)  :: To       ! Number of To zone

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
      REAL(r64) zF(2),zT(2) ! Startheights of layers in FROM-, TO-zone
      REAL(r64) zStF(2),zStT(2) ! Startheights of layers within the LO, starting with the actual startheight of the LO.
                           ! The values in the arrays DpF, DpT, BetaF, BetaT, RhoStF, RhoStT are calculated at these heights.
      REAL(r64) hghtsFR, hghtsTR
      REAL(r64) hghtsF(NrInt+2),hghtsT(NrInt+2) ! Heights of evaluation points for pressure and density profiles
      REAL(r64) Interval ! Distance between two evaluation points
      REAL(r64) delzF,delzT ! Interval between actual evaluation point and startheight of actual layer in FROM-, TO-zone
      INTEGER AnzLayF,AnzLayT ! Number of layers in FROM-, TO-zone
      INTEGER lF,lT ! Actual index for DpF/T, BetaF/T, RhoStF/T, zStF/T
      INTEGER n,i,k


      ! FLOW:
      ! Initialization
      delzF=0.0d0
      delzT=0.0d0
      Interval=ActLh*OwnHeightFactor/NrInt

      DO n=1,NrInt
        hghtsF(n+1)=AirflowNetworkLinkageData(il)%NodeHeights(1)+Interval*(n-0.5d0)
        hghtsT(n+1)=AirflowNetworkLinkageData(il)%NodeHeights(2)+Interval*(n-0.5d0)
      END DO
      hghtsF(1)=AirflowNetworkLinkageData(il)%NodeHeights(1)
      hghtsT(1)=AirflowNetworkLinkageData(il)%NodeHeights(2)
      hghtsF(NrInt+2)=AirflowNetworkLinkageData(il)%NodeHeights(1)+ActLh*OwnHeightFactor
      hghtsT(NrInt+2)=AirflowNetworkLinkageData(il)%NodeHeights(2)+ActLh*OwnHeightFactor

      lF=1
      lT=1
      If (From .eq. 0) then
        AnzLayF=1
      Else
        AnzLayF=0
      End if
      If (To .eq. 0) then
        AnzLayT=1
      Else
        AnzLayT=0
      End if

      IF (AnzLayF.GT.0) THEN
        DO n=1,AnzLayF
          zF(n)=0.0d0
          If (hghtsF(1) .LT. 0.0d0) zF(n)=hghtsF(1)
        ENDDO
      ENDIF

      IF (AnzLayT.GT.0) THEN
        DO n=1,AnzLayT
          zT(n)=0.0d0
          If (hghtsT(1) .LT. 0.0d0) zT(n)=hghtsT(1)
        ENDDO
      ENDIF

      zStF(1)=AirflowNetworkLinkageData(il)%NodeHeights(1)
      i=2
      k=1

      DO WHILE (k.LE.AnzLayF)
        IF (zF(k).GT.zStF(1)) Exit
        k=k+1
      END DO

      DO WHILE (k.LE.AnzLayF)
        IF (zF(k).GT.hghtsF(NrInt)) EXIT
        zStF(i)=zF(k) !Objexx:BoundsViolation zStF(i) @ i>2 and zF(k) @ k>2
        i=i+1
        k=k+1
      END DO

      zStF(i)=AirflowNetworkLinkageData(il)%NodeHeights(1)+ActLh*OwnHeightFactor !Objexx:BoundsViolation zStF(i) @ i>2
      zStT(1)=AirflowNetworkLinkageData(il)%NodeHeights(2)
      i=2
      k=1

      DO WHILE (k.LE.AnzLayT)
        IF (zT(k).GT.zStT(1)) EXIT
        k=k+1
      END DO

      DO WHILE (k.LE.AnzLayT)
        IF (zT(k).GT.hghtsT(NrInt)) EXIT !Objexx:BoundsViolation zT(k) @ k>2
        zStT(i)=zT(k) !Objexx:BoundsViolation zStF(i) @ i>2 and zT(k) @ k>2
        i=i+1
        k=k+1
      END DO

      zStT(i)=AirflowNetworkLinkageData(il)%NodeHeights(2)+ActLh*OwnHeightFactor !Objexx:BoundsViolation zStT(i) @ i>2

      ! Calculation of DpProf, RhoProfF, RhoProfT
      DO i=1,NrInt+2
        hghtsFR=hghtsF(i)
        hghtsTR=hghtsT(i)

        DO
          IF (hghtsFR.GT.zStF(lF+1)) THEN
            IF (lF .gt. 2) EXIT
            lF=lF+1
          END IF
          IF (hghtsFR.LE.zStF(lF+1)) Exit
        END DO

        DO
          IF (hghtsTR.GT.zStT(lT+1)) THEN
            lT=lT+1
          END IF
          IF (hghtsTR.LE.zStT(lT+1)) Exit
        END DO

        delzF=hghtsF(i)-zStF(lF)
        delzT=hghtsT(i)-zStT(lT)

        RhoProfF(i+Pprof)=RhoStF(lF)+BetaF(lF)*delzF
        RhoProfT(i+Pprof)=RhoStT(lT)+BetaT(lT)*delzT

        DpProf(i+Pprof)=DpF(lF)-DpT(lT) &
            -G*(RhoStF(lF)*delzF+BetaF(lF)*delzF**2/2.0d0) &
            +G*(RhoStT(lT)*delzT+BetaT(lT)*delzT**2/2.0d0)
      END DO

      RETURN
    END SUBROUTINE PresProfile

    SUBROUTINE PStack

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Oct. 2005
          !       MODIFIED       na
          !       RE-ENGINEERED  This subroutine is revised based on PresProfile subroutine from COMIS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine calculates the stack pressures for a link between two zones

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
          ! Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          ! na

          ! SUBROUTINE PARAMETER DEFINITIONS:
          REAL(r64), PARAMETER :: PI = 3.14159265358979323846d0
          REAL(r64), PARAMETER :: PSea=101325.0d0

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
!      REAL(r64) RhoOut ! air density outside [kg/m3]
      REAL(r64) G ! gravity field strength [N/kg]
      REAL(r64) RhoL1,RhoL2 ! Air density [kg/m3]
      REAL(r64) Pbz ! Pbarom at entrance level [Pa]
      REAL(r64) RhoDrL(2,NumOfLinksMultiZone) ! dry air density on both sides of the link [kg/m3]
      REAL(r64) TempL1,TempL2 ! Temp in From and To zone at link level [C]
!      REAL(r64) Tout ! outside temperature [C]
      REAL(r64) Xhl1,Xhl2 ! Humidity in From and To zone at link level [kg/kg]
!      REAL(r64) Xhout ! outside humidity [kg/kg]
      REAL(r64) Hfl(NumOfLinksMultiZone) ! Own height factor for large (slanted) openings
      INTEGER Nl ! number of links

      REAL(r64) DpF(2),DpP,DpT(2)
      REAL(r64) H
      REAL(r64) RhoStF(2),RhoStT(2),RhoDrDummi
      REAL(r64) BetaStF(2),BetaStT(2)
      REAL(r64) T,X,HSt(2)
      REAL(r64) TzFrom,XhzFrom
      REAL(r64) TzTo, XhzTo
      REAL(r64) ActLh,ActLOwnh
      REAL(r64) Pref
      REAL(r64) PzFrom,PzTo
      REAL(r64) RhoLd(2),RhoStd
      INTEGER From,To,Fromz,Toz,Ltyp
      INTEGER I,ll,j,k,Pprof
      INTEGER ilayptr,OpenNum

      REAL(r64) RhoREF, CONV

      ! FLOW:
      RhoREF = PsyRhoAirFnPbTdbW(Psea,OutDryBulbTemp,OutHumRat)

      CONV=Latitude*2.d0*PI/360.d0
      G=9.780373d0*(1.d0+.0052891d0*(SIN(CONV))**2-.0000059d0*(SIN(2.d0*CONV))**2)

      Hfl = 1.0d0
      Pbz = OutBaroPress
      Nl = NumOfLinksMultiZone
      OpenNum=0
      RhoLd(1)=1.2d0
      RhoLd(2)=1.2d0
      RhoStd=1.2d0

      DO I=1,Nl
        ! Check surface tilt
        If (AirflowNetworkLinkageData(i)%DetOpenNum > 0 .and. &
            Surface(MultizoneSurfaceData(I)%SurfNum)%Tilt .LT. 90) then
          Hfl(i) = Surface(MultizoneSurfaceData(I)%SurfNum)%SinTilt
        End If
        ! Initialisation
        From=AirflowNetworkLinkageData(i)%NodeNums(1)
        To=AirflowNetworkLinkageData(i)%NodeNums(2)
        if (AirflowNetworkNodeData(From)%EPlusZoneNum > 0 .and.AirflowNetworkNodeData(To)%EPlusZoneNum > 0) then
          LL=0
        else if(AirflowNetworkNodeData(From)%EPlusZoneNum == 0 .and.AirflowNetworkNodeData(To)%EPlusZoneNum > 0) then
          LL=1
        else
          LL=3
        End if

        Ltyp=AirflowNetworkCompData(AirflowNetworkLinkageData(i)%CompNum)%CompTypeNum
        IF (Ltyp.EQ.CompTypeNum_DOP) THEN
          ActLh=MultizoneSurfaceData(i)%Height
          ActLOwnh=ActLh*1.0d0
        ELSE
          ActLh=0.0d0
          ActLOwnh=0.0d0
        ENDIF

        TempL1=Tz(From)
        XhL1=Wz(From)
        TzFrom=Tz(From)
        XhzFrom=Wz(From)
        RhoL1=RhoZ(From)
        IF (LL.eq.0 .OR. LL.EQ.3) THEN
          PzFrom=Pz(From)
        ELSE
          PzFrom=0.0d0
          From=0
        ENDIF

        ilayptr = 0
        if (from .eq. 0) ilayptr = 1
        IF (ilayptr.EQ.0) THEN
          Fromz=0
        ELSE
          Fromz=From
        ENDIF

        TempL2=Tz(To)
        XhL2=Wz(To)
        TzTo=Tz(To)
        XhzTo=Wz(To)
        RhoL2=RhoZ(To)

        IF (LL.LT.3) THEN
          PzTo=Pz(To)
        ELSE
          PzTo=0.0d0
          To=0
        ENDIF
        ilayptr = 0
        if (To .eq. 0) ilayptr = 1
        IF (ilayptr.EQ.0) THEN
          Toz=0
        ELSE
          Toz=To
        ENDIF

        ! RhoDrL is Rho at link level without pollutant but with humidity
        RhoDrL(1,i) = PsyRhoAirFnPbTdbW(OutBaroPress+PzFrom,TempL1,XhL1)
        RhoDrL(2,i) = PsyRhoAirFnPbTdbW(OutBaroPress+PzTo,TempL2,XhL2)

        ! End initialisation

        ! calculate DpF the difference between Pz and P at Node 1 height
        ilayptr = 0
        If (fromz .eq. 0) ilayptr = 1
        J=ilayptr
        k=1
        CALL LClimb(G,Rhold(1),AirflowNetworkLinkageData(i)%NodeHeights(1),TempL1,XhL1, &
                   DpF(k),Toz,PzTo,Pbz,RhoDrL(1,I))
        RhoL1=Rhold(1)
        ! For large openings calculate the stack pressure difference profile and the
        ! density profile within the the top- and the bottom- height of the large opening
        IF (ActLOwnh.GT.0.0d0) THEN
          HSt(k)=AirflowNetworkLinkageData(i)%NodeHeights(1)
          RhoStF(k)=RhoL1
          k=k+1
          Hst(k)=0.0d0
          If (Hst(k-1) .LT. 0.0d0) Hst(k) = Hst(k-1)

          ! Search for the first startheight of a layer which is within the top- and the
          ! bottom- height of the large opening.
          Do
            ilayptr = 0
            If (Fromz .eq. 0) ilayptr = 9
            IF ((J.GT.ilayptr).OR.(Hst(k).GT.AirflowNetworkLinkageData(i)%NodeHeights(1))) Exit
            J=J+9
            Hst(k)=0.0d0
            If (Hst(k-1) .LT. 0.0d0) Hst(k) = Hst(k-1)
          End Do

          ! Calculate Rho and stack pressure for every StartHeight of a layer which is
          ! within the top- and the bottom-height of the  large opening.
          Do
            ilayptr = 0
            If (Fromz .eq. 0) ilayptr = 9
            IF ((J.GT.ilayptr).OR. &
               (Hst(k).GE.(AirflowNetworkLinkageData(i)%NodeHeights(1)+ActLOwnh))) Exit !Objexx:BoundsViolation HSt(k) @ k>2
            T=TzFrom
            X=XhzFrom
            CALL LClimb(G,RhoStd,HSt(k),T,X,DpF(k),Fromz,PzFrom,Pbz,RhoDrDummi) !Objexx:BoundsViolation HSt(k) and DpF(k) @ k>2
            RhoStF(k)=Rhostd !Objexx:BoundsViolation RhoStF(k) @ k>2
            J=J+9
            k=k+1 !Objexx k>2 now
            Hst(k)=0.0d0 !Objexx:BoundsViolation @ k>2
            If (Hst(k-1) .LT. 0.0d0) Hst(k) = Hst(k-1) !Objexx:BoundsViolation @ k>2
          End Do
          ! Stack pressure difference and rho for top-height of the large opening
          HSt(k)=AirflowNetworkLinkageData(i)%NodeHeights(1)+ActLOwnh !Objexx:BoundsViolation k>2 poss
          T=TzFrom
          X=XhzFrom
          CALL LClimb(G,RhoStd,HSt(k),T,X,DpF(k),Fromz,PzFrom,Pbz,RhoDrDummi) !Objexx:BoundsViolation k>2 poss
          RhoStF(k)=RhoStd !Objexx:BoundsViolation k >= 3 poss

          DO J=1,(k-1)
            BetaStF(J)=(RhoStF(J+1)-RhoStF(J))/(HSt(J+1)-HSt(J))
          END DO
        ENDIF

        ! repeat procedure for the "To" node, DpT
        ilayptr = 0
        if (Toz .eq. 0) ilayptr = 1
        J = ilayptr
        ! Calculate Rho at link height only if we have large openings or layered zones.
        k=1
        CALL LClimb(G,RhoLd(2),AirflowNetworkLinkageData(i)%NodeHeights(2),TempL2,XhL2, &
                   DpT(k),Toz,PzTo,Pbz,RhoDrL(2,I))
        RhoL2=RhoLd(2)

        ! For large openings calculate the stack pressure difference profile and the
        ! density profile within the the top- and the bottom- height of the large opening
        IF (ActLOwnh.GT.0.0d0) THEN
          HSt(k)=AirflowNetworkLinkageData(i)%NodeHeights(2)
          RhoStT(k)=RhoL2
          k=k+1
          Hst(k)=0.0d0
          If (Hst(k-1) .LT. 0.0d0) Hst(k) = Hst(k-1)
          Do
            ilayptr = 0
            if (Toz .eq. 0) ilayptr = 9
            IF ((J.GT.ilayptr).OR.(Hst(k).GT.AirflowNetworkLinkageData(i)%NodeHeights(2))) Exit
            J=J+9
            Hst(k)=0.0d0
            If (Hst(k-1) .LT. 0.0d0) Hst(k) = Hst(k-1)
          End Do
          ! Calculate Rho and stack pressure for every StartHeight of a layer which is
          ! within the top- and the bottom-height of the  large opening.
          Do
            ilayptr = 0
            if (Toz .eq. 0) ilayptr = 9
            IF ((J.GT.ilayptr).OR. (Hst(k).GE.(AirflowNetworkLinkageData(i)%NodeHeights(2)+ActLOwnh))) Exit !Objexx:BoundsViolation Hst(k) @ k>2
            T=TzTo
            X=XhzTo
            CALL LClimb(G,RhoStd,HSt(k),T,X,DpT(k),Toz,PzTo,Pbz,RhoDrDummi) !Objexx:BoundsViolation HSt(k) and DpT(k) @ k>2
            RhoStT(k)=RhoStd !Objexx:BoundsViolation RhoStT(k) @ k>2
            J=J+9
            k=k+1 !Objexx k>2 now
            Hst(k)=0.0d0 !Objexx:BoundsViolation @ k>2
            If (Hst(k-1) .LT. 0.0d0) Hst(k) = Hst(k-1) !Objexx:BoundsViolation @ k>2
          End Do
          ! Stack pressure difference and rho for top-height of the large opening
          HSt(k)=AirflowNetworkLinkageData(i)%NodeHeights(2)+ActLOwnh !Objexx:BoundsViolation k>2 poss
          T=TzTo
          X=XhzTo
          CALL LClimb(G,RhoStd,HSt(k),T,X,DpT(k),Toz,PzTo,Pbz,RhoDrDummi) !Objexx:BoundsViolation k>2 poss
          RhoStT(k)=RhoStd !Objexx:BoundsViolation k>2 poss

          DO J=1,(k-1)
            BetaStT(J)=(RhoStT(J+1)-RhoStT(J))/(HSt(J+1)-HSt(J))
          END DO
        ENDIF

        ! CALCULATE STACK PRESSURE FOR THE PATH ITSELF for different flow directions
        H=REAL(AirflowNetworkLinkageData(i)%NodeHeights(2),r64) -   &
                REAL(AirflowNetworkLinkageData(i)%NodeHeights(1),r64)
        IF (LL.eq.0 .OR. LL.EQ.3 .OR. LL.EQ.6) THEN
          H=H-AirflowNetworkNodeData(From)%NodeHeight
        ENDIF
        IF (LL.LT.3) THEN
          H=H+AirflowNetworkNodeData(To)%NodeHeight
        ENDIF

        ! IF AIR FLOWS from "From" to "To"
        Pref=Pbz+Pzfrom+DpF(1)
        DpP=psz(Pref,RhoLd(1),0.0d0,0.0d0,H,G)
        DpL(1,I)=(DpF(1)-DpT(1)+DpP)

        ! IF AIR FLOWS from "To" to "From"
        Pref=Pbz+Pzto+DpT(1)
        DpP=-psz(Pref,RhoLd(2),0.0d0,0.0d0,-H,G)
        DpL(2,I)=(DpF(1)-DpT(1)+DpP)

        IF (Ltyp.EQ.CompTypeNum_DOP) THEN
          Pprof=OpenNum*(NrInt+2)
          CALL PresProfile(i,Pprof,G,DpF,DpT,BetaStF,BetaStT,RhoStF,RhoStT,from,To,ActLh,Hfl(I))
          OpenNum = OpenNum+1
        ENDIF

      END DO

      RETURN

    END SUBROUTINE PStack

     REAL(r64) FUNCTION psz(Pz0,Rho0,beta,z0,z,g)

          ! FUNCTION INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Oct. 2005
          !       MODIFIED       na
          !       RE-ENGINEERED  This subroutine is revised based on psz function from COMIS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This function determines the pressure due to buoyancy in a stratified density environment

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          REAL(r64), INTENT(IN)     :: Pz0   ! Pressure at altitude z0 [Pa]
          REAL(r64), INTENT(IN)     :: Rho0  ! density at altitude z0 [kg/m3]
          REAL(r64), INTENT(IN)     :: beta  ! density gradient [kg/m4]
          REAL(r64), INTENT(IN)     :: z0    ! reference altitude [m]
          REAL(r64), INTENT(IN)     :: z     ! altitude[m]
          REAL(r64), INTENT(IN)     :: g     ! gravity field strength [N/kg]

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
      REAL(r64) dz,rho

      ! FLOW:
      dz=z-z0
      rho=(Rho0+beta*dz/2.0d0)
      Psz=-Pz0*(1.d0-EXP(-dz*rho*g/Pz0)) ! Differential pressure from z to z0 [Pa]

      RETURN
      END FUNCTION psz

      SUBROUTINE LClimb(G,Rho,Z,T,X,Dp,zone,Pz,Pbz,RhoDr)

          ! FUNCTION INFORMATION:
          !       AUTHOR         Lixing Gu
          !       DATE WRITTEN   Oct. 2005
          !       MODIFIED       na
          !       RE-ENGINEERED  This subroutine is revised based on subroutine IClimb from COMIS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This function the differential pressure from the reflevel in a zone To Z, the level of a link

          ! METHODOLOGY EMPLOYED:
          ! na

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
          ! na

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
          REAL(r64), INTENT(IN)                 :: g     ! gravity field strength [N/kg]
          REAL(r64), INTENT(INOUT)  :: Rho   ! Density link level (initialized with rho zone) [kg/m3]
          REAL(r64), INTENT(IN)                 :: Z     ! Height of the link above the zone reference [m]
          REAL(r64), INTENT(INOUT)              :: T     ! temperature at link level [C]
          REAL(r64), INTENT(INOUT)              :: X     ! absolute humidity at link level [kg/kg]
          REAL(r64), INTENT(OUT)    :: Dp    ! Stackpressure to the linklevel [Pa]
          INTEGER, INTENT(IN)              :: zone  ! Zone number
          REAL(r64), INTENT(IN)     :: Pz    ! Zone Pressure (reflevel) [Pa]
          REAL(r64), INTENT(IN)                 :: Pbz   ! Barometric pressure at entrance level [Pa]
          REAL(r64), INTENT(OUT)                :: RhoDr ! Air density of dry air on the link level used
                                                    ! for the concentration routine [kg/m3]

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na


          ! INTERFACE BLOCK SPECIFICATIONS
          ! na


          ! DERIVED TYPE DEFINITIONS
          ! na


          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
      REAL(r64) H         ! Start Height of the layer
      REAL(r64) BetaT     ! Temperature gradient of this layer
      REAL(r64) BetaXfct  ! Humidity gradient factor of this layer
      REAL(r64) BetaCfct  ! Concentration 1 gradient factor of this layer
      REAL(r64) X0
      REAL(r64) P,Htop,Hbot
      REAL(r64) Rho0,Rho1,Betarho
      INTEGER :: L=0
      INTEGER :: ilayptr=0

      ! FLOW:
      Dp=0.0d0
      Rho0=Rho
      x0=x
      IF (Z.GT.0.0d0) THEN
        ! initialize start values
        H=0.0d0
        BetaT=0.0d0
        BetaXfct=0.0d0
        BetaCfct=0.0d0
        BetaRho=0.0d0
        Hbot=0.0d0

        DO WHILE (H.LT.0.0d0)
          ! loop until H>0 ; The start of the layer is above 0
          BetaT=0.0d0
          BetaXfct=0.0d0
          BetaCfct=0.0d0
          L=L+9
          ilayptr = 0
          if (zone .eq. 0) ilayptr = 9
          IF (L.GE.ilayptr) THEN
            H=Z+1.0d0
          ELSE
            H=0.0d0
          ENDIF
        END DO

        ! The link is in this layer also if it is on the top of it.

        DO
          IF (H.GE.Z) THEN
            ! the link ends in this layer , we reached the final Dp and BetaRho
            Htop=Z
            P=Pz+Dp
            IF (Htop.NE.Hbot)THEN
              Rho0 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              T=T+(Htop-Hbot)*BetaT
              X=X+(Htop-Hbot)*BetaXfct*X0
              Rho1 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              BetaRho=(Rho1-Rho0)/(Htop-Hbot)
              Dp=Dp+psz(Pbz+p,Rho0,BetaRho,Hbot,Htop,G)
            ENDIF
            Rhodr = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)
            Rho = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)
            RETURN

          ELSE
            ! bottom of the layer is below Z  (Z above ref)
            Htop=H
            ! P is the pressure up to the start height of the layer we just reached
            P=Pz+dp
            IF (Htop.NE.Hbot)THEN
              Rho0 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              T=T+(Htop-Hbot)*BetaT
              X=X+(Htop-Hbot)*BetaXfct*X0
              Rho1 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              BetaRho=(Rho1-Rho0)/(Htop-Hbot)
              Dp=Dp+psz(Pbz+p,Rho0,BetaRho,Hbot,Htop,G)
            ENDIF

            Rhodr = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)
            Rho = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)

            ! place current values Hbot and Beta's
            Hbot=H
            BetaT=0.0d0
            BetaXfct=0.0d0
            BetaCfct=0.0d0
            L=L+9
            ilayptr = 0
            if (zone .eq. 0) ilayptr = 9
            IF (L.GE.ilayptr) THEN
              H=Z+1.0d0
            ELSE
              H=0.0d0
            ENDIF
          ENDIF
        END DO

      ELSE
        ! This is the ELSE for negative linkheights Z below the refplane
        H=0.0d0
        BetaT=0.0d0
        BetaXfct=0.0d0
        BetaCfct=0.0d0
        BetaRho=0.0d0
        Htop=0.0d0
        DO WHILE (H.GT.0.0d0)
          ! loop until H<0 ; The start of the layer is below the zone refplane
          L=L-9
          ilayptr = 0
          if (zone .eq. 0) ilayptr = 1
          IF (L.LT.ilayptr) THEN
            ! with H=Z (negative) this loop will exit, no data for interval Z-refplane
            H=Z
            BetaT=0.0d0
            BetaXfct=0.0d0
            BetaCfct=0.0d0
            BetaRho=0.0d0
          ELSE
            H=0.0d0
            BetaT=0.0d0
            BetaXfct=0.0d0
            BetaCfct=0.0d0
          ENDIF
        END DO

        ! The link is in this layer also if it is on the bottom of it.
        DO
          IF (H.LE.Z) THEN
            Hbot=Z
            P=Pz+dp
            IF (Htop.NE.Hbot)THEN
              Rho1 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              T=T+(Hbot-Htop)*BetaT
              X=X+(Hbot-Htop)*BetaXfct*X0
              Rho0 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              BetaRho=(Rho1-Rho0)/(Htop-Hbot)
              Dp=Dp-psz(Pbz+p,Rho0,BetaRho,Hbot,Htop,G)
            ENDIF
            Rhodr = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)
            Rho =  PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)
            RETURN
          ELSE
            ! bottom of the layer is below Z  (Z below ref)
            Hbot=H
            P=Pz+dp
            IF (Htop.NE.Hbot)THEN
              Rho1 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              ! T,X,C calculated for the lower height
              T=T+(Hbot-Htop)*BetaT
              X=X+(Hbot-Htop)*BetaXfct*X0
              Rho0 = PsyRhoAirFnPbTdbW(Pbz+p,t,x)
              BetaRho=(Rho1-Rho0)/(htop-hbot)
              Dp=Dp-psz(Pbz+p,Rho0,BetaRho,Hbot,Htop,G)
            ENDIF
            Rhodr = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)
            Rho = PsyRhoAirFnPbTdbW(Pbz+pz+dp,t,x)

            ! place current values Hbot and Beta's
            Htop=H
            L=L-9
            ilayptr = 0
            if (zone .eq. 0) ilayptr = 1
            IF (L.LT.ilayptr) THEN
              H=Z-1.0d0
              BetaT=0.0d0
              BetaXfct=0.0d0
              BetaCfct=0.0d0
            ELSE
              H=0.0d0
              BetaT=0.0d0
              BetaXfct=0.0d0
              BetaCfct=0.0d0
            ENDIF
          ENDIF
        ! ENDIF H<Z
        END DO
      ENDIF

      RETURN

      END SUBROUTINE LClimb

!*****************************************************************************************
!     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 AirflowNetworkSolver

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