CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C  TYPE 202 - MULTIPLE-ZONE INDOOR CONTAMINANT TRANSPORT
C
C  CALCULATES THE CONCENTRATION OF A CONTAMINANT IN A MULTIPLE-ZONE 
C  VENTILATED VOLUME AS A FUNCTION OF SIMULATION TIME
C 
C  SAMPLE INPUT AND OUTPUT FOLLOW THE NOMENCLATURE
C
C  CREATED BY:
C  	PAUL KNOESPEL
C 
C  CREATED ON:
C	MAR  2, 1990
C 
C  MODIFIED ON / BY / REASON:
C
C
C  SUBROUTINE TYPE60 CALLS:
C	SUBROUTINE SHORTCALC
C	SUBROUTINE CALC
C
C  VARIABLE LISTING
C
C  I - INPUT VALUE
C  C - CALCULATED VALUE
C
C	A1	    C	ZONE SOURCES INDEPENDENT OF ZONE CONCENTRATION
C	A2	    C	ZONE SOURCES DEPENDENT ON ZONE CONCENTRATION
C	A3	    C	ZONE SINKS DEPENDENT ON ROOM RETURN DUCT CONCENTRATION
C	C	    C	ZONE AVERAGE CONTAMINANT VOLUME CONCENTRATION
C	C1	    C	ITERATIVE CONCENTRATION VALUE AT TIME STEP I+1
C	CAHU	    C	COMBINED RETURN CONTAMINANT CONCENTRATION AT THE AHU
C	COA	    I	CONTAMINANT CONCENTRATION IN OUTSIDE AIR
C	CRET	    C	RETURN DUCT CONTAMINANT CONCENTRATION AT THE ZONE
C	CSUP	    C	SUPPLY DUCT CONTAMINANT CONCENTRATION 
C	DELTAT	    I	TRNSYS SIMULATION TIME STEP
C	DT	    I	TIME STEP FOR EACH CALCULATION (S)
C	DTDT	    C	TRNSYS VARIABLE - NOT USED
C	ERR	    C	ERROR TOLERANCE FOR CHANGE BETWEEN ITERATION STEPS
C	ERRFLAG	    C	ERROR FLAG SET WHEN ERROR TOLERANCE IS NOT MET
C	F	    C	FUNCTION FOR RHS OF CONCENTRATION DIFFERENTIAL EQUATION
C	FILTEFF	    I	FRACTIONAL FILTER EFFICIENCY
C	HOD	    C	CURRENT HOUR OF THE DAY BASED ON SIMULATION TIME
C	HOW	    C	CURRENT HOUR OF THE WEEK BASED ON SIMULATION TIME
C	I	    C	TIME STEP INDEX
C	IDELAY	    C	NUMBER OF TIME STEPS FOR AIR TO FLOW THROUGH RETURN PATH
C	INFO	    C	TRNSYS TYPE INFORMATION
C	IZF	    I	INTER-ZONE FLOW (KG/HR, LBM/HR)
C	J	    C	ZONE NUMBER INDEX
C	K	    C	SOURCE SCHEDULE INDEX
C	L	    C	ZONE NUMBER INDEX IN 2D IZF(J,L) ARRAY
C	LUB	    I	FORTRAN LOGICAL UNIT FOR BUILDING DATA FILE
C	LUS	    I	FORTRAN LOGICAL UNIT FOR SOURCE DATA FILE
C	ND	    C	NUMBER OF DERIVATIVES
C	NI	    C	NUMBER OF INPUTS
C	NP	    C	NUMBER OF PARAMETERS
C	NSTEPS	    C	NUMBER OF TIME STEPS
C	NZONES	    I	NUMBER OF ZONES
C	OUT	    C	OUTPUT ARRAY TO TRNSYS
C	PAR	    I	CONSTANT INPUT ARRAY FROM TRNSYS
C	QCIZF	    C	TOTAL INTER-ZONE CONTAMINANT FLOW INTO A ZONE (ACS)
C	QIZF	    C	INTER-ZONE AIR FLOW (M3/S, FT3/S)
C	QIZFNET	    C	NET INTER-ZONE FLOW FOR A ZONE (M3/S, FT3/S)
C	QZIFOUT	    C	INTER-ZONE FLOW LEAVING A ZONE (M3/S, FT3/S)
C	QOA	    I	OUTSIDE AIR VOLUME FLOW RATE INTO A ZONE (M3/S, FT3/S)
C	QOAPU	    C	ZONE FLOW RATE OF OUTSIDE AIR PER SOURCE UNIT (L3/S-U,
C			FT3/MIN-U)
C	QOATOT	    C	TOTAL OUTSIDE AIR VOLUME FLOW RATE (M3/S, FT3/S)
C	QRET	    C	VENTILATION VOLUME FLOW RATE LEAVING A ZONE (M3/S, FT3/S)
C	QRETTOT	    C	TOTAL RETURN DUCT VOLUME FLOW RATE (M3/S, FT3/S)
C	QSUP	    I	CIRCULATION AIR VOLUME FLOW RATE INTO A ZONE (M3/S, FT3/S)
C	QCIRC	    C	TOTAL CIRCULATION AIR VOLUME FLOW RATE (M3/S, FT3/S)
C	RE	    I	ZONE POLLUTANT REMOVAL EFFECTIVENESS WITH HVAC ON
C	REMEFF	    C	ZONE REMOVAL EFFECTIVENESS USED IN CALCULATIONS
C	RETVOL	    I	VOLUME OF THE ZONE'S RETURN AIR PATH TO THE AHU (M3, FT3)
C	RHOAIR	    C	CONSTANT AIR DENSITY (KG/M3, LBM/FT3)
C	SOURCE	    I	CONTAMINANT SOURCE VOLUME FLOW RATE IN A ZONE (M3/S, FT3/S)
C	SMULT	    C	SOURCE MULTIPLIER FOR THE CURRENT TIME STEP
C	SMULTWD	    I	WEEKDAY SCHEDULED SOURCE MULTIPLIER
C	STIMEWD	    I	WEEKDAY SCHEDULE TIME THAT SOURCE MULTIPLIER IS APPLIED
C	SMULTWE	    I	WEEKEND SCHEDULED SOURCE MULTIPLIER
C	STIMEWE	    I	WEEKEND SCHEDULE TIME THAT SOURCE MULTIPLIER IS APPLIED
C	SUNIT	    I	UNIT CONTAMINANT SOURCE VOLUME FLOW RATE (M3/S, FT3/S)
C	T	    C	TRNSYS VARIABLE - NOT USED
C	TIME	    C	SIMULATION TIME (HR)
C	TIMEPREV	    C	TIME OF THE PREVIOUS SIMULATION STEP
C	UNITFLAG	    I	FLAG TO INDICATE SI OR ENGLISH UNITS (1=SI, 2=ENGLISH)
C	XIN	    I	INPUT ARRAY FROM TRNSYS
C	ZONEVOL	    C	ZONE VOLUME (M3, FT3)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       SAMPLE INPUT AND OUTPUT VALUES
C
C       PARAMETERS              SAMPLE VALUES
C       PAR(1)  NZONES              2.0
C       PAR(2)  DT(S)              30.0
C       PAR(3)  DELTAT(HR)           .25
C       PAR(4)  COA                  .000310
C       PAR(5)  LUS                13
C       PAR(6)  LUB                14
C       PAR(7)  FILTEFF1            0.011
C       PAR(8)  FILTEFF2            0.002
C       PAR(9)  FILTEFF3            0.040
C       PAR(10) UNITS               2.0
C
C       INPUTS
C       XIN(1)                      1.7
C       XIN(2)                      2.4
C       XIN(3)                      3.6
C       XIN(4)                      0.5
C
C       BUILDING DATA FILE (LOGICAL UNIT 14)
C
C       119000.0	2850.0          : ZONEVOL(J)
C	0.9		1.0             : REMEFF(J)
C	21750.0		507.0           : RETVOL(J)
C	0.2		0.0             : INFILT(J)
C                                       :
C	0.0		0.0             : IZF(1,1)      ... IZF(1,NZONES)
C	107.0		0.0             : IZF(NZONES,1) ... IZF(NZONES,NZONES)
C
C       SOURCE DATA FILE (LOGICAL UNIT 13)
C
C	1.77E-4                                       : SUNIT(J)
C	6.0	8.0	12.0	13.0	17.0	19.0  : STIMEWD(J,1-6)
C	20.0	100.0	50.0	100.0	30.0	0.0   : SMULTWD(J,1-6)
C	6.0	8.0	12.0	13.0	17.0	19.0  : STIMEWE(J,1-6)
C	0.0	10.0	5.0	10.0	0.0	0.0   : SMULTWE(J,1-6)
C
C       1.77E-4	                                      : SUNIT(J)
C       8.0	9.0	12.0	13.0	16.0	17.0  : STIMEWD(J,1-6)
C       2.0	10.0	0.0	10.0	3.0	0.0   : SMULTWD(J,1-6)
C       6.0	8.0	12.0	13.0	17.0	19.0  : STIMEWE(J,1-6)
C	0.0	0.0	0.0	0.0	0.0	0.0   : SMULTWE(J,1-6)
C
C       OUTPUTS          SAMPLE VALUES
C       NOTE: THIS IS THE LAST SET OF OUTPUTS AFTER A NINE HOUR SIMULATION
C       WITH A ONE HOUR TIMESTEP.  THIS SET OF INPUTS HAD NO ZERO VALUES.
C       OUT(1)              21.4
C       OUT(2)              11.88
C       OUT(3)               4.504E-4
C       OUT(4)               3.721E-4
C       OUT(5)               4.053E-4
C       OUT(6)               3.721E-4
C       OUT(7)               3.781E-4
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

