      SUBROUTINE DUCTIN(IDBG)
C
C SP43 DUCT MODEL ROUTINE. JULY84-FEJ
C     ORIGINAL CARRIER CODE HAS BEEN EXTRACTED FROM THE MAIN
C     ROUTINE AND PLACED IN THIS SUBMODEL. MOST OF THE SUBROUTINE
C     PARAMETERS RELATE TO THE ORIGINAL MODEL (IE INPUT IS TRANSFERRED
C     FROM THE MAIN ROUTINE TO THIS ONE). THE BATTELLE DUCT MODEL
C     READS ITS OWN DATA FROM THE INPUT FILE.
C
C THIS SUBMODEL CONSISTS OF THE FOLLOWING ROUTINES
C     1) DUCTIN.....PROCESS INPUT DATA. LOAD COMMON BLOCKS
C     2) RETDUCT....RETURN DUCT MODEL. CONDITIONED SPACE TO FURNACE.
C     3) SUPDUCT....SUPPLY DUCT MODEL. FURNACE TO CONDITIONED SPACE.
C
C
C COMMON BLOCK :
C     /CARDUCT/...TRANSFERS INPUT DATA FOR THE ORIGINAL (CARRIER)
C                 DUCT MODEL TO THE DUCT SUBROUTINES
C     /DUCTS/.....TRANSFERS TEMPERATURE,HEAT FLOW, AND MASS FLOW
C                 INFORMATION BETWEEN THE DUCT SUBROUTINES AND THE
C                 MAIN PROGRAM
C                 ZONE ASSIGNMENTS: 1)ZONE 1, 2)BASEMENT,
C                      3)CRAWLSPACE, 4)ATTIC 1, 5)FURNACE ENCLOSURE,
C                      6)GARAGE, 7)OUTDOORS, 8)ATTIC 2, 9- Z-2,3...
C------------------------------------------------------------------
C *** DUCTIN: DUCT INPUT DATA PROCESSOR
C
CMDK BTUKWH
CMDK CPAIR
CMDK NZN
CMDK NZN7
      LOGICAL FIRST
      REAL TADAVE(NZN)
CMDK BAR
CMDK BLKGS2
CMDK CPDUCTM      
CMDK CARDUCT
CMDK DUCTSC
CMDK DUCTS1
CMDK DUCTS2
CMDK DUCTS4
CMDK FANBLK
CMDK HUMIDC
CMDK TIMEB
C
      DATA FIRST/.TRUE./
      IF(FIRST)THEN
          DO 10 J=1,NZN7
          IF(J.LE.NZN)QDCTBA(J)=0.
          WDUCTZ(J)=0.
   10     QDUCTZ(J)=0.
C     SET CONSTANTS
          C1=U45*BARRAT*REDC
          C2=U45*BARRAT*REDH
C  INITIAL CALL TO BCLSD AND DCTFLW
          IF(SUPDTYP.EQ.'BATTELLE')CALL BCLSD(DTIME,TOTFLOW,TPL,TADAVE,
     +                             QSUP,QSUPB,IDBG)
          ECFMC4=ECFMC2*U45*BARRAT
          FEAFOF=ECFMC2/ECFMC1
          CONS10(1)=ECFMC4
          CONS14(1)=ECFMC4
          C3=BTUKWH/(U45*BARRAT*REDC*CPAIR)
          C4=BTUKWH/(U45*BARRAT*REDH*CPAIR)
          R1=EFSC/ECFC1S
          R2=EFSH/ECFH1S
          ENDIF
C
C  ENTRY FROM DCTFLW
      ENTRY DUCTN(IDBG)
C
C  FLOW CONSTANTS
      EFANC=R1*REDC*ECFMC1
      EFANH=R2*REDH*ECFMH1
      EFANCC(1)=0.
      EFANCC(2)=EFANC*C3/ECFMC1
      EFANCC(3)=EFANH*C4/ECFMH1
      EFANCC(4)=EFANCC(2)
      CFM(1)=ECFMH1*REDH
      CFM(2)=ECFMC1*REDC
C
C *** CALCULATE CONSTANTS FOR DUCT MODELS
C
C        FURNACE AIR FLOW CONSTANTS
C
      CONS9C= (ECFMC1-CFMBP)*C1
      CONS9H= (ECFMH1-CFMBP)*C2
      CONS9S= ECFMC1*C1
      CONS9T= ECFMH1*C2
      CONS10(2)= CONS9C
      CONS10(3)= CONS9H
      CONS10(4)= CONS9C
      DO 110 I=1,4
  110 CONS11(I)=CONS10(I)*CPAIR
      CONS14(2)= CONS9S
      CONS14(3)= CONS9T
      CONS14(4)= CONS9S
