IMPLICIT REAL*8 (A-H,O-Z)
        character*20 mess
        CHARACTER*20 PW
	INTEGER flood,fcond,ftest,team
	COMMON fv,tfeed,dv,tcond,treb,reflux,conc,prof,pa,pc
	REAL*8 DES(200,6),Y(200),P(200)
	REAL*8 XL(6),XU(6)
	DATA XL/75.0,70.0,2.0,100.0,250.0,10.0/
	DATA XU/125.0,300.0,20.0,300.0,500.0,100.0/
	data nmax/999/
	ftest=0
        write(6,10)
 10	format('data')
	read(*,*) team,fv,tfeed,dv,tcond,treb,reflux
	if(team.eq.1) OPEN(UNIT=7,FILE='team1.data')
        if(team.eq.2) open(unit=7,file='team2.data')
        if(team.eq.3) open(unit=7,file='team3.data')
        if(team.eq.4) open(unit=7,file='team4.data')
        if(team.eq.5) open(unit=7,file='team5.data')
        if(team.eq.6) open(unit=7,file='team6.data')
        if(team.eq.7) open(unit=7,file='team7.data')
        if(team.eq.8) open(unit=7,file='team8.data')
        if(team.eq.9) open(unit=7,file='team9.data')
        READ(7,20) PW
 20	FORMAT(A20)
	READ(7,30) NPREV,PBAL
 30	FORMAT(I3,2X,F10.2)
	DO 120 I = 1,NPREV
	   READ(7,110) (DES(I,J),J=1,6),Y(I),P(I)
 110	   FORMAT(F10.2,6F9.4,F10.2)
 120	CONTINUE
	NPREV = NPREV + 1
	IF(fv .LT. XL(1)) fv = XL(1)
	IF(fv .GT. XU(1)) fv = XU(1)
	IF(tfeed .LT. XL(2)) tfeed = XL(2)
	IF(tfeed .GT. XU(2)) tfeed = XU(2)
	IF(dv .LT. XL(3)) dv = XL(3)
	IF(dv .GT. XU(3)) dv = XU(3)
	IF(tcond .LT. XL(4)) tcond = XL(4)
	IF(tcond .GT. XU(4)) tcond = XU(4)
	IF(treb .LT. XL(5)) treb = XL(5)
	IF(treb .GT. XU(5)) treb = XU(5)
	IF(reflux .LT. XL(6)) reflux = XL(6)
	IF(reflux .GT. XU(6)) reflux = XU(6)
	DES(NPREV,1) = fv
	DES(NPREV,2) = tfeed
	DES(NPREV,3) = dv
	DES(NPREV,4) = tcond
	DES(NPREV,5) = treb
	DES(NPREV,6) = reflux
	call column(nprev,flood,fcond,ftest)
	IF(conc.GT.0.999) conc = 0.999
	IF(conc.LT.0.0) conc = 0.0
	Y(NPREV) = conc
	P(NPREV) = prof
	PBAL = PBAL + prof
	if(flood .eq. 0 .and. fcond .eq. 0) mess = 'OK'
	if(flood .eq. 1) mess = 'column flooded'
	if(fcond .eq. 1) mess = 'no condensate'
        write(6,130) fv
 130	format('Actual Feed Rate:',5x,f6.1)
        write(6,140) tfeed
 140	format('Actual Feed Temp:',5x,f6.1)
        write(6,150) dv
 150	format('Actual Dist Rate:',5x,f6.1)
        write(6,160) tcond
 160	format('Actual Cond Temp:',5x,f6.1)
        write(6,170) treb
 170	format('Actual Reboil Rate:',3x,f6.1)
        write(6,180) reflux
 180	format('Actual Reflux Rate:',3x,f6.1)
        write(6,190) conc
 190	format(///'Concentration:',5x,f8.2)
        write(6,200) prof
 200	format('Profit:',13x,f8.2)
        write(6,210) pa
 210	format('Column Pressure:',3x,f6.1)
        write(6,220) pc
 220	format('Cond Pressure:',5x,f6.1)
        write(6,230) mess
 230	format('Message:',12x,a20)
        write(6,240) nprev
 240	format(//'Total Runs to Date:',5x,i4)
        write(6,250) pbal
 250	format('Current Balance:',8x,f10.2)
        write(6,260)
 260	format('done')
 	REWIND 7
        write(7,20) PW
	WRITE(7,30) NPREV,PBAL
	DO 310 I = 1,NPREV
	   WRITE(7,110) (DES(I,J),J=1,6),Y(I),P(I)
 310	CONTINUE
	STOP
	END
C
	subroutine column(ns,flood,fcond,ftest)
	implicit real*8 (a-h,o-z)
	real*8 cmb(50,50),cinv(50,50)
	real*8 x(50),y(50),tl(50),tv(50),al(50),v(50),a(50),e(50),val(12)
	real*8 ka,ltop,lbot
	integer fstage,fstop,seed,flood,fcond,ftest
	common fv,tfeed,dv,tcond,treb,reflux,conc,prof,pa,pc
	if(ftest .eq. 1) open(unit=8,file='column.out')
	data fstage/20/
	data nstage/41/
	data ep/.0001/
	data tref/70.0/
	data cwater/0.04/
	data csteam/0.08/
	data value1/0.50/
	data value2/1.00/
	data value3/1.25/
	data feedc/0.015/
	data time/24.0/
	data ld/50/
	data seed/-9999/ldr/12/nrand/12/
	data qrmax/1300000.00/qcmax/1200000.00/
	seed = seed + ns
	call random(seed,nrand,ldr,val)
	sum = -6.0
	DO 165 j = 1,12
	   sum = sum + val(j)
 165	continue
	rerr = .002*sum
	xvfeed = .06 + rerr
	feed = fv*(0.143*xvfeed + 0.521*(1.0-xvfeed))
	xfeed = 0.143*xvfeed/(.143*xvfeed + 0.521*(1.0 - xvfeed))
	dinit = 0.143*dv
	wv = fv - dv
	winit = 0.521*wv
	xw = 0
	call press(xw,yw,treb,p)
	pa = p/760.0
	if(pa .gt. 6.0) then
	   conc = 0.0
	   prof = -50000.00
	   return
	endif
c
c
c                        initalize column
c                       for Thiele-Geddes
c
c
	xc = 1.0
	call bubble(tc,xc,p,treb)
	call bubble(tbf,xfeed,p,treb)
	call vol1(ka,tfeed,p)
	yfeed = ka*xfeed
	call dew(tdf,yfeed,p,treb)
	vtop = dinit*(reflux+1.0)
	vmin = dinit
	ltop = dinit*reflux
	if(tfeed .lt. tdf) then
	   lbot = ltop + feed
	   vbot = vtop
	else
	   z = 1.0
	   vbot = vtop - feed
	   if(vbot .lt. vmin) then
	      z = (vtop - vmin)/feed
	      vbot = vmin
	   endif
	   lbot = ltop + (1.0 - z)*feed
	endif
	ttop = tcond
	call cpl(cpc,xc)
	hd0 = dinit*cpc*(tcond - tref)
	call heat(hvap,xc)
	qc0 = dinit*(reflux+1.0)*hvap
	call cpl(cpr,xfeed)
	hb0 = al(1)*cpr*(treb - tref)
	call cpl(cpf,xfeed)
	hf0 = feed*cpf*(tfeed - tref)
	if(tfeed .ge. tdf) then 
	   call heat(hvap,xfeed)
	   hf0 = hf0 + feed*hvap
	endif
	qr0 = hd0 + hb0 + qc0 - hf0
	call heat(hvap,xfeed)
	vtemp = qr0/hvap
	if(vtemp .lt. vmin) then
	   qr0 = vmin*hvap
	   ttop = ttop+(qr0+hf0-hd0-hb0-qc0)/((reflux+1.0)*dinit*cpc)
	endif
	if(ftest .eq. 1) write(8,2) ttop
 2	format(' ttop:',f16.6)
	flood = 0
	if(ttop .gt. treb) then
	   flood = 1
	   conc = 0
	   pc = 0
	   prof = -25000.00
	   return
	endif
	del = (treb - ttop)/(nstage-1)
	t = treb + del
	do 10 i = 1,fstage-1
	   t = t - del
	   tl(i) = t
	   tv(i) = t
	   al(i) = lbot
	   v(i) = vbot
	   y(i) = xfeed
	   call eff(e(i),v(i),y(i),tl(i),p,al(i))
 10	continue
	al(1) = feed - dinit
	t = t - del
	tl(fstage) = t
	tv(fstage) = t
	al(fstage) = lbot
	if(tfeed .lt. tdf) then
	   v(fstage) = vbot
	else
	   v(fstage) = vtop
	endif
	y(fstage) = xfeed
	call eff(e(fstage),v(fstage),y(fstage),tfeed,p,al(fstage))
	nbegin = fstage + 1
	ntop = nstage - 1
	do 20 i = nbegin,nstage
	   t = t - del
	   tl(i) = t
	   tv(i) = t
	   al(i) = ltop
	   if(i .eq. nstage .and. tfeed .ge. tdf) al(i) = reflux*dinit
	   v(i) = vtop
	   if(i .eq. ntop .and. tfeed .ge. tdf) v(i)=(reflux+1.0)*dinit
	   y(i) = xfeed
	   call eff(e(i),v(i),y(i),tl(i),p,al(i))
 20	continue
	v(nstage) = dinit
	if(ftest .eq. 1) then 
	   write(8,21) tl(nstage)
 21	   format(' tl',f16.6)
	   do 25 i = 1,nstage
	      write(8,22) i,al(i),v(i)
 22	      format(i5,2f16.6)
 25	   continue
	endif
	fstop = 0
	pconc = 0
	index = 0
	fcond = 0
	do 200 while(fstop .eq. 0)
	   index = index + 1
	   call vol1(ka,tl(1),p)
	   cmb(1,1) = al(1) + ka*v(1)
	   cmb(1,2) = -al(2)
	   do 30 j = 3,nstage
	      cmb(1,j) = 0
 30	   continue
	   do 45 i = 2,nstage-1
	      do 40 j = 1,i-1
		 call vol1(ka,tl(j),p)
		 prod = ka*e(j)
		 do 35 jp = j+1,i-1
		    prod = prod*(1.0 - e(jp))
 35		 continue
		 cmb(i,j) = prod*(v(i)*(1.0-e(i)) - v(i-1))  
 40	      continue
	      call vol1(ka,tl(i),p)
	      cmb(i,i) = al(i) + ka*e(i)*v(i)
	      cmb(i,i+1) = -al(i+1)
	      do 43 j = i+2,nstage
		 cmb(i,j) = 0
 43	      continue
 45	   continue
	   do 50 j = 1,nstage-1
	      call vol1(ka,tl(j),p)
	      prod = ka*e(j)
	      do 48 jp = j+1,nstage-1
		 prod = prod*(1.0 - e(jp))
 48	      continue
	      cmb(nstage,j) = -prod*vtop
 50	   continue
	   cmb(nstage,nstage) = vtop
	   call invert(ld,nstage,cmb,cinv)
	   do 55 i = 1,nstage
	      a(i) = 0
 55	   continue
	   a(fstage) = xfeed*feed
	   do 65 i = 1,nstage
	      x(i) = 0
	      do 60 j = 1,nstage
		 x(i) = x(i) + cinv(i,j)*a(j)
 60	      continue
	      call bubble(tl(i),x(i),p,treb)
	      tv(i) = tl(i)
	      call vol1(ka,tl(i),p)
	      y(i) = ka*x(i)
 65	   continue
c       
c
c
c
c
c
c
	   call press(x(nstage),y0,tcond,pcond)
	   pc = pcond/760.0
	   if(pcond .gt. 3040.0) then
	      conc = 0.0
	      prof = -25000.00
	      return
	   endif
	   tl(1) = treb
	   tv(1) = treb
	   call cpl(cpc,x(nstage))
	   hd = dinit*cpc*(tcond - tref)
	   call heat(hvap,x(nstage))
	   call dew(tdc,x(nstage),p,treb)
	   if(tcond .ge. tdc) then
	      conc = 0
	      prof = -5000.00
	      fcond = 1
	      return
	   else
	      qc = v(nstage-1)*hvap
	      if(tcond .lt. tl(nstage)) then
		 qc = qc + v(nstage-1)*cpc*(tl(nstage)-tcond)
		   if(qc .gt. qcmax) qc = qcmax
	      endif
	   endif
	   costc = 60.0*time*cwater*qc/100000.0
	   call cpl(cpr,x(1))
	   hb = al(1)*cpr*(treb - tref)
	   call cpl(cpf,xfeed)
	   hf = feed*cpf*(tfeed - tref)
	   if(tfeed .ge. tdf) then 
	      call heat(hvap,x(fstage))
	      hf = hf + feed*hvap
	   endif
	   costf = 60.0*time*csteam*hf/100000.0
	   qr = hd + hb + qc - hf
	    if(qr .gt. qrmax) qr = qrmax  
	   costr = 60.0*time*csteam*qr/100000.0
	   do 80 i = 1,fstage-1
	      i2 = i+1
	      call cpl(cp,x(i2))
	      h = cp*(tl(i2)-tref)
	      call vol1(ka,tl(i),p)
	      y0 = ka*x(i)
	      call cpl(cp,y0)
	      call heat(hvap,y0)
	      hv = cp*(tl(i) - tref) + hvap
	      v(i) = (qr + al(1)*h - hb)/(hv - h)
 80	   continue
	   call cpl(cp,x(fstage+1))
	   h = cp*(tl(fstage+1) - tref)
	   call vol1(ka,tl(fstage),p)
	   y0 = ka*x(fstage)
	   call cpl(cp,y0)
	   call heat(hvap,y0)
	   hv = cp*(tl(fstage) - tref) + hvap
	   call cpl(cp,xfeed)
	   hf = cp*(tfeed - tref)
	   v(fstage) = (qr + feed*(hf-h) + al(1)*h - hb)/(hv-h)
	   nstop = nstage-2
	   do 90 i = nbegin,nstop
	      i2 = i+1
	      call cpl(cp,x(i2))
	      h = cp*(tl(i2)-tref)
	      call vol1(ka,tl(i),p)
	      y0 = ka*x(i)
	      call cpl(cp,y0)
	      call heat(hvap,y0)
	      hv = cp*(tl(i) - tref) + hvap
	      v(i) = (qc + hd - dinit*h)/(hv - h)
 90	   continue
	   i = nstage
	   al(nstage) = reflux*dinit
	   do 100 idex = nbegin,nstage-1
	      ia = i
	      i = i-1
	      ib = i-1
	      al(i) = al(ia) + v(ib) - v(i)
	      call eff(e(i),v(i),y(i),tl(i),p,al(i))
 100	   continue
	   al(fstage) = feed + al(nbegin) + v(fstage-1) - v(fstage)
	   i = fstage
	   do 110 idex = 2,fstage-1
	      ia = i
	      i = i-1
	      ib = i-1
	      al(i) = al(ia) + v(ib) - v(i)
 110	   continue
	   if(ftest .eq. 1) then
	      write(8,120) index
 120	      format(' Thiele Geddes Results for iter:',i5)
	      do 140 i = 1,nstage
		 write(8,130) x(i),y(i),al(i),v(i),tl(i),tv(i)
 130		 format(6f10.4)
 140	      continue
	   endif
	   diff = x(nstage) - pconc
	   if(dabs(diff) .lt. ep) fstop = 1
	   if(index .ge. 30) fstop = 1
	   pconc = x(nstage)
 200	continue
	conc = pconc
	value = 0.0
	if(conc .ge. 0.80 .and. conc .lt. 0.90) value = value1
	if(conc .ge. 0.90 .and. conc .lt. 0.95) value = value2
	if(conc .ge. 0.95) value = value3
	prof = 60.0*time*(dv*value - fv*feedc)
	prof = prof - costc - costf - costr - 1000.0
	return
	end
c
	subroutine press(x,y,t,p)
	implicit real*8 (a-h,o-z)	
	vpa = dexp(20.9 - 5021.994/((t-32.0)/1.8 + 273))
	vpb = dexp(20.5 - 5175.996/((t-32.0)/1.8 + 273))
	p = vpa*x + vpb*(1.0 - x)
	y = vpa*x/p
	return
	end
c
	subroutine vol1(k,t,p)
	implicit real*8 (a-h,o-z)	
	real*8 k
	vpa = dexp(20.9 - 5021.994/((t-32.0)/1.8 + 273))
	k = vpa/p
	return
	end
c
	subroutine vol2(k,t,p)
	implicit real*8 (a-h,o-z)	
	real*8 k
	vpb = dexp(20.5 - 5175.996/((t-32.0)/1.8 + 273))
	k = vpb/p
	return
	end
c
	subroutine eff(e,v,y,t,p,l)
	implicit real*8 (a-h,o-z)
	real*8 l
	e1 = dexp(-0.00475*v)
	ua = 0.25*v*(t + 460.0)/p
	rhog = (46.0*y + 18.0*(1.0-y))*p/(550.0*(t + 460.0))
	fg = ua*dsqrt(rhog)
	if(fg .lt. 0.5) e1 = 0.83 - 1.66*(0.5 - fg)
	if(fg .ge. 0.5 .and. fg .lt. 1.0) e1 = 0.83
	if(fg .gt. 1.0 .and. fg . lt. 5.0) e1 = 0.83 - 0.17*(fg - 0.5)
	if(fg .gt. 5.0) e1 = 0.065
	e2 = dexp(-0.001*l)
	if(e2 .gt. 1.0) e2 = 1.0
	e = e1*e2
	return
	end
c
	subroutine heat(hvap,x)
	implicit real*8 (a-h,o-z)
	hvap = 17000.0*x + 19000.0*(1.0-x)
	return
	end
c
	subroutine cpl(cp,x)
	implicit real*8 (a-h,o-z)
        cp = 31.3*x + 18.0*(1.0-x)
	return
	end
c
	subroutine bubble(btemp,x,p,tinit)
	implicit real*8 (a-h,o-z)
	real*8 ka,kw
	t = tinit
	call vol1(ka,t,p)
	call vol2(kw,t,p)
	sum = ka*x + kw*(1-x)
	do 10 while(sum .gt. 1.0)
	   tprev = t
	   sump = sum
	   t = t - 10.0
	   if(t.lt.50.0) return
	   call vol1(ka,t,p)
	   call vol2(kw,t,p)
	   sum = ka*x + kw*(1-x)
 10	continue
	t = tprev
	sum = sump
	do 20 while(sum .gt. 1.0)
	   tprev  = t
	   sump = sum
	   t = t - 1.0
	   call vol1(ka,t,p)
	   call vol2(kw,t,p)
	   sum = ka*x + kw*(1-x)
 20	continue
	t = tprev
	sum = sump
	do 30 while(sum .gt. 1.0)
	   tprev = t
	   t = t - 0.1
	   call vol1(ka,t,p)
	   call vol2(kw,t,p)
	   sum = ka*x + kw*(1-x)
 30	continue
	btemp = tprev
	return
	end
c
	subroutine dew(dtemp,y,p,tinit)
	implicit real*8 (a-h,o-z)
	real*8 ka,kw
	t = tinit
	call vol1(ka,t,p)
	call vol2(kw,t,p)
	sum = y/ka + (1.0 - y)/kw
	do 10 while(sum .lt. 1.0)
	   tprev = t
	   sump = sum
	   t = t - 10.0
	   call vol1(ka,t,p)
	   call vol2(kw,t,p)
	   sum = y/ka + (1.0 - y)/kw
 10	continue
	t = tprev
	sum = sump
	do 20 while(sum .lt. 1.0)
	   tprev = t
	   sump = sum
	   t = t - 1.0
	   call vol1(ka,t,p)
	   call vol2(kw,t,p)
	   sum = y/ka + (1.0 - y)/kw
 20	continue
	t = tprev
	sum = sump
	do 30 while(sum .lt. 1.0)
	   tprev = t
	   t = t - 0.1
	   call vol1(ka,t,p)
	   call vol2(kw,t,p)
	   sum = y/ka + (1.0 - y)/kw
 30	continue
	dtemp = tprev
	return
	end
c
	SUBROUTINE INVERT(LD,N,A,AINV)
	REAL*8 A(LD,LD),AINV(LD,LD)
	REAL*8 MAX,TEMP
	DO 110 I = 1,N
          DO 100 J = 1,N
             AINV(I,J) = 0
 100      CONTINUE
          AINV(I,I) = 1.0
 110	CONTINUE
	DO 200 I = 1,N-1
          MAX = DABS(A(I,I))
          DO 130 I2 = I+1,N
             IF(DABS(A(I2,I)) .LE. MAX) GO TO 130
             MAX = DABS(A(I2,I))
             DO 120 J = 1,N
                TEMP = A(I,J)
                A(I,J) = A(I2,J)
                A(I2,J) = TEMP
                TEMP = AINV(I,J)
                AINV(I,J) = AINV(I2,J)
                AINV(I2,J) = TEMP   
 120	     CONTINUE
 130	   CONTINUE
          MAX = A(I,I)
          DO 150 I2 = I+1,N
             FACT = A(I2,I)/MAX
             DO 140 J = 1,N
                A(I2,J) = A(I2,J) - FACT*A(I,J)
                AINV(I2,J) = AINV(I2,J) - FACT*AINV(I,J)
 140         CONTINUE
 150      CONTINUE
 200	CONTINUE
	DO 170 I2 = 1,N
          MAX = A(I2,I2)
	  IF(MAX.EQ.0) WRITE(6,151)
           IF(MAX.EQ.0) STOP
           IF(DABS(MAX).LT.0.000000001) WRITE(6,152)
 151       FORMAT(1X,'MATRIX IS SINGULAR')
 152       FORMAT(1X,'MATRIX IS NEAR SINGULAR')
	   DO 160 J = 1,N
	      A(I2,J) = A(I2,J)/MAX
	      AINV(I2,J) = AINV(I2,J)/MAX
 160      CONTINUE
 170    CONTINUE
        N2 = N+1
        DO 220 I = 1,N-1
           I2 = N2-I
           AINV(I2,I2) = AINV(I2,I2)/A(I2,I2)
           A(I2,I2) = 1.0
           DO 215 I3 = 1,I2-1
	      DO 210 J = 1,N     
		 AINV(I3,J) = AINV(I3,J) - A(I3,I2)*AINV(I2,J)
 210       CONTINUE
 215	CONTINUE
 220    CONTINUE
        AINV(1,1) = AINV(1,1)/A(1,1)
        RETURN
        END
c
	SUBROUTINE RANDOM(SEED,N,LD,VALUE)
	IMPLICIT REAL*8 (A-H,O-Z)
	INTEGER SEED
	REAL*8 R(97),VALUE(LD)
	M1 = 259200
	M2 = 134456
	M3 = 243000
	IA1 = 7141
	IA2 = 8121
	IA3 = 4561
	IC1 = 54773
	IC2 = 28411
	IC3 = 51349
	RM1 = 1.0/M1
	RM2 = 1.0/M2
	IF(SEED.GE.0) GO TO 150
	IX1 = MOD(IC1 - SEED,M1)
	IX1 = MOD(IX1,M2)
	IX2 = MOD(IX1,M2)
	IX1 = MOD(IA1*IX1+IC1,M1)
	IX3 = MOD(IX1,M3)
	DO 100 J = 1,97
	   IX1 = MOD(IA1*IX1+IC1,M1)
	   IX2 = MOD(IA2*IX2+IC2,M2)
	   R(J) = (IX1 + IX2*RM2)*RM1
 100	CONTINUE
 150	DO 200 I = 1,N
	   IX1 = MOD(IA1*IX1+IC1,M1)
	   IX2 = MOD(IA2*IX2+IC2,M2)
	   IX3 = MOD(IA3*IX3+IC3,M3)
	   J = 1 + (97*IX3)/M3
	   VALUE(I) = R(J)
	   R(J) = (IX1 + IX2*RM2)*RM1
 200	CONTINUE
	SEED = IX1
	RETURN
	END