!  This component has been assigned Type Number 202. If that number conflicts with 
!  another user Type number, you will need to change it and recompile the appropriate
!  dll.

	SUBROUTINE TYPE202(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)

!DEC$ATTRIBUTES DLLEXPORT :: TYPE202



	IMPLICIT NONE

	INTEGER INFO(15), NI, NP, ND, LUB, LUS, LUR, LUW, IFORM, LUK
	INTEGER NSTEPS, NZONES, IDELAY(25), I, J, K, L, UNITFLAG

	DOUBLE PRECISION XIN(50),OUT(77)
	REAL TIME, T, DTDT, ICNTRL, PAR(10)
	REAL C(-60:60,25), CAHU, CSUP, COA, FILTEFF(3)
	REAL SOURCE(25), SUNIT(25), SMULT(25), RHOAIR
	REAL STIMEWD(25,6), SMULTWD(25,6), STIMEWE(25,6), SMULTWE(25,6)
	REAL ZONEVOL(25), RE(25), REMEFF(25), RETVOL(25), DT, DELTAT
	REAL QSUP(25), QRET(25), QOA(25), INFILT(25), IZF(25,25) 
	REAL QCIRC, QOATOT, QRETTOT, QCRETTOT, QOAPU
	REAL QIZF(25,25), QIZFNET, QIZFOUT, TIMEPREV, HOW, HOD
	REAL F, A1(25), A2(25), A3(25)

	COMMON A1, A2, A3
	COMMON /LUNITS/LUR,LUW,IFORM,LUK