C
C           SUPPLY DUCTS
C
      IF(SUPDTYP.EQ.'CARRIER')THEN
          IF(FIRST)THEN
         CONS1=HAD*ASUP*(1.-FSDIB)
         CONS2=HOD*ASUP*(1.-FSDIB)
         CONS3=WSUP*CPDUCTM*(1.-FSDIB)
         CONS4 = 0.
         IF (ABS(CONS3) .GT. 0.) CONS4 = DTIME/CONS3
         CONS1B=HAD*ASUP*FSDIB
         CONS2B=HOD*ASUP*FSDIB
         CONS3B=WSUP*CPDUCTM*FSDIB
         CONS4B = 0.
         IF (ABS(CONS3B) .GT. 0.) CONS4B = DTIME/CONS3B
          ENDIF
         DO 210 I=1,4
            FHADFO=1.0
            IF(I.EQ.1)FHADFO=FEAFOF
            C11=CONS11(I)*(1.0-PSUPL)
            SUPCN(I)= 1.-EXP(-CONS1*FHADFO/C11)
            SUPCNB(I)= 1.-EXP(-CONS1B*FHADFO/C11)
  210    CONTINUE
      ELSE IF(SUPDTYP.EQ.'BATTELLE')THEN
C  THIS IS TO CHECK SUPDTYP INPUT
      ELSE
         WRITE(60,*) 'BAD SUPTYP VALUE IN SUB.DUCTIN...',SUPDTYP
         STOP 'IN SUB.DUCTIN'
      ENDIF
C
C           RETURN DUCTS
C
      IF(RETDTYP.EQ.'CARRIER')THEN
          IF(FIRST)THEN
         CONS5=HAD*ARET*(1.-FRDIB)
         CONS6=HOD*ARET*(1.-FRDIB)
         CONS7=WRET*CPDUCTM*(1.-FRDIB)
         CONS8 = 0.
         IF (ABS(CONS7) .GT. 0.) CONS8 = DTIME/CONS7
         CONS5B=HAD*ARET*FRDIB
         CONS6B=HOD*ARET*FRDIB
         CONS7B=WRET*CPDUCTM*FRDIB
         CONS8B = 0.
         IF (ABS(CONS7B) .GT. 0.) CONS8B = DTIME/CONS7B
          ENDIF
         DO 410 I=1,4
            FHADFO=1.0
            IF(I.EQ.1)FHADFO=FEAFOF
            C11=CONS11(I)*(1.0-PRETL)
            RETCN(I)= 1.-EXP(-CONS5*FHADFO/C11)
            RETCNB(I)= 1.-EXP(-CONS5B*FHADFO/C11)
  410    CONTINUE
      ELSE IF(RETDTYP.EQ.'BATTELLE')THEN
      WRITE(60,*) 'BATTELLE VERSION OF RETURN DUCT MODEL NOT'
     +           ,' IMPLEMENTED'
         WRITE(60,*)  'CARRIER MODEL IS SUFFICIENT'
         STOP 'IN SUB.DUCTIN USE CARRIER RETURN DUCT MODEL'
      ELSE
         WRITE(60,*) 'BAD RETDTYP VALUE IN SUB.DUCTIN...',RETDTYP
         STOP 'IN SUB.DUCTIN'
      ENDIF
      FIRST=.FALSE.
      RETURN
      END
      SUBROUTINE RETDUCT(        I10, TI, IZLVG, TWRT,TWRTB,
     +   TRET,TRETB,QRET,QRETB, FLOW10,FLOW11,FLOW14)
C
C *** RETDUCT:  RETURN DUCT MODEL
C        TI  MEAN AIR TEMP OF ZONES BEING SUPPLIED (SET IN LOOP)
C        IZLVG  LIVING SPACE ZONE WHERE MOST OF RETURN DUCTS ARE
C
CMDK NZN
CMDK NZN7
CMDK CARDUCT
CMDK DUCTSC
CMDK DUCTS1
CMDK DUCTS4
CMDK IZZQ
C
C *** NOTE THAT THE CALLING SEQUENCE FOR THE DUCT MODEL IS
C     RETURN DUCTS, (FURNACE), SUPPLY DUCTS. THE QDUCTZ AND WDUCTZ
C     ARRAYS ARE INITIALIZED IN THIS ROUTINE AND ADDED TO IN SUPDUCT
C
C-----------------------------------------------------------------
C        FLOWS NEEDED BY MAIN PROGRAM ARE DEFINED HERE
      FLOW10=CONS10(I10)
      FLOW11=CONS11(I10)
      FLOW14=CONS14(I10)
C-----------------------------------------------------------------
C        ORIGINAL SP43 RETURN DUCT ROUTINE
      IF(RETDTYP.EQ.'CARRIER')THEN
         TDUCTEB=TZONE(IFLOC)
         TRET=TI+(TWRT-TI)*RETCN(I10)
         C11=CONS11(I10)*(1.0-PRETL)
         TWRTN=CONS8*(C11*(TI-TRET)-CONS6*(TWRT-TI))+TWRT
         TRETB = (TRET+(TWRTB-TRET)*RETCNB(I10))*(1-IBRP)+TDUCTEB*IBRP
         TRETB=((1.0-PRETL)*TRETB+PRETL*TDUCTEB)
         TWRTBN=CONS8B*(C11*(TRET-TRETB)-CONS6B*(TWRTB-TDUCTEB))+TWRTB
         QRET=CONS6*(TWRT-TI)
         QRETB=CONS6B*(TWRTB-TDUCTEB)*(1-IBRP)+
     +               CONS11(I10)*(TI-TDUCTEB)*IBRP
C
C  INITIALIZE DUCT LEAKAGE ENERGY
C
      DO 10 J=1,NZN7
      WDUCTZ(J)=0.
   10 QDUCTZ(J)=0.
