
		!============================================================================================
		PROGRAM DESIGN
		
		!THIS WILL BE USED TO DESIGN A NETWORK HAVING IN MIND THE INLET OF ALL THE STREAMS.
		
		USE Prop_mod
		!INCLUDE "mkl_lapack95.f90"
		USE blas95
		USE f95_precision
		USE lapack95	

		IMPLICIT NONE	

		!DECLARATIONS OF VARIABLES
		INTEGER :: j, k, x, y, ij, ik, ii, info
		INTEGER :: HTMIX, CDMIX, HEX, HSID, CSID, num, fhid, fcid
		INTEGER :: read_status, cbranch, hbranch
		
		REAL :: hdata, cdata, names, streams, ht1, ht2, mh, ct1, ct2, mc, cpc, cph, mhb, mcb
		REAL :: cbmass, hbmass, bmass, cbypass, hbypass, mcbypass, mhbypass
		REAL :: t, rhom, vism, tcm, hcm
		REAL :: X1, X2, X3, X4, X5   !DUMMY VARIABLES
		
		REAL, DIMENSION(:), ALLOCATABLE :: RefTh, Rhflow1, Rhhtc1, Rhflow2, Rhhtc2, RefTc, Rcflow1, Rchtc1, Rcflow2, Rchtc2
		REAL, DIMENSION(:), ALLOCATABLE :: area, tubeod, tubeid,  kex, kpipe, ffc, ffh
		REAL, DIMENSION(:), ALLOCATABLE :: htmflow, cdmflow, htemp, ctemp, temp
		REAL, DIMENSION(:), ALLOCATABLE :: cprref, hprref
		REAL, DIMENSION(:), ALLOCATABLE :: aexc, aemxc, aexh, aemxh
		REAL, DIMENSION(:), ALLOCATABLE :: hbm, cbm
		
		REAL, DIMENSION(:,:), ALLOCATABLE :: xtemp, nxtemp
		
		INTEGER, DIMENSION(:), ALLOCATABLE :: tubes, passes, coldside, insert, shells
		INTEGER, DIMENSION(:), ALLOCATABLE :: htp, ctp, ztemp, ipiv

		LOGICAL :: branchc, branchh, bypassc, bypassh
		
		CHARACTER :: STID*2, dataname*8, stname*20, branching*1
		CHARACTER*10, DIMENSION(:), ALLOCATABLE :: CID, HID
		
	! FILES
		OPEN(UNIT = 21, FILE = 'HexData.csv', STATUS = 'old', ACTION = 'read')	! Input Data
		OPEN(UNIT = 11, FILE = 'InputData.csv', STATUS = 'old', ACTION = 'read')	! Input Data	
		OPEN(UNIT = 30, FILE = 'MATRIX.csv', STATUS = 'unknown')
		OPEN(UNIT = 31, FILE = 'MATRIXinv.csv', STATUS = 'unknown')
		
	!INPUT DATA			
		CALL Props(hrhoa, hrhob, hvisa, hvisb, htca, htcb, hhcapa, hhcapb, crhoa, crhob, cvisa, cvisb, ctca, ctcb, chcapa, chcapb, cst, hst, cname, hname)

		WRITE(*, 501) 'ENTER THE NUMBER OF MIXING POINTS FOR ...'
		WRITE(*,501, ADVANCE = 'no') '... HOT STREAM  ' 
		READ(*, *) HTMIX
		WRITE(*, *)	
		
		WRITE(*,501, ADVANCE = 'no') '... COLD STREAM  ' 
		READ(*, *) CDMIX
		WRITE(*, *)

		WRITE(*,501, ADVANCE = 'no') 'ENTER THE NUMBER OF HEAT EXCHANGERS  ' 
		READ(*, *) HEX
		WRITE(*, *) 
		
		num = hst + cst + 2 * HEX + HTMIX + CDMIX
		ALLOCATE(temp(num))
		ALLOCATE(ztemp(num))
		ALLOCATE(xtemp(num, num))
		
	!STREAM DATA
		ALLOCATE(htmflow(hst))
		ALLOCATE(htemp(hst))
		ALLOCATE(HID(hst))
		WRITE(*, *)
		WRITE(*, *) 'ENTER THE INPUT DATA FOR THE HOT STREAMS'
		WRITE(*, *)
		
		DO i = 1, hst
			WRITE(STID,'(I2)') i
			WRITE(*,*) ADJUSTL('HOT STREAM ')//ADJUSTL(STID),' - ' , hname(i)
			WRITE(*, 501, ADVANCE = 'no') 'INLET FLOWRATE (kg/s) =  '
			READ(*, *) hdata
			htmflow(i) = hdata
			WRITE(*, 501, ADVANCE = 'no') 'INLET TEMPERATURE (C) =  '
			READ(*, *) hdata
			htemp(i) = hdata
			HID(i) = ADJUSTL('STREAM H')//ADJUSTL(STID) 
		ENDDO 
		
		ALLOCATE(cdmflow(cst))
		ALLOCATE(ctemp(cst))
		ALLOCATE(CID(cst))
		WRITE(*, *)
		WRITE(*, *) 'ENTER THE INPUT DATA FOR THE COLD STREAMS'
		WRITE(*, *)
		
		DO i = 1, cst
			WRITE(STID,'(I2)') i
			WRITE(*, *) ADJUSTL('COLD STREAM ')//ADJUSTL(STID), ' - ', cname(i)
			WRITE(*, 501, ADVANCE = 'no') 'INLET FLOWRATE (kg/s) =  '
			READ(*, *) cdata
			cdmflow(i) = cdata
			WRITE(*, 501, ADVANCE = 'no') 'INLET TEMPERATURE (C) =  '
			READ(*, *) cdata
			ctemp(i) = cdata
			CID(i) = ADJUSTL('STREAM C')//ADJUSTL(STID) 
		ENDDO
		
	!HEAT EXCHANGER DATA
		WRITE(*, *)
		WRITE(*, *) 'HEAT EXCHANGER INFORMATION FOR EACH STREAM'
		WRITE(*, *)
		DO i = 1, 7
			READ(21, *) stname
		ENDDO

		ALLOCATE(shells(HEX))
		ALLOCATE(RefTh(HEX))
		ALLOCATE(Rhflow1(HEX))
		ALLOCATE(Rhhtc1(HEX))
		ALLOCATE(Rhflow2(HEX))
		ALLOCATE(Rhhtc2(HEX))
		ALLOCATE(RefTc(HEX))
		ALLOCATE(Rcflow1(HEX))
		ALLOCATE(Rchtc1(HEX))
		ALLOCATE(Rcflow2(HEX))
		ALLOCATE(Rchtc2(HEX))
		ALLOCATE(area(HEX))
		ALLOCATE(tubeid(HEX))
		ALLOCATE(tubeod(HEX))
		ALLOCATE(tubes(HEX))
		ALLOCATE(passes(HEX))
		ALLOCATE(coldside(HEX))
		ALLOCATE(insert(HEX))
		ALLOCATE(kex(HEX))
		ALLOCATE(kpipe(HEX))
		ALLOCATE(ffc(HEX))
		ALLOCATE(ffh(HEX))
		
	!HOT STREAM DATA
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, shells(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, RefTh(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rhflow1(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rhhtc1(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rhflow2(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rhhtc2(1:HEX)
		
	!COLD STREAM DATA
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, RefTc(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rcflow1(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rchtc1(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rcflow2(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, Rchtc2(1:HEX)

	!HEAT EXCHANER SPECIFICATIONS
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, area(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, tubeod(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, tubeid(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, tubes(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, passes(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, coldside(1:HEX)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, insert(1:HEX)
		READ(21, *)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, kex(1:HEX) 
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, kpipe(1:HEX) 
		READ(21, *)
		READ(21, *)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, ffh(1:HEX) 
		READ(21, *)
		READ(21, *)
		READ(21, *, ERR = 901, IOSTAT = read_status) stname, ffc(1:HEX) 
	
		WRITE(*, *)
		WRITE(*, *) 'CALCULATIONS FOR EACH HEAT EXCHANGER'
	! SPECIFY THE ID FOR THE EXCHANGER 1, 2, .... BASED ON HOW THEY WERE ENTERED INTO THE PROGRAM
	! THE HOT STREAM ID AND COLD STREAM NEEDS TO BE SPECIFIED.
	! CALCULATIONS WILL THEN BE DONE FOR EACH 

		ALLOCATE(hprref(hst))
		DO i = 1, hst
			t = RefTh(i)
			rhoa = hrhoa(i)
			rhob = hrhob(i)
			visa = hvisa(i)
			visb = hvisb(i)
			tca = htca(i)
			tcb = htcb(i)
			hca = hhcapa(i)
			hcb = hhcapb(i)
			CALL Property(t, rhoa, rhob, visa, visb, tca, tcb, hca, hcb, rhom, vism, tcm, hcm)		
			hprref(i) = hcm * vism / tcm
		ENDDO

		ALLOCATE(cprref(cst))
		DO i = 1, cst 
			t = RefTc(i)
			rhoa = crhoa(i)
			rhob = crhob(i)
			visa = cvisa(i)
			visb = cvisb(i)
			tca = ctca(i)
			tcb = ctcb(i)
			hca = chcapa(i)
			hcb = chcapb(i)
			CALL Property(t, rhoa, rhob, visa, visb, tca, tcb, hca, hcb, rhom, vism, tcm, hcm)
			cprref(i) = hcm * vism / tcm 
		ENDDO

		ALLOCATE(aexh(HEX))
		ALLOCATE(aemxh(HEX))
		CALL ExFIT(HEX, Rhhtc1, Rhhtc2, Rhflow1, Rhflow2, aexh, aemxh)

		ALLOCATE(aexc(HEX))
		ALLOCATE(aemxc(HEX))
		CALL ExFIT(HEX, Rchtc1, Rchtc2, Rcflow1, Rcflow2, aexc, aemxc)

	!===========================================================================================
	!HEAT EXCHANGER ANALYSIS
	!===========================================================================================
		WRITE(*, *)
		WRITE(*, *) 'STARTING CALCULATION ...'
		
		branchh = .TRUE.
		branchc = .TRUE. 
		bypassc = .FALSE.
		bypassh = .FALSE.
		hbmass = 0.0
		cbmass = 0.0 			
		
		x = 1
		y = 2 * HEX + 1 	
		xtemp = 0 
		ztemp = 0
		ij = 1 
		fhid = 0 
		fcid = 0
		
		DO i = 1, HEX
			WRITE(*, *)
			WRITE(*, *)
			WRITE(*, *) 'FOR HEAT EXCHANGER ', i
			WRITE(*, *)
		
		! HOT STREAM
			WRITE(*, 501, ADVANCE = 'no') 'ENTER THE HOT STREAM ID '
			READ(*, *) HSID
			ht1 = htemp(HSID)
			mh = htmflow(HSID)			
 
			IF (branchh == .TRUE.) THEN
				WRITE(*, 501, ADVANCE = 'no') 'IS THERE BRANCHING? Y or N '
				READ(*, *) branching
				IF (branching .EQ. 'Y' .OR. branching .EQ. 'y') THEN
					WRITE(*, 502, ADVANCE = 'no') 'ENTER NUMBER OF BRANCHES FOR HOT ID ', HSID
					READ(*, *) hbranch
					ALLOCATE(htp(hbranch))
					ALLOCATE(hbm(hbranch))
					ELSEIF (branching .EQ. 'N' .OR. branching .EQ. 'n') THEN 
					hbranch = 0
				ENDIF
				j = 1
			ENDIF
			
			IF (hbranch .EQ. 0) THEN
				bmass = 100
				branchh = .TRUE.
				j = 1
			ELSEIF (j .LT. hbranch) THEN 
				WRITE(*, 502, ADVANCE = 'no') '% OF MASS FLOW IN BRANCH ', j
				READ(*, *) bmass
				hbm(j) = bmass / 100
				hbmass = hbmass + bmass
				branchh = .FALSE.
			ELSE 
				bmass = 100 - hbmass
				hbm(j) = bmass / 100
				branchh = .TRUE.
				hbmass = 0.0
			ENDIF
			mhb = (bmass /100) * mh
			WRITE(*, 501, ADVANCE = 'no') 'BYPASS FRACTION (BTN 0 AND 1) '
			READ(*, *) hbypass
			IF (hbypass .GT. 0.0) THEN 
				mhbypass = hbypass * mhb
				mhb = (1 - hbypass) * mhb
				hbm(j) = hbypass * (bmass / 100)
				bypassh = .TRUE.
			ENDIF
			
		! COLD STREAM		
			WRITE(*, *)
			WRITE(*, 501, ADVANCE = 'no') 'ENTER THE COLD STREAM ID '
			READ(*, *) CSID
			ct1 = ctemp(CSID)
			mc = cdmflow(CSID)

			IF (branchc .EQ. .TRUE.) THEN
				WRITE(*, 501, ADVANCE = 'no') 'IS THERE BRANCHING? Y or N '
				READ(*, *) branching
				IF (branching .EQ. 'Y' .OR. branching .EQ. 'y') THEN 
					WRITE(*, 502, ADVANCE = 'no') 'ENTER NUMBER OF BRANCHES FOR FOR COLD ID ', CSID
					READ(*, *) cbranch
					ALLOCATE(ctp(cbranch))
					ALLOCATE(cbm(cbranch))
				ELSEIF (branching .EQ. 'n' .OR. branching .EQ. 'N') THEN 
					cbranch = 0
				ENDIF
				k = 1
			ENDIF 
			
			IF (cbranch .EQ. 0 ) THEN
				bmass = 100 
				branchc = .TRUE.
			ELSEIF (k .LT. cbranch) THEN 
				WRITE(*, 502, ADVANCE = 'no') '% OF MASS FLOW IN BRANCH ', k
				READ(*, *) bmass
				cbm(k) = bmass / 100
				cbmass = cbmass + bmass
				branchc = .FALSE.
			ELSE 
				bmass = 100 - cbmass
				cbm(k) = bmass / 100
				branchc = .TRUE.
				cbmass = 0.0
			ENDIF
			mcb = (bmass /100) * mc 
			WRITE(*, 501, ADVANCE = 'no') 'BYPASS FRACTION (BTN 0 AND 1) '
			READ(*, *) cbypass
			IF (cbypass .GT. 0.0) THEN 
				mcbypass = cbypass * mcb
				mcb = (1 - cbypass) * mcb
				cbm(k) = cbypass * (bmass / 100) 
				bypassc = .TRUE.
			ENDIF
			
			IF (i .GT. 1) THEN 
				CALL INPUT1(y, CSID, fcid, k)
				CALL INPUT1(x, HSID, fhid, j)
			ENDIF

		!	WRITE(*, *)
		!	PRINT *, 'INLET TEMPERATURES: COLD = ', y, ' = ', ct1, 'HOT = ', x, ' = ', ht1 
			
			CALL EXUNIT(i, HSID, CSID, X1, X2)
		
		!	WRITE(*, *)
		!	PRINT *, 'OUTLET PARAMETERS: X1 = ', X1, 'X2 = ', X2
		!	PAUSE
			
			X3 = 1 - X1
			X4 = 1 - X2
			
		!ARRAY OF COEFFICIENTS
			xtemp(ij, x) = X1
			xtemp(ij, y) = X3
			xtemp(ij, y + k) = -1

			ij = ij + 1
			xtemp(ij, x) = X4
			xtemp(ij, x + j) = -1
			xtemp(ij, y) = X2
			
			IF (HSID .NE. fhid) THEN 
				ij = ij + 1
				ztemp(ij) = htemp(HSID)
				xtemp(ij, x) = 1 !check this again
			ENDIF
			
			IF (CSID .NE. fcid) THEN 
				ij = ij + 1
				ztemp(ij) = ctemp(CSID)
				xtemp(ij, y) = 1 !check this again
			ENDIF

			IF (cbypass .GT. 0) THEN
				ij = ij + 1
				xtemp(ij, y) = cbypass
				xtemp(ij, y + 1) = (1 - cbypass)
				xtemp(ij, y + 2) = -1
			ENDIF
			
			IF (hbypass .GT. 0) THEN
				ij = ij + 1 
				xtemp(ij, x) = hbypass
				xtemp(ij, x + 1) = (1 - hbypass)
				xtemp(ij, x + 2) = -1
			ENDIF
		
			IF (cbypass .GT. 0 .AND. cbranch .GT. 0) THEN 
				ctp(k) = y + 2 * k
			ELSEIF (cbypass .GT. 0) THEN 
				ctp(k) = y + 2
			ELSEIF (bypassc .EQ. .TRUE.) THEN 
				ctp(k) = y + k + 1
				bypassc = .FALSE.
			ELSEIF (cbranch .GT. 0) THEN 
				ctp(k) = y + k
			ENDIF
				
			IF (hbypass .GT. 0 .AND. hbranch .GT. 0) THEN 
				htp(j) = x + 2 * j
			ELSEIF (hbypass .GT. 0) THEN 
				htp(j) = x + 2
			ELSEIF (bypassh .EQ. .TRUE.) THEN 
				htp(j) = x + j + 1
				bypassh = .FALSE.
			ELSEIF (hbranch .GT. 0) THEN 
				htp(j) = x + j
			ENDIF	

			
			IF (j .EQ. hbranch) THEN
				ij = ij + 1 
				DO ii = 1, hbranch
					xtemp(ij, htp(ii)) = hbm(ii)
				ENDDO
				x = htp(ii - 1) + 1
				xtemp(ij, x) = -1
				DEALLOCATE(htp)
				DEALLOCATE(hbm)
			ENDIF

			IF (k .EQ. cbranch) THEN
				ij = ij + 1 
				DO ii = 1, cbranch
					xtemp(ij, ctp(ii)) = cbm(ii)
				ENDDO
				y = ctp(ii - 1) + 1
				xtemp(ij, y) = -1
				DEALLOCATE(ctp)
				DEALLOCATE(cbm)
			ENDIF

			
			IF (k .LT. cbranch) k = k + 1
			IF (j .LT. hbranch) j = j + 1
			
			IF (HSID .GT. fhid) THEN 
				fhid = HSID
			ENDIF
		
			IF (CSID .GT. fcid) THEN 
				fcid = CSID
			ENDIF
			
			ij = ij + 1
			IF (i .EQ. 1) THEN 
				y = y + k
				x = x + j 
			ENDIF 
			
						
			WRITE(*, *)

			
		ENDDO
		
		ALLOCATE(nxtemp(num, 1))
		PRINT *, ztemp
		WRITE(*, *)
		
		PRINT *, xtemp
		WRITE(*, *)
		
		DO i = 1, num 
			WRITE(30, '(*(G0.4, ",", :))') xtemp(i, :)
		ENDDO
		
		ALLOCATE(ipiv(num))
		nxtemp = RESHAPE(ztemp, (/num, 1/))

		CALL gesv(xtemp, nxtemp, ipiv, info)
		PRINT *, 'info = ', info, 'if = 0 then good  else problem'
		WRITE(*, *)

		!CALL getrf(xtemp, ipiv)	
		
		PRINT *, xtemp
		PAUSE
		
		!CALL getrs(xtemp, ipiv, nxtemp)
		!CALL sgetrs('N', num, 0, xtemp, num, ipiv, nxtemp, num, info)   

		WRITE(30, *)
		WRITE(30, *)
	
		
		DO i = 1, num 
			WRITE(31, '(*(G0.4, ",", :))') xtemp(i, :)
		ENDDO		

		
		WRITE(*, *)
		WRITE(*,*) ztemp
		WRITE(*, *)
		
		
		WRITE(*, *) nxtemp

		WRITE(*, *)
		WRITE(*, *)			'WELL DONE'
		
		
501		FORMAT(A50)	
502		FORMAT(A48, I2)	
	
	!HANDLING ERRORS
901		IF (read_status > 0) THEN 
			WRITE(*,*)
			WRITE(*,*) 'Data COULD NOT be read from file', read_status
			STOP
		END IF
		
		CONTAINS
			
			SUBROUTINE Property(t, rdata, rdatb, vdata, vdatb, tdata, tdatb, hdata, hdatb, rhom, vism, tcm, hcm)
				REAL, INTENT(IN) :: rdata, rdatb, vdata, vdatb, tdata, tdatb, hdata, hdatb, t
				REAL, INTENT(OUT) :: rhom, vism, tcm, hcm
				REAL :: tin
				
				tin = t + 273.15
				vism = vdata * Exp(vdatb / tin)
				tcm = tdata + tdatb * tin
				hcm = hdata + hdatb * tin			
				CALL Den(tin, rdata, rdatb, rhom)
			END SUBROUTINE Property
			
			SUBROUTINE Den(tin, sumx, sumy, rhom)
				REAL, INTENT(IN) :: sumx, sumy, tin
				REAL, INTENT(OUT) :: rhom
				
				rhom = sumx + sumy * tin
			END SUBROUTINE Den
			
			SUBROUTINE ExFIT(n, ha, hb, fra, frb, aex, aemx)
				INTEGER :: n
				REAL, PARAMETER :: eps = 0.001
				REAL, DIMENSION(n) :: test1, fra, frb, ha, hb, aex, aemx
				
				test1 = Abs(fra - frb)
				WHERE (test1 < 0.01)
					aex = ha
					aemx = 0.0
				ELSEWHERE (hb < eps)
					aex = ha
					aemx = 0.0
				ELSEWHERE
					aemx = Log(hb) - Log(ha) 
					aex = Log(fra) - Log(frb)
					aemx = aemx /aex
					aex = ha / (fra ** aemx)
				END WHERE				
			END SUBROUTINE ExFIT
			
			SUBROUTINE EXUNIT(n, nh, nc, x1, x2)
				INTEGER :: n, nh, nc, nshell
				REAL :: flow1, flow2, A, U
				REAL, INTENT(OUT) :: x1, x2
				
				flow1 = mhb
				flow2 = mcb
				nshell = shells(n)
				A = area(n)
				
				CALL OVLHTC(n, nh, nc, flow1, flow2, U, Cpc, Cph)
				
				CALL EXTEMP(nshell, flow1, flow2, Cpc, Cph, x1, x2, A, U)		
			
			END SUBROUTINE EXUNIT
			
			SUBROUTINE OVLHTC(n, ij, ik, flow1, flow2, U, Cpc, Cph)
				INTEGER :: n, ik, ij, side
				REAL :: x1, x2, x3, x4, x5, x6, tid, tod, Univ, trw
				REAL :: hpr, hh, flow1, hfoul
				REAL :: cpr, hc, flow2, cfoul
				REAL, INTENT(OUT) :: U, Cpc, Cph

				tod = tubeod(n)
				tid = tubeid(n)
				cfoul = ffc(n)
				hfoul = ffh(n)
				side = coldside(n)

			!HOT STREAM PROPERTIES
				x1 = aexh(ij)
				x2 = aemxh(ij)
				x3 = Rhflow1(ij)
				x4 = hprref(ij)
				x5 = x1 * (flow1 ** x2)

				rhoa = hrhoa(ij)
				rhob = hrhob(ij)
				visa = hvisa(ij)
				visb = hvisb(ij)
				tca = htca(ij)
				tcb = htcb(ij)
				hca = hhcapa(ij)
				hcb = hhcapb(ij)
				
				CALL Property(ht1, rhoa, rhob, visa, visb, tca, tcb, hca, hcb, rhom, vism, tcm, hcm)
				hpr = hcm * vism / tcm
				Cph = hcm
				hh = x5 * (hpr / x4) ** (1.0/3.0)
				
			!COLD STREAM PROPERTIES
				x1 = aexc(ik)
				x2 = aemxc(ik)
				x3 = Rcflow1(ik)
				x4 = cprref(ik)
				x5 = x1 * (flow2 ** x2)
				
				rhoa = crhoa(ik)
				rhob = crhob(ik)
				visa = cvisa(ik)
				visb = cvisb(ik)
				tca = ctca(ik)
				tcb = ctcb(ik)
				hca = chcapa(ik)
				hcb = chcapb(ik)
				
				CALL Property(ct1, rhoa, rhob, visa, visb, tca, tcb, hca, hcb, rhom, vism, tcm, hcm)
				cpr = hcm * vism / tcm
				Cpc = hcm
				hc = x5 * (cpr / x4) ** (1.0/3.0)
				
				trw = 0.5 * 0.001 * (tod - tid) / 48.0
				
				IF (side .GT. 0) THEN
					Univ = 1 / hc + cfoul + (tod / tid) * (1 / hh + hfoul) + trw
				ELSE
					Univ = 1 / hh + hfoul + (tod / tid) * (1 / hc + cfoul) + trw
				ENDIF
    
				U = 1 / Univ			
			END SUBROUTINE OVLHTC
			
			SUBROUTINE EXTEMP(nshell, flow1, flow2, Cpc, Cph, x1, x2, A, U)
				INTEGER :: idc, nshell
				REAL :: U, A, R, CPmin, Cpc, Cph, NTU, eff
				REAL, INTENT(IN) :: flow1, flow2
				REAL, INTENT(OUT) :: x1, x2
				
				Cph = flow1 * Cph
				Cpc = flow2 * Cpc
				
				R = Cph / Cpc
				CPmin = Cph
				idc = 1 
				IF (R .GT. 1) THEN 
					R = 1 / R
					CPmin = Cpc
					idc = 2 
				ENDIF 
				
				NTU = U * A / CPmin
				x1 = (1 + R * R) ** (0.5)
				x2 = -NTU * x1
				eff = (1 + R + x1 * ((1 + Exp(x2)) / (1 - Exp(x2))))
				eff = 2 / eff
				
				IF (nshell .GT. 1) THEN 
					x1 = ((1 - eff * R) / (1 - eff)) ** nshell
					eff = (x1 - 1) / (x1 - R)
				ENDIF
		
				IF (idc .EQ. 1) THEN
					x1 = eff * R
					x2 = eff
				ELSE
					x1 = eff
					x2 = eff * R
			   ENDIF 
			END SUBROUTINE EXTEMP
			
			SUBROUTINE INPUT1(Z1, Z2, Z3, Z4)
				INTEGER, INTENT(OUT) :: Z1
				INTEGER, INTENT(IN) :: Z2, Z3, Z4
				
				IF (Z2 .EQ. Z3) THEN 
					IF (Z4 .LT. 1) THEN
						Z1 = Z1 + 1
					ENDIF
				ELSE
					Z1 = Z1 + 1
				ENDIF
			END SUBROUTINE INPUT1
			
		END PROGRAM DESIGN