!    Set the version information for TRNSYS
      IF (INFO(7).EQ.-2) THEN
	  INFO(12) = 15
	  RETURN 1
	ENDIF


C
C  SET PARAMETERS, DO CHECKS, SET ALL INITIAL ZONE CONCENTRATIONS EQUAL
C  TO THE OUTSIDE AIR CONCENTRATION, AND READ IN SOURCE AND ZONE INFORMATION 
C  ON THE FIRST CALL OF THE SIMULATION
C

	IF(INFO(7).EQ.-1)THEN
	  NZONES = NINT(PAR(1))
	  DT = PAR(2)
	  DELTAT = PAR(3)
	  COA = PAR(4)
	  LUS = NINT(PAR(5))
	  LUB = NINT(PAR(6))
	  FILTEFF(1) = PAR(7)
	  FILTEFF(2) = PAR(8)
	  FILTEFF(3) = PAR(9)
	  UNITFLAG = NINT(PAR(10))

	  NSTEPS = NINT(DELTAT*3600.0/DT)
	  TIMEPREV = -1.0

C
C  USE ENGLISH UNITS IF THE UNITS FLAG IS SET; OTHERWISE USE SI UNITS
C

	  IF (UNITFLAG .EQ. 2) THEN
	    RHOAIR = 0.07517
	  ELSE
	    RHOAIR = 1.204
	  ENDIF