C
C   IZLVG IS LVG SPC ZONE WITH THE MAJORITY OF RETURN DUCTS.
C    IT IS SET IN LOOP.
         QDUCTZ(IZLVG)=QRET
         QDUCTZ(IFLOC)=QRETB
         TWRT=TWRTN
         TWRTB=TWRTBN
         RETURN
      ENDIF
C
C-----------------------------------------------------------------
C        BATTELLE RETURN DUCT MODEL
C
      IF(RETDTYP.EQ.'BATTELLE')THEN
      WRITE(60,*) 'BATTELLE VERSION OF RETURN DUCTS WERE NEVER CODED.'
         WRITE(60,*) 'THE CARRIER MODEL IS SUFFICIENT FOR THESE LOW'
         WRITE(60,*) 'TEMPERATURE SIMPLE DUCTS'
         STOP 'SUB.RETDUCT USE CARRIER RETURN DUCT MODEL'
      ENDIF
      STOP 'BAD RETDTYP IN SUB.RETDUCT'
      END
C
      SUBROUTINE SUPDUCT(        I10, TI, IZLVG, TPL, TWSP,TWSPB,
     +    TSUP, TSUPB, QSUP,QSUPB,IDBG)
C
C *** SUPDUCT: SUPPLY DUCT MODEL
C       TI   MEAN AIR TEMP OF ZONES BEING SUPPLIED (SET IN LOOP)
C       IZLVG  LIVING SPACE ZONE WHERE MOST OF RETURN DUCTS ARE
CMDK NZN
CMDK CARDUCT
CMDK DUCTSC
CMDK DUCTS1
CMDK DUCTS4
CMDK IZZQ
CMDK TIMEB
      DIMENSION TSUP(NZN)
C
C------------------------------------------------------------------
C        ORIGINAL SP43 SUPPLY DUCT MODEL
C           (THIS MODEL IS SET FOR SINGLE-ZONE USE)
C
      IF(SUPDTYP.EQ.'CARRIER')THEN
         TODDB=TZONE(7)
         TDUCTEB=TZONE(IFLOC)
         TSUPB=TPL+(TWSPB-TPL)*SUPCNB(I10)
         C11=CONS11(I10)*(1.0-PSUPL)
         TWSPBN=CONS4B*(C11*(TPL-TSUPB)-CONS2B*(TWSPB-TDUCTEB))+TWSPB
         TSUP(1)=TSUPB+(TWSP-TSUPB)*SUPCN(I10)
         TWSPN=CONS4*(C11*(TSUPB-TSUP(1))-CONS2*(TWSP-TI))+TWSP
C
C        EVALUATE EFFECTIVE TEMPERATURES FOR DUCT LEAKAGE
         IF(PSUPL.LT.PRETL) THEN
C           VACUUM ON BASEMENT, PRESSURIZING UPSTAIRS
            TEFFN=TI
            TEFFL=PEDLO*TODDB+(1.0-PEDLO)*TI
         ELSE IF(PSUPL.GT.PRETL) THEN
C           PRESSURIZING BASEMENT, VACUUM UPSTAIRS
            TEFFN=PEDLO*TODDB + (1.0-PEDLO)*TDUCTEB
            TEFFL=TDUCTEB
         ELSE
C           BALANCED LEAKAGE
            TEFFN=0.
            TEFFL=0.
         ENDIF
         QNET(1)=CONS11(I10)*( (1.0-PSUPL)*TSUP(1)- (1.0-PRETL)*TI +
     +                      (PSUPL-PRETL)*TEFFN )
         QLEAK=CONS11(I10)*(PSUPL*TPL-PRETL*TDUCTEB-(PSUPL-PRETL)*TEFFL)
         QSUP=CONS2*(TWSP-TI)
         QSUPB=CONS2B*(TWSPB-TDUCTEB)
         QDUCTZ(IZLVG)= QDUCTZ(IZLVG)+ QSUP
         QDUCTZ(IFLOC)= QDUCTZ(IFLOC)+ QLEAK+QSUPB
         TWSP=TWSPN
         TWSPB=TWSPBN
         RETURN
      ENDIF
C
C------------------------------------------------------------------
C        BATTELLE SUPPLY DUCT MODEL
C
      IF(SUPDTYP.EQ.'BATTELLE')THEN
         TOTFLOW=CONS10(I10)
          QSUP=0.
          QSUPB=0.
         CALL BCLSD(DTIME,TOTFLOW,TPL,TSUP,QSUP,QSUPB,IDBG)
         RETURN
      ENDIF
      STOP 'BAD SUPDTYP IN SUB.SUPDTYP'
      END
      SUBROUTINE BCLSD(DT,TOTFLOW,TPL,TADAVE,QSUP,QSUPB,IDBG)
