c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
c **************************************************************
c
double precision function avg(Y,f,N)
implicit double precision(a-h,o-z)
dimension f(N), Y(N)
avg = 0.d0
do i = 1, N
avg = avg + Y(i)*f(i)
enddo
return
end
c **************************************************************
c
double precision function avgtabip(T,Y,tab,N)
implicit double precision(a-h,o-z)
include "ck.i"
dimension f(lensp), Y(lensp), tab(lensp,TabS:TabE)
itemp = int(T*TABFAC)
if (itemp.lt.TabS .or. itemp+1.gt.TabE) then
write(6,*) 'Error: ',itemp,' not in avgtabip'
end if
call tabintp(T,f,tab,N)
avgtabip = avg(Y,f,N)
return
end
c **************************************************************
c
subroutine tabintp(T,f,tab,N)
implicit double precision(a-h,o-z)
include "ck.i"
dimension f(N), tab(lensp,TabS:TabE)
itemp = int(T*TABFAC)
if (itemp.lt.TabS .or. itemp+1.gt.TabE) then
write(6,*) 'Error: ',itemp,' not in tabular'
goto 200
end if
alpha = (T - itemp/TABFAC) * TABFAC
do i = 1, N
f(i) = (1.d0-alpha)*tab(i,itemp) + alpha*tab(i,itemp+1)
enddo
200 continue
return
end
c **************************************************************
c
subroutine tabsinit(it)
implicit double precision(a-h,o-z)
include "ck.i"
do itemp = TabS, TabE
delt = dble(itemp-it)/TABFAC
do k = 1, Nsp
cpk(k,itemp) = cpinit(k)
hms(k,itemp) = hinit(k) + delt * cpinit(k)
enddo
enddo
call tabrange()
210 continue
return
end
c **************************************************************
c
subroutine tabrange()
implicit double precision(a-h,o-z)
common /vtabrange/ Trmin, Trmax
include "ck.i"
dimension Yk(leNsp)
Trmin = TabS/TABFAC
Trmax = (TabE-1)/TABFAC
do m = 1, Nsp
do k = 1, Nsp
Yk(k) = 0.d0
enddo
Yk(m) = 1.d0
dso = dsign(1.d0, RU/Wk(m)-avgtabip(Trmin,Yk,cpk,Nsp))
do itemp = TabS+1, TabE-1
T = itemp/TABFAC
ds = dsign(1.d0, RU/Wk(m)-avgtabip(T,Yk,cpk,Nsp))
if (ds.ne.dso) Trmax = dmin1(Trmax,(itemp-1)/TABFAC)
dso = ds
enddo
enddo
if (Trmax.le.(TabE-2)/TABFAC)
& write (6,600) Trmin,Trmax,Trmin,(TabE-1)/TABFAC
600 format('Using temperature range ',f9.2,'-',f9.2,
& ' instead of ',f9.2,'-',f9.2)
return
end