C
C  STOP THE SIMULATION IF NSTEPS IS GREATER THAN 60
C

	  IF ((NSTEPS .LT. 1) .OR. (NSTEPS .GT. 60)) THEN
	    WRITE (LUW,*) 'THE TRNSYS SIMULATION TIME STEP PAR(3) IN'
	    WRITE (LUW,*) 'TYPE 60 MUST BE BETWEEN 1 AND 60 TIMES THE'
	    WRITE (LUW,*) 'TYPE 60 TIME STEP PAR(2)'
	    WRITE (LUW,*) '-- SIMULATION STOPPED --' 
	    CALL MYSTOP(1001)
		RETURN 1
	  ENDIF

	  NI = 2*NZONES
	  NP = 10
	  ND = 0
	  INFO(6) = 3*NZONES+2
	  INFO(9) = 1
	  CALL TYPECK(1,INFO,NI,NP,ND)

	  DO 5 J = 1, NZONES
	    DO 5 I = -60, 60
	      C(I,J) = COA
5	  CONTINUE
 

C
C  READ INPUT VALUES FOR THE SOURCE AND ZONE PHYSICAL DATA
C

	  OPEN (LUS, STATUS = 'OLD')
	  REWIND (LUS)
	  DO 10 J = 1, NZONES
	    READ (LUS,*) SUNIT(J)
	    READ (LUS,*) (STIMEWD(J,K), K = 1, 6)
	    READ (LUS,*) (SMULTWD(J,K), K = 1, 6)
	    READ (LUS,*) (STIMEWE(J,K), K = 1, 6)
	    READ (LUS,*) (SMULTWE(J,K), K = 1, 6)
10	  CONTINUE
	  CLOSE (LUS)

	  OPEN (LUB, STATUS = 'OLD')
	  REWIND (LUB)
	  READ (LUB,*) (ZONEVOL(J), J = 1, NZONES)
	  READ (LUB,*) (RE(J), J = 1, NZONES)
	  READ (LUB,*) (RETVOL(J), J = 1, NZONES)
	  READ (LUB,*) (INFILT(J), J = 1, NZONES)

	  DO 20 J = 1, NZONES
	  READ (LUB,*) (IZF(J,L), L = 1, NZONES)
	    DO 15 L = 1, NZONES
	      QIZF(J,L) = IZF(J,L) / (3600.0 * RHOAIR)
15	    CONTINUE
20	  CONTINUE
	  CLOSE (LUB)

	  ENDIF

C
C  IF THE TIME HAS NOT CHANGED SINCE THE PREVIOUS SUBROUTINE CALL, RETURN TO
C  TRNSYS TO PREVENT CHANGING THE CONCENTRATIONS TWICE DURING A SINGLE TIME STEP
C

	IF (TIME .EQ. TIMEPREV) THEN
	  RETURN 1
	ELSE
	  TIMEPREV = TIME
	ENDIF

C
C  FIND THE CURRENT SOURCE VALUE BASED ON THE SOURCE MULTIPLIER APPROPRIATE
C  FOR THE CURRENT SIMULATION TIME; IF THERE IS NO SOURCE IN A ZONE, SET THE 
C  ZONE REMOVAL EFFECTIVENESS TO 1.0
C

	HOW = AMOD(TIME,168.0)
	HOD = AMOD(HOW,24.0)


	DO 25 J = 1, NZONES

	  IF (HOW .LE. 120.0) THEN
	    SMULT(J) = SMULTWD(J,6)
	    DO 22 K = 1, 5
	      IF ((HOD .GT. STIMEWD(J,K)) . AND.
     .	      (HOD .LE. STIMEWD(J,K+1))) THEN
	        SMULT(J) = SMULTWD(J,K)
	        GOTO 24
	      ENDIF
22	    CONTINUE

	  ELSE
	    SMULT(J) = SMULTWE(J,6)
	    DO 23 K = 1, 5
	      IF ((HOD .GT. STIMEWE(J,K)) .AND.
     .	      (HOD .LE. STIMEWE(J,K+1))) THEN
	        SMULT(J) = SMULTWE(J,K)
	        GOTO 24
	      ENDIF
23	    CONTINUE

	  ENDIF

24	  SOURCE(J) = SUNIT(J)*SMULT(J)
	  IF (SOURCE(J) .EQ. 0.0) THEN
	    REMEFF(J) = 1.0
	  ELSE
	    REMEFF(J) = RE(J)
	  ENDIF