C
C BCLSD= (B)ATTELLE (C)OLUMBUS (L)ABS (S)UPPLY (D)UCT SUBMODEL
C
C SUPPLY DUCT MODEL SUBROUTINE
C        BHF ORIGINAL MODIFIED BY FEJ (AUG83, JUL84, SEPT84)
C                     MODIFIED BY MTL (SEPT86)
C                     DEBUGGED BY KEH (OCT 86)
C                     MODIFIED FOR MULTIZONE BY RDF (APR 87)
C
C  TZONE(I) DEFINITION (SET IN LOOP AT LINE LOOP117)
C  I= 1  ZONE NO. 1 AIR TEMP
C     2  BASEMENT AIR TEMP
C     3  CRWLSPC AIR TEMP
C     4  ATTIC NO. 1 AIR TEMP
C     5  FURNACE ENCLOSURE AIR TEMP
C     6  GARAGE AIR TEMP
C     7  OUTDOOR AIR TEMP
C     8  ATTIC NO. 2 AIR TEMP
C     9-12  ZONES 2-5 AIR TEMPS
C
CMDK CPAIR
CMDK MAXNODE
CMDK MAXSECT
CMDK NWL
CMDK NZN
CMDK NZW
CMDK RHOAIR
CMDK ZERO
      PARAMETER (EXP1=0.8,EXP2=0.2,SMALL=0.01)
C
CMDK BLKICL
CMDK DUCTS1
CMDK DUCTS2
CMDK DUCTS3
CMDK ENCBK1
CMDK TIMEB
      REAL  A(MAXSECT), D(MAXSECT), DLLAST(MAXSECT), FP1(MAXSECT),
     +      FP2(MAXSECT), HE1(MAXSECT), HE2(MAXSECT),
     +      MCPD(MAXSECT), PCTDZ(16), PCTL1(MAXSECT),
     +      PCTL2(MAXSECT), PE(MAXSECT), PI(MAXSECT), RE(MAXSECT),
     +      RI(MAXSECT), TAD(MAXSECT,MAXNODE), TADAVE(NZN),
     +      TDM1(MAXSECT,MAXNODE), TDM2(MAXSECT,MAXNODE),
     +      XLS(MAXSECT)
      CHARACTER*3 IEND
      INTEGER I2S(MAXSECT), I2SI(MAXSECT), ITYPIN(MAXSECT),
     +      IZNFLR(MAXSECT), NFEED(MAXSECT),
     +      NODES(MAXSECT), NZONE1(MAXSECT), NZONE2(MAXSECT)
      DIMENSION FLWW(5)
      LOGICAL LPRINT,FIRST,IREAD
C
      DATA FIRST/.TRUE./, LPRINT/.FALSE./ ,IREAD/.TRUE./
      DATA PCTDZ/16*1./
      DATA I2S/MAXSECT*0/
C-------------------------------------------------------------
C     FIRST TIME THROUGH (CALLED FROM SUB.DUCTIN)
C     READ DATA, CONVERT PERCENT TO DECIMALS, DETERMINE NO. OF NODES
C     AND LENGTH OF LAST NODE, INITIALIZE TEMPERATURES
C
      IF(.NOT.FIRST)GO TO 55
         FIRST=.FALSE.
         DTIME=DT
         OPEN(2,FILE='TAPE2',STATUS='OLD',IOSTAT=IO2)
         IF(IO2.NE.0)THEN
           WRITE(60,*)' BCLSD: CANT OPEN TAPE2 WITH DUCT DATA'
           STOP ' BCLSD: CANT OPEN TAPE2 WITH DUCT DATA!'
           ENDIF
C
C  NSECT  - TOTAL NO. OF DUCT SECTIONS IN THE SUPPLY DUCTS
C  DL     - LENGTH OF EACH NODE IN EACH DUCT SECTION, FT
C  NSCT2S - TOTAL NO. OF DUCT SECTIONS WITH TWO SURFACES HAVING
C           HEAT LOSS TO SEPARATE ZONES
C  HIFAC  - MULTIPLIER ON CONVECTIVE COEF INSIDE DUCT
C           (IF INPUT AS NEG. VALUE, ORIG H AS IN SINGLE-ZONE MODEL
C            WILL BE USED)
C  ILEAK  - >0 TO READ PCTDZ VALUES (PERCENTAGE MULTIPLIERS ON
C           LEAKAGE FLOWS SPECIFIED BY PCTL1 AND PCTL2.  PCTDZ VALUES
C           ARE FOR EACH OF THE POSSIBLE ZONE DAMPER CONFIGURATIONS,
C           WHICH IS = ILEAK--WHICH IS THE SAME AS NSUM IN DCTFLW)
         READ(2,*,END=900) NSECT,DL,NSCT2S,HIFAC,ILEAK
         IOLDH=0
         IF(HIFAC.LT.0.)THEN
            IOLDH=1
            HIFAC=ABS(HIFAC)
            WRITE(60,500)
         ELSE
            WRITE(60,501)
            ENDIF
         IF(ILEAK.GT.0)THEN
C  PCTDZ - SEE ILEAK ABOVE         
            READ(2,*)(PCTDZ(I),I=1,ILEAK)
            DO 2 I=1,ILEAK
   2        PCTDZ(I)=PCTDZ(I)*.01
            WRITE(60,1040)(PCTDZ(I),I=1,ILEAK)
            ENDIF
         WRITE(60,1000) NSECT,DL,NSCT2S,HIFAC
C
         IF(NSECT.GT.MAXSECT) GO TO 900
