dimension x(100),xf(100),xm(100),xp(100),ap(100),ae(100),aw(100)
dimension cn(100),gm(100),t(100),ta(100),p(100),q(100),r(100),rm(100)
write(*,*)'do you want input data from tabekey(:1) or file(:2)'
read(*,*) ns
if(ns.eq.2) then
open(unit=11,file='onet.dat')
read(11,*) m1,xi,xe,ls,ki,ke,ai,bi,ae0,be,ti,te
read(11,*) sc,sp
read(11,*) ck
read(11,*) mode,xpr1,xpr2
close(11)
else
write(*,*) 'input m1,xi,xe,ls,ki,ke,ai,bi,ae0,be,ti,te'
read(*,*) m1,xi,xe,ls,ki,ke,ai,bi,ae0,be,ti,te
write(*,*) 'input Sc,Sp,ck'
read(*,*) Sc,Sp,ck
write(*,*) 'input mode'
read(*,*) mode
write(*,*)'input xpr1,xpr2'
read(*,*) xpr1,xpr2
endif
!c ************ grid************
m2=m1-1
do 10 i=2,m1
10 xf(i)=xi+(xe-xi)*(i-2)/(m2-1)
! *********************)/(m********************************
! ************** !单向加密-或变稀 ********************)
! do 10 i=2,m1
! 10 xf(i)=xi+(xe-xi)*(float(i-2)/float(m2-1))**xpr1
! *****************************************************
! ************** 双向加密或变稀 ***********************
! do 10 i=2,m1/2
! 10 xf(i)=xi+(xe-xi)*(float(i-2)/float(m2-1))**xpr1
! do 15 i=m1/2+1,m1
! 15 xf(i)=xf(m1/2)+(xe-xf(m1/2))*(float(i-m1/2)/float(m1-m1/2))**xpr2
! *****************************************************************
x(1)=xi
x(m1)=xe
r(1)=x(1)
r(m1)=x(m1)
rm(m1)=xf(m1)
do 20 i=2,m2
x(i)=0.5*(xf(i+1)+xf(1))
r(i)=x(i)
rm(i)=xf(i)
xp(i)=xf(i+1)-x(i)
20 xm(i)=x(i)-xf(i)
if(mode.eq.1) then
do 22 i=1,m1
r(i)=1.
rm(i)=1.
22 continue
endif
open(unit=12,file='grit.dat')
y=1.
Xf(1)=0
do 25 i=1,m1
write(12,*) xf(i),y
25 write(*,*)'x(',i,')= ',x(i)
close(12)
!c ************ initial vailu**********
do 30 i=2,m1
ta(i)=20.0
30 t(i)=20.0
!c ************bundry condition********
b=1
99 if(ls.eq.1) goto 200
do 50 i=1,m1
gm(i)=ck+b*ta(i)
cn(i)=Sc
50 ap(i)=Sp
goto 210
200 continue
do 100 i=1,m1
gm(i)=ck
cn(i)=Sc
ap(i)=Sp
100 continue
210 continue
ex=1
if(mode.eq.3) ex=2
if(ki.gt.1) goto 120
ae(1)=0.
aw(1)=0.
ap(1)=1.
cn(1)=ti
goto 130
120 ae(1)=gm(1)/xm(2)
aw(1)=0.
ap(1)=ae(1)+aw(1)+bi
cn(1)=ai
130 if(ke.gt.1) goto 160
ae(m1)=0.
aw(m1)=0.
ap(m1)=1
cn(m1)=te
goto 180
160 aw(m1)=gm(m1)/xp(m2)
ae(m1)=0.
cn(m1)=ae0
ap(m1)=ae(m1)+aw(m1)+be
180 continue
!c ********** INTERNAL POINT COEFFICIENTS ******************
aw(2)=gm(2)/xm(2)*(0.5*(r(1)+r(2)))**ex
ae(m2)=gm(m2)/xp(m2)*rm(m1)**ex
do 150 i=2,m2-1
ae(i)=1/(xp(i)/gm(i)+xm(i+1)/gm(i+1))*rm(i+1)**ex
aw(I+1)=ae(I)
150 continue
do 250 i=2,m2
cn(i)=cn(i)*(xf(I+1)-xf(I))*r(i)**ex
ap(i)=ae(I)+aw(I)-ap(I)*(xf(I+1)-xf(I))*r(i)**ex
250 continue
!c **************TDMA**********
p(1)=ae(1)/ap(1)
q(1)=cn(1)/ap(1)
do 300 i=2,m1
p(i)=ae(i)/(ap(i)-aw(i)*p(i-1))
q(i)=(cn(i)+aw(i)*q(i-1))/(ap(i)*p(i-1))
300 continue
t(m1)=q(m1)
do 400 i=1,m2
k=m1-i
t(k)=p(k)*t(k+1)+q(k)
400 continue
!c ****************************
if(ls.eq.1) goto 500
df=0.
do 410 i=1,m1
ds=abs((t(i)-ta(i))/t(i))
if(ds.ge.df) df=ds
410 continue
if(df.lt.0.001) goto 500
do 450 i=1,m1
450 ta(i)=ta(i)+0.5*(t(i)-ta(i))
goto 99
500 continue
do 600 i=1,m1
600 write(*,555) i,t(i)
555 format(5x,'T(',i2,')=',f10.3)
stop
end