25	CONTINUE

	QCIRC = 0.0
	QOATOT = 0.0

	DO 30 J = 1, NZONES

C
C  READ INPUT DATA; AIR FLOWS ARE CONVERTED FROM AIR CHANGES PER HOUR
C  TO M3/S (OR FT3/S)
C

	  QSUP(J) = XIN(J)*ZONEVOL(J) / 3600.0
	  QCIRC = QCIRC + QSUP(J)
	  QOA(J) = XIN(J+NZONES)*ZONEVOL(J) / 3600.0
	  QOATOT = QOATOT + QOA(J)

C
C  CALCULATE THE RETURN AIR TIME STEP DELAY
C

	  IF (QSUP(J) .EQ. 0.0) THEN
	    IDELAY(J) = 60
	  ELSE
	    IDELAY(J) = NINT( RETVOL(J) / (QSUP(J)*DT) )
	  ENDIF


	  IF (IDELAY(J) .GT. 60) IDELAY(J) = 60

30	CONTINUE

	QRETTOT = 0.0

	DO 35 J = 1, NZONES

C
C  CALCULATE THE RETURN AIR FLOW BY A FLOW BALANCE; IF THE VENTILATION 
C  SYSTEM IS OFF, INFILTRATION CONTINUES BUT INTER-ZONE FLOW IS ASSUMED 
C  TO BE ZERO 
C

	  QIZFNET = 0.0
	  QIZFOUT = 0.0

	  IF (QCIRC .EQ. 0.0) THEN
	    QRET(J) = INFILT(J)*ZONEVOL(J) / 3600.0
	  ELSE
	    DO 33 L = 1, NZONES
	      QIZFNET = QIZFNET + QIZF(L,J) - QIZF(J,L)
	      QIZFOUT = QIZFOUT + QIZF(J,L)
33	    CONTINUE
	    QRET(J) = QSUP(J) + QIZFNET + INFILT(J)*ZONEVOL(J) / 3600.0
	  ENDIF

	  QRETTOT = QRETTOT + QRET(J)

C
C  CALCULATE CONSTANTS FOR THE dC/dT FUNCTION F
C

	  A1(J) = SOURCE(J) / ZONEVOL(J) + ( INFILT(J) * COA ) / 3600.0
	  A2(J) = QSUP(J) / ZONEVOL(J)

	  IF (QCIRC .EQ. 0.0) THEN
	    A3(J) = QRET(J) * REMEFF(J) / ZONEVOL(J)
	  ELSE
	    A3(J) = ( QRET(J) * REMEFF(J) + QIZFOUT ) / ZONEVOL(J)
	  ENDIF

35	CONTINUE

C
C  CALL THE MAIN CALCULATION SUBROUTINE; IF THE VENTILATION SYSTEM IS OFF,
C  A SHORTENED CALCULATION WITHOUT RECIRCULATION IS PERFORMED
C

	IF (QCIRC .EQ. 0.0) THEN
	  CALL SHORTCALC (C,COA,DT,NSTEPS,NZONES)
	ELSE
	  CALL CALC (C,COA,DT,IDELAY,QRET,QRETTOT,QCIRC,QOATOT,
     .	             QIZF,REMEFF,NSTEPS,NZONES,ZONEVOL,FILTEFF)
	ENDIF


C
C  CALCULATE VALUES FOR OUTPUT; AIR FLOWS ARE CONVERTED BACK TO AIR
C  CHANGES PER HOUR
C

	QCRETTOT = 0.0

	DO 40 J = 1, NZONES

C
C  CALCULATE THE VOLUME FLOW RATE OF OUTSIDE AIR PER SOURCE UNIT;
C  ENGLISH UNITS IN FT3/MIN, SI UNITS IN L/S
C  

	  IF (SMULT(J) .LE. 0.0) THEN
	    QOAPU = 0.0
	  ELSEIF (UNITFLAG .EQ. 2) THEN
	    QOAPU = QOA(J)*60.0/SMULT(J)
	  ELSE
	    QOAPU = QOA(J)*1000.0/SMULT(J)
	  ENDIF 

	  QCRETTOT = QCRETTOT + QRET(J)*C(-IDELAY(J),J)*REMEFF(J)

	  OUT(J) = QOAPU
	  OUT(NZONES+J) = C(0,J)
	  OUT(2*NZONES+J) = C(0,J)*REMEFF(J)