C
C !!REQD ORDER OF DUCT SECTIONS READ ON TAPE2::!!!!!
C  ENTER DATA FOR A DUCT SET FOLLOWED BY THAT FOR ANOTHER SET, IF ANY.
C  A DUCT SET CONSISTS OF TRUNK DUCTS FOLLOWED BY BRANCH DUCTS THAT
C  CONNECT TO THESE TRUNK DUCTS.  ENTER TRUNK DUCTS IN ORDER OF
C  CONNECTION, WITH TRUNK DUCT THAT IS CONNECTED TO PLENUM ENTERED
C  FIRST.  NUMBER ISECT = SAME AS ORDER OF READING ON TAPE2, FOR
C  CONVENIENCE.  IF YOU DONT OBSERVE THIS ORDER, CALCULATION OF A
C  DOWNSTREAM SEGMENT WILL USE UPSTREAM AIR TEMPS FROM THE PREVIOUS
C  TIMESTEP RATHER THAN FROM THE PRESENT TIMESTEP.
C    I FROM 1 TO NSECT!!
C  ISECT(I)  - DUCT SECTION OR INDEX NO.(SEE ABOVE)
C  NFEED(I)  - NUMBER OR INDEX OF DUCT SECTION FEEDING THIS DUCT SECTION
C  NZONE1(I) - ZONE IN WHICH DUCT IS LOCATED(SEE DEF. OF I FOR TZONE ABOVE)
C              USUALLY WILL BE 2 BUT COULD BE 1 TO 12!
C  IZNSUP(I) - INDICATOR WHERE AIR FROM THIS DUCT SECTION GOES
C              =0 INTO ANOTHER DUCT SECTION;
C              =I INTO LVG SPC ZONE I (I=1 TO 5)
C  IZNFLR(I) - INDICATOR WHERE HEAT LOSS FROM THIS DUCT SECTION GOES:
C              =0  NOT UP THROUGH FLOOR ELEMENT
C              =I  UP THROUGH FLOOR ELEMENT INTO ZONE I (I=1 TO 5)
C  ITYPIN(I) - =0, NO X/D CORRECTION IN TURBULENT FLOW
C              =1, SQ EDGED ORIFICE AT INLET TO DUCT SECTION
C              =2, 90 DEG ELBOW AT INLET (BOTH 1 AND 2 ARE FOR TURB FLOW)
C  PCTL1(I) - NOTE...PCTL1 IS LEAKAGE FLOW IN PCT OF FLOW
C           IN DUCT RATHER THAN PCT OF TOTAL FURNACE FLOW!!!!!!! ASSIGN
C           PCTL1 TO SECTION OF DUCT RUN CONNECTING TO REGISTER
C           AND NOT TO ANY UPSTREAM DUCT RUNNER SECTIONS. LEAKAGE WILL
C           BE CALCULATED AT FIRST NODE OF THIS ASSIGNED DUCT RUN.
C  XLS(I) - LENGTH OF DUCT SECTION I,FT
C  D(I)   - INTERNAL DIAMETER OF DUCT SECTION I, FT
C  A(I)   - FLOW AREA OF DUCT SECTION I, FT2
C  PI(I)  - INTERIOR PERIMETER OF DUCT SECTION I, FT
C  PE(I)  - EXTERIOR PERIMETER OF DUCT SECTION I, FT
C  RI(I)  - THERMAL RESISTANCE OF INSULATION ON INSIDE OF DUCT SECTION I,
C           FT2-F-HR/BTU
C  RE(I)  - THERMAL RESISTANCE OF INSULATION ON OUTSIDE OF DUCT SECTION I,
C           FT2-F-HR/BTU
C  HE1(I) - EXTERIOR CONVECTIVE HEAT TRANSFER COEF FOR DUCT SECTION I,
C           BTU/HR-FT2-F
C  MCPD(I)- PRODUCT OF MASS AND CP OF DUCT SECTION I, BTU/F
C  
         DO 30 I=1,NSECT
            IF(NSCT2S.EQ.0) THEN
               READ(2,*,END=900) ISECT(I),NFEED(I),NZONE1(I),
     +            IZNSUP(I), IZNFLR(I), ITYPIN(I), PCTL1(I), XLS(I),
     +            D(I),A(I),PI(I),PE(I),RI(I),RE(I),HE1(I),MCPD(I)
               FP1(I)=1.
               FP2(I)=0.
               HE2(I)=HE1(I)
               PCTL2(I)=0.
               NZONE2(I)=NZONE1(I)
            ELSE
C
C    INDENTIFY SECTIONS WITH 2 DUCT WALL SURFACES
               IF(IREAD) THEN
                  IREAD=.FALSE.
                  IF(NSCT2S.NE.0)THEN
C I2SI IS DUCT SECT NO. OF DUCTS EXPOSED TO TWO ENVIRONMENTS
                    READ(2,*,END=900) (I2SI(K), K=1, NSCT2S)
                    DO 10 K=1,NSCT2S
                    I2SQ=I2SI(K)
   10               I2S(I2SQ)=1
                    ENDIF
               ENDIF
C
                  IF(I2S(I).NE.0) THEN
