fortran语言 一维稳态 程序编写

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='one-t1.dat')

write(*,*)'pass'

read(11,*) m1,xi,xe,ls,ki,ke,ai,bi,ae0,be,ti,te

read(11,*) sc,sp,st

read(11,*) ck,b

read(11,*) mode,xpr1,xpr2

close(11)

write(*,*)'pass1'

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,st,ck,b'

read(*,*) Sc,Sp,st,ck,b

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)

! *****************************************************

! ************** !单向加密或变稀********************

! 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)=xf(i)+0.5*(xe-xi)/(m1-2)

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.

do 25 i=1,m1

!write(12,*) xf(i),y

write(12,*)'x(',i,')= ',x(i)

25 write(*,*)'x(',i,')= ',x(i)

close(12)

!c ************ initial vailu**********

do 30 i=1,m1

ta(i)=20.0

30 t(i)=20.0

write(*,*)'pass'

!c ************bundry condition******** 99 if(ls.eq.1) goto 200

do 50 i=1,m1

gm(i)=ck +b*ta(i)

cn(i)=Sc+st*ta(i)

50 ap(i)=Sp! s*

goto 210

200 continue

do 100 i=1,m1

gm(i)=ck

cn(i)=Sc

ap(i)=Sp

100 continue

210 continue

if(mode.eq.2) ex=1

if(mode.eq.1) ex=0

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)-aw(i)*p(i-1))

300 continue

t(m1)=q(m1)

do 400 j=1,m2

i=m1-j

400 t(i)=p(i)*t(i+1)+q(i)

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