40	CONTINUE

	CAHU = QCRETTOT/QRETTOT

	OUT(3*NZONES+1) = CAHU
	IF (QCIRC .NE. 0.0) THEN
	  OUT(3*NZONES+2) = ((QCIRC-QOATOT)*CAHU*(1-FILTEFF(1))+
     .	                    QOATOT*COA*(1-FILTEFF(2)))*
     .	                    (1-FILTEFF(3))/QCIRC
	ENDIF

	RETURN 1
	END



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  SUBROUTINE SHORTCALC
C
C  USES AN ITERATIVE MODIFIED EULER (SECOND-ORDER RUNGE-KUTTA) 
C  METHOD FOR SOLVING THE CONCENTRATION DIFFERENTIAL EQUATION WHEN
C  THE VENTILATION SYSTEM IS OFF
C
C  SUBROUTINE SHORTCALC IS CALLED BY:
C  	SUBROUTINE TYPE60
C
C  SUBROUTINE SHORTCALC CALLS:
C	FUNCTION F
C

	SUBROUTINE SHORTCALC (C,COA,DT,NSTEPS,NZONES)

	IMPLICIT NONE
	INTEGER I, J, NSTEPS, NZONES, ERRFLAG
	REAL F, C(-60:60,NZONES), COA, DT, C1, ERR


	DO 30 I = 0, NSTEPS - 1

C
C  ESTIMATE CONCENTRATION AT STEP I+1
C

	  DO 10 J = 1, NZONES
	    C(I+1,J) = C(I,J) + DT*F(0.0,C(I,J),0.0,J)
10	  CONTINUE

C
C  CALCULATE CONCENTRATION AT STEP I+1 USING THE SLOPE ESTIMATE AT STEP I+1/2
C

15	  ERRFLAG = 0
	  DO 20 J = 1, NZONES
 	    C1 = C(I,J) + DT*F(0.0,(C(I,J)+C(I+1,J))/2.0,0.0,J)
	    ERR = ABS( C1 - C(I+1,J) )
	    IF (ERR .GE. 1.0E-8) ERRFLAG = 1
	    C(I+1,J) = C1
20	  CONTINUE
	  IF (ERRFLAG .EQ. 1) GOTO 15

30	CONTINUE

C
C  RESET THE CURRENT CONCENTRATION TIME STEP TO 0 TO GET READY FOR THE NEXT
C  ENTRY INTO THIS SUBROUTINE
C

	DO 40 J = 1, NZONES
	  DO 40 I = -60+NSTEPS, NSTEPS
	    C(I-NSTEPS,J) = C(I,J)
40	CONTINUE

	RETURN
	END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  SUBROUTINE CALC
C  
C  USES AN ITERATIVE MODIFIED EULER (SECOND-ORDER RUNGE-KUTTA) 
C  METHOD FOR SOLVING THE CONCENTRATION DIFFERENTIAL EQUATION AND
C  INCLUDES A TIME DELAY IN THE HVAC AIR PATH FOR EACH ZONE.
C 
C  SUBROUTINE CALC IS CALLED BY:
C  	SUBROUTINE TYPE60
C
C  SUBROUTINE CALC CALLS:
C	FUNCTION F
C

	SUBROUTINE CALC (C,COA,DT,IDELAY,QRET,QRETTOT,QCIRC,QOATOT,        
     .	                 QIZF,REMEFF,NSTEPS,NZONES,ZONEVOL,FILTEFF)

	IMPLICIT NONE

	INTEGER I, J, L, NSTEPS, NZONES
	INTEGER IDELAY(NZONES), ERRFLAG

	REAL F, C(-60:60,NZONES), COA, DT, QRET(NZONES), QRETTOT, QCRETTOT
	REAL QCIRC, QOATOT, QIZF(25,NZONES), REMEFF(NZONES)
	REAL ZONEVOL(NZONES), FILTEFF(3)
	REAL QCIZF(25), CSUP, CAHU, C1, ERR

C
C  BEGIN THE MAIN CALCULATION LOOP BY CALCULATING A NEW CONCENTRATION
C  FOR EACH ROOM AND THEN PROCEEDING TO THE NEXT TIME STEP 
C

	DO 50 I = 0, NSTEPS-1