C   FP1, FP2 ARE FRACTION OF PERIMETER EXPOSED TO TZONE(NZONE1)
C            AND TZONE2(NZONE2), RESPECTIVELY
C   PCTL2(I) - SEE PCTL1(I) ABOVE BUT FOR LEAKAGE TO 2ND ENVIRONMENT
C   HE2(I)  - SEE HE1(I) ABOVE BUT FOR HEAT TRANSFER COUPLING WITH
C             2ND ENVIRONMENT
                     READ(2,*,END=900) ISECT(I), NFEED(I), NZONE1(I),
     +                  NZONE2(I), IZNSUP(I), IZNFLR(I), ITYPIN(I),
     +                  PCTL1(I), PCTL2(I), XLS(I), D(I),
     +                  A(I), PI(I), PE(I), FP1(I), FP2(I),
     +                  RI(I), RE(I), HE1(I), HE2(I), MCPD(I)
                  ELSE
                     READ(2,*,END=900) ISECT(I),NFEED(I),NZONE1(I),
     +                  IZNSUP(I), IZNFLR(I), ITYPIN(I), PCTL1(I),
     +                  XLS(I),D(I),A(I), PI(I),PE(I),RI(I),RE(I),
     +                  HE1(I),MCPD(I)
                     FP1(I)=1.
                     FP2(I)=0.
                     HE2(I)=HE1(I)
                     PCTL2(I)=0.
                     NZONE2(I)=NZONE1(I)
                  ENDIF
            ENDIF
            WRITE(60,1010)ISECT(I), NFEED(I), NZONE1(I), NZONE2(I),
     +         IZNSUP(I), IZNFLR(I),ITYPIN(I), PCTL1(I), PCTL2(I),
     +         XLS(I),D(I), A(I), PI(I), PE(I), FP1(I), FP2(I),
     +         RI(I), RE(I), HE1(I), HE2(I), MCPD(I)
            PCTL1(I)=PCTL1(I)*0.01
            PCTL2(I)=PCTL2(I)*0.01
C  NODES(I) IS NO. OF DL-LENGTH NODES PER DUCT SECTION
            NODES(I)=INT(XLS(I)/DL)
C  DLLAST(I) IS LENGTH OF LAST NODE IN DUCT SECTION I
            DLLAST(I)= XLS(I) - NODES(I)*DL
            IF(DLLAST(I).LT. SMALL ) THEN
               DLLAST(I)=DL
            ELSE
               NODES(I)=NODES(I)+1
            ENDIF
            IF(NODES(I).GT.MAXNODE) THEN
               NODES(I)=MAXNODE
               DLLAST(I)= XLS(I) - NODES(I)*DL
               WRITE(60,*) ' BCLSD: WARNING! SECTION ',I,' IS TOO LONG.'
              WRITE(60,*) ' THE LAST SECTION IS ',DLLAST(I),' FT LONG!!'
               WRITE(60,*) ' SO, INCREASE VALUE OF PARAMETER MAXNODE'
            ENDIF
            DO 20 J=1,NODES(I)
            NZ1=NZONE1(I)
            TDM1(I,J)=TZONE(NZ1)
            IF(I2S(I).NE.0)THEN
               NZ2=NZONE2(I)
               TDM2(I,J)=TZONE(NZ2)
               TAD(I,J)=0.5*(TZONE(NZ1)+TZONE(NZ2))
            ELSE
               TAD(I,J)=TZONE(NZ1)
               ENDIF
   20       CONTINUE
            HE1(I)=1./(RE(I)+1./HE1(I))
            HE2(I)=1./(RE(I)+1./HE2(I))
   30    CONTINUE
          TPLS=TPL
C  READ LINES UNTIL ENCOUNTER  END  IN 1ST THREE COLUMNS
   32 READ(2,502)IEND
      IF(IEND.NE.'END')GO TO 32          
C  INITIAL CALL TO DCTFLW
          CALL DCTFLW(IDBG)
         GO TO 200
C
C-------------------------------------------------------------
C     FOR EACH TIMESTEP, RUN THROUGH EACH SECTION (CALLED FROM MAIN PGM
C
   55  CONTINUE
          DO 60 I=1,NZN
          TADAVE(I)=0.
          FLWW(I)=0.
   60     QNET(I)=0.
      QLEAK=0.0
C  ICOL IS SET IN DCTFLW
          PCTD=PCTDZ(ICOL)
      DO 120 I=1,NSECT
          IZN=IZNSUP(I)
          NZ1=NZONE1(I)
          NZ2=NZONE2(I)
          NFD=NFEED(I)
          NDS=NODES(I)
          ISCT=ISECT(I)
          DI=D(I)
          ITYPD=ITYPIN(I)
          FLOW=TOTFLOW*PCTF(I)
          IF(FLOW.LT.0.)THEN
             WRITE(60,1050)ISCT,PCTF(I)
             STOP 'BCLSD:FLOW NEGATIVE'
             ENDIF
          IF(NFD.EQ.0)THEN
              NDSP=1
          ELSE
              NDSP=NODES(NFD)
              ENDIF
C
         SUMQD=0.
         XLSUM=0.
         DO 110 J=1,NDS
            QLEAKI=ZERO
            WLEAKI=ZERO
            QLEAK1=ZERO
            WLEAK2=ZERO
C
C           FOR EACH NODE@D SET INLET TEMP = OUTLET OF PREVIOUS NODE
            IF(J.EQ.1) THEN
               IF(NFD.EQ.0) THEN
                  TAD(I,J)=TPL
               ELSE
                  TAD(I,J)=TAD(NFD,NDSP)
               ENDIF
C
C     CALCULATE ENERGY LOSS DUE TO AIR LEAKAGE FROM DUCT TO ZONE
C        (FIRST NODE OF SECTION ONLY)
C
               IF(PCTL1(I).GT.ZERO.OR.PCTL2(I).GT.ZERO) THEN
                  CTD1= TAD(I,J)-TZONE(NZ1)
                  WLEAK1= FLOW*PCTL1(I)*PCTD
                  QLEAK1= WLEAK1*CPAIR*CTD1
                  WLEAK2=0.
                  QLEAK2=0.
                  CTD2=CTD1
                  IF(I2S(I).NE.0) THEN
                     CTD2=TAD(I,J)-TZONE(NZ2)
                     WLEAK2= FLOW*PCTL2(I)*PCTD
                     QLEAK2= WLEAK2*CPAIR*CTD2
                     ENDIF
                  WLEAKI=WLEAK1+WLEAK2
                  FLOW=FLOW-WLEAKI
                  QLEAKI=QLEAK1+QLEAK2
                  QLEAK=QLEAK+QLEAKI
               ELSE
                  QLEAK2=0.
                  WLEAK1=0.
               ENDIF
            ELSE
               TAD(I,J)=TAD(I,J-1)
            ENDIF
C
C     DETERMINE HEAT FLOWS (Q IN BTUH)
C        FOR DUCT SECTIONS WITH 2 HEAT TRANSFER SURFACES, ASSUME
C        HI AND RI ARE THE SAME FOR BOTH SURFACES
C  GET CONVECTIVE HEAT TRANSFER COEFFICIENT FOR INSIDE DUCT SECTION
            TA=TAD(I,J)
               DLEN=DL
               IF(J.EQ.NDS)DLEN=DLLAST(I)
            IF(TOTFLOW.LT.10.)THEN
               HI=0.
            ELSE
               XLSUM=XLSUM+DLEN
               XOD=XLSUM/DI
C        DETERMINE MASS VELOCITY IN DUCT AND CONVECTIVE H
               G=FLOW/A(I)
               CALL HDUCT(ISCT,TA,G,DI,XOD,ITYPD,IOLDH,HI,IDBG)
               HI=HIFAC/(RI(I)+1./HI)
               ENDIF
            TM1=TDM1(I,J)
            QDTA1=HI    *PI(I)*DLEN*FP1(I)*(TM1-TA)
            QDTZ1=HE1(I)*PE(I)*DLEN*FP1(I)*(TM1-TZONE(NZ1))
            IF(I2S(I).EQ.0) THEN
C              SECTIONS WITH A SINGLE DUCT WALL
               QDTA2=0.
               QDTZ2=0.
            ELSE
C              SECTIONS WITH TWO DUCT WALL SURFACES
               TM2=TDM2(I,J)
               QDTA2=HI    *PI(I)*DLEN*FP2(I)*(TM2-TA)
               QDTZ2=HE2(I)*PE(I)*DLEN*FP2(I)*(TM2-TZONE(NZ2))
               ENDIF
            TDM1(I,J)=TM1-(QDTZ1+QDTA1)*DTIME/
     +                (FP1(I)*MCPD(I)*DLEN/XLS(I))
            IF(TDM1(I,J).LT.TZONE(NZ1))THEN
               QDTZ1=0.
               TDM1(I,J)=TZONE(NZ1)
               TM1=TDM1(I,J)
               ENDIF
            IF(FP2(I).NE.0.) THEN
               TDM2(I,J)=TM2-(QDTZ2+QDTA2)*DTIME/
     +                (FP2(I)*MCPD(I)*DLEN/XLS(I))
               IF(TDM2(I,J).LT.TZONE(NZ2))THEN
                  QDTZ2=0.
                  TDM2(I,J)=TZONE(NZ2)
                  TM2=TDM2(I,J)
                  ENDIF
            ELSE
               TDM2(I,J)=TDM1(I,J)
               ENDIF
C
C           DETERMINE THE NEW TEMPERATURES (FOR NEXT TIME STEP)
            IF(HI.EQ.0.)THEN
               TAD(I,J)=TA
            ELSE
               TAD(I,J)=TA+(QDTA1+QDTA2)/(FLOW*CPAIR)
C              IF(FP2(I).EQ.0.)THEN
C                IF(TAD(I,J).LT.TDM1(I,J))THEN
C                   TAD(I,J)=TDM1(I,J)
C                   QDTA1=0.
C                   GO TO 90
C                   ENDIF
C             ELSE
C                IF(TDM1(I,J).LT.TDM2(I,J))THEN
C                   IF(TAD(I,J).LT.TDM1(I,J))THEN
C                      TAD(I,J)=TDM1(I,J)
C                      QDTA1=0.
C                      GO TO 90
C                      ENDIF
C                ELSE
C                   IF(TAD(I,J).LT.TDM2(I,J))THEN
C                      TAD(I,J)=TDM2(I,J)
C                      QDTA2=0.
C                      GO TO 90
C                      ENDIF
C                   ENDIF
C                ENDIF
              ENDIF
C  90      CONTINUE
C
C      SUM HEAT&MASS FLOWS FOR SUPPLY AIR AND EXTERNAL ZONES
C         BTU/HR AND LBM/HR RESPECTIVELY
          SUMQD=SUMQD+QDTZ1+QLEAK1
            QDUCTZ(NZ1)=QDUCTZ(NZ1)+QDTZ1 + QLEAK1
            QDUCTZ(NZ2)=QDUCTZ(NZ2)+QDTZ2 + QLEAK2
            WDUCTZ(NZ1)=WDUCTZ(NZ1)+WLEAK1
            WDUCTZ(NZ2)=WDUCTZ(NZ2)+WLEAK2
C  SUM DUCT HEAT LOSS TO LVG SPC (QSUB) AND TO BSMT(QSUPB)
          IF(NZONE1(I).EQ.1.OR.NZONE1(I).GT.8)QSUP=QSUP+QDTZ1
          IF(NZONE1(I).EQ.2)QSUPB=QSUPB+QDTZ1
  110    CONTINUE
      IF(IZNFLR(I).GT.0)THEN
          IZNF=IZNFLR(I)
          QDCTBA(IZNF)=QDCTBA(IZNF)+SUMQD
          ENDIF
      IF(IZN.NE.0)THEN
          IQZ=1
          IF(IZN.GT.1)IQZ=IZN+7
          DQNET=FLOW*CPAIR*(TAD(I,NDS)-TZONE(IQZ))
          QNET(IZN)=QNET(IZN)+ DQNET
          FLWW(IZN)=FLWW(IZN)+FLOW
          RAT=PCTF(I)/RWSDCT(IZN)
          DTAV=TAD(I,NDS)*RAT
          TADAVE(IZN)=TADAVE(IZN)+DTAV
          ENDIF
  120 CONTINUE
C      WRITE(20,567)TIME*60.,TOTFLOW/20.,TPL,FLWW(1)/20.,QNET(1),
C     +FLWW(2)/20.,QNET(2),FLWW(3)/20.,QNET(3)
C  567 FORMAT(1X,F8.3,2X,F7.2,2X,F6.1,2X,3(F7.2,2X,F8.1,4X))
C-------------------------------------------------------------
C     PRINT OUT RESULTS WHEN LPRINT=TRUE
C
      IF(LPRINT)THEN
         DO 190 I=1,NSECT
          NDS=NODES(I)
              WRITE(60,1030)ISECT(I),TAD(I,1),TDM1(I,2),TDM2(I,2),
     +                     TAD(I,NDS), TDM1(I,NDS), TDM2(I,NDS)
  190    CONTINUE
         TIN=TPL*0.505+TPLS*0.495
C  NOTE::DEFF BELOW ONLY VALID FOR ZONE NO. 1
         DEFF=(TADAVE(1)-TZONE(1))/(TIN-TZONE(1))
         WRITE(60,3000) TADAVE
          WRITE(60,3010)DEFF
      ENDIF
      TPLS=TPL
C     STOP 'BCLSD: STOPPED AT FIRST RETURN'
  200 CONTINUE
      RETURN
C
C-------------------------------------------------------------
C        ERROR EXIT
C
  900 CONTINUE
      WRITE(60,*) 'BCLSD: ??? IN READING DUCT DATA ON TAPE2'
      STOP 'TROUBLE IN SUB.DUCT'
C
  500 FORMAT(1H0,'ORIGINAL TURB CONVECT COEF USED WITHIN DUCTS')
  501 FORMAT(1H0,'LAM AND TURB CONVECT COEF USED WITHIN DUCTS'/
     +'  WITH ENTRANCE LENGTH CORRECTION IF ITYPIN > 0 FOR SECTION')
  502 FORMAT(A3)
 1000 FORMAT('0BATTELLE SUPPLY DUCT MODEL.',I5,' SECTIONS. ',
     +'DIFFERENTIAL LENGTH=',F5.2,' FT.',/,
     +'0NO. OF SECTIONS WITH TWO SURFACES=',I3,/,
     +'0 HIFAC = ',G13.5,' MULTIPLIER ON ALL DUCT INTERIOR H VALUES'//
     +'0ISECT NFED  NZ1  NZ2  SUP  FLR ITYP PCTL1 PCTL2  XLS  D   ',
     +'A   PI    PE    FP1    FP2    RI    RE   HE1   HE2   MCP',
     +'   HIFAC')
 1010 FORMAT(2I6,5I4,F6.1,2F6.2,F6.1,F5.2,F6.3,2F6.2,2F7.4,4F6.3,F7.1,
     +       F7.1)
 1020 FORMAT('0FOR ZONE ',I3,' FLOW FRACTION SUPPLIED IS=',G13.5)
 1030 FORMAT(1X,'ISECT= ',I4,' INLET TEMPS(DEG F)--AIR,METAL(1),',
     +          'METAL(2)=',3G13.5/
     +12X,'OUTLET TEMPS(DEG F)--AIR,METAL(1),METAL(2)=',3G13.5)
 1040 FORMAT(1X,'MULTIPLIERS ON LEAKAGE FLOWS FOR EACH SETTING OF '
     +,'ZONE DAMPERS= '/5X,16F6.3)
 1050 FORMAT(1X,'BCLSD: FLOW FOR DUCT NO. ',I5,' NEG. PCTF(I)= ',
     +G13.5)
 3000 FORMAT(1X,'PREDICTED DUCT EXIT T(AVE FOR EA ZONE)= ',5F6.1)
 3010 FORMAT(1X,'PREDICTED DUCT EFFICIENCY  =   ',F6.3)
 3020 FORMAT(1H0,'SEE SUGGESTED ORDER OF DUCT-SECTION DATA IN ',
     +'LISTING FOR SUBPROGRAM BCLSD')
      END