C
C  CALCULATE THE SUPPLY AIR CONCENTRATION AND THE INTER-ZONE 
C  CONTAMINANT FLOWS FOR TIME STEP I
C

	  QCRETTOT = 0.0

	  DO 10 J = 1, NZONES
	    QCRETTOT = QCRETTOT + QRET(J)*C(I-IDELAY(J),J)*REMEFF(J)
	    QCIZF(J) = 0.0

	    DO 5 L = 1, NZONES
	      QCIZF(J) = QCIZF(J) + QIZF(L,J)*C(I,L)/ZONEVOL(J)
5	    CONTINUE
10	  CONTINUE

	  CAHU = QCRETTOT / QRETTOT
	  CSUP = ((QCIRC-QOATOT)*CAHU*(1-FILTEFF(1))+
     .	            QOATOT*COA*(1-FILTEFF(2)))*(1-FILTEFF(3))/QCIRC

C
C  ESTIMATE CONCENTRATION AT STEP I+1
C


	  DO 20 J = 1, NZONES
	    C(I+1,J) = C(I,J) + DT*F(CSUP,C(I,J),QCIZF(J),J)
20	  CONTINUE

C
C  CALCULATE THE SUPPLY AIR CONCENTRATION AND THE INTER-ZONE 
C  CONTAMINANT FLOWS FOR TIME STEP I+1/2
C

22	  QCRETTOT = 0.0

	  DO 30 J = 1, NZONES
	    IF (IDELAY(J) .LE. 0 ) THEN
	      QCRETTOT = QCRETTOT + QRET(J)*REMEFF(J)*
     .	                                   (C(I,J)+C(I+1,J))/2.0
	    ELSE
	      QCRETTOT = QCRETTOT + QRET(J)*REMEFF(J)*
     .	               (C(I-IDELAY(J),J)+C(I-IDELAY(J)+1,J))/2.0
	    ENDIF

 	    QCIZF(J) = 0.0

	    DO 25 L = 1, NZONES
	      QCIZF(J) = QCIZF(J) + QIZF(L,J)*(C(I,L)+C(I+1,L))/
     .	                                           (2.0*ZONEVOL(J))
25	    CONTINUE
30	  CONTINUE

	  CAHU = QCRETTOT / QRETTOT
	  CSUP = ((QCIRC-QOATOT)*CAHU*(1-FILTEFF(1))+
     .	             QOATOT*COA*(1-FILTEFF(2)))*(1-FILTEFF(3))/QCIRC

C
C  CALCULATE CONCENTRATION AT STEP I+1 USING THE SLOPE ESTIMATE AT STEP I+1/2
C

	  ERRFLAG = 0
	  DO 40 J = 1, NZONES
	    C1 = C(I,J) + DT*F(CSUP,(C(I,J)+C(I+1,J))/2.0,QCIZF(J),J)
	    ERR = ABS( C1 - C(I+1,J) )
	    IF (ERR .GE. 1.0E-8) ERRFLAG = 1
	    C(I+1,J) = C1
40	  CONTINUE
	  IF (ERRFLAG .EQ. 1) GOTO 22

50	CONTINUE

C
C  RESET THE CURRENT CONCENTRATION TIME STEP TO 0 TO GET READY FOR THE NEXT
C  ENTRY INTO THIS SUBROUTINE
C

	DO 60 J = 1, NZONES
	  DO 60 I = -60+NSTEPS, NSTEPS
	    C(I-NSTEPS,J) = C(I,J)
60	CONTINUE

	RETURN
	END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C	 
C  FUNCTION F	 
C 
C  THIS FUNCTION IS FOR A ZONE WITH A SOURCE TERM, A VENTILATION 
C  AIR TERM, AN INTER-ZONE AIR FLOW TERM, AND A REMOVAL TERM
C 

	FUNCTION F (CSUP, C, QCIZF, J)

	IMPLICIT NONE

	INTEGER J, NZONES
	REAL F, A1(25), A2(25), A3(25), CSUP, C, QCIZF

	COMMON A1, A2, A3

	F = A1(J) + A2(J)*CSUP + QCIZF - A3(J)*C 

	RETURN
	END