forrt1:severe(161)-Fortran语法讨论-专业Fortran论坛

Nimo 发表于 2017-10-25 10:34:04

forrt1:severe(161)

program main
implicit none
real p,V0,f
!real RHS
real ox,ot
real De,Dp,qe,bt,KB
!real RS,git,rit
real miue,miup
real icjifen,vol,gama
real s0, sb,d1,d2,l,Dg
integer i,j,ii,k,error
!real am,bm,cm,dm,em

real xo(1001,1)
real ne(1001,2),np(1001,2),E(1000,2)
real g(1000,2),R(1000,2),gg(1000,2),RR(1000,2) ,gte(1000,2),Rte(1000,2)
real je(1000,1),jp(1000,1),jq(1000,1),XL(1000,2),XLL(1000,2),Xte(1000,2)
real IC(1000,1)
real B(999,1),A(998,1),C(998,1),S(999,1)
real Bte(999,1),Ate(998,1),Cte(998,1),Ste(999,1)

real BB(999,1),AA(998,1),CC(998,1),SS(999,1),T(1,1),TT(1,1),ITz(1,1) ,Tte(1,1)
real vg(1,1),va(1,1),bdh1(2,1),bdh2(2,1),bdh11(2,1)
real Te(1000,1)

V0=1.5e3
P=760.0
F=3.3e4

s0=8.85E-14
bt=1.389e-28
sb=7.5
D1=0.1
D2=0.1
Dg=0.2
L=0.2267
KB=1.380605e-23
gama=0.01
ox=2.0E-4
ot=1.0E-9

miue=4.863e2
De=7.2944e2
miup=10.4
Dp=0.53947
qe=1.6E-19

!!!!!!!!!!!!!!!!!!!!!!空间格点的坐标
j=1
xo(1,1)=0.1d0
do i=2,1001
xo(i,1)=xo(i-1,1)+ox
end do
!!!!!!!!!!!!!!!!!!!!!给ne np赋初值，即t=0时刻各点ne的值
do i=1,1001
ne(i,1)=1.0d7
np(i,1)=1.0d7
te(i,1)=0.026
end do

!!!!!!!!!!!!!!!!!!!!!给E赋值，即t=0时刻各点E的值
do i=1,1000
E(i,1)= 1.0d7*qe
end do
do while(.true.) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!
va(1,1)=V0*sin(2.0d0*3.14d0*f*real(j)*ot)
do i=1,1000
g(i,1)=miue*ox*E(i,1)/De
R(i,1)=(exp(g(i,1))-1)/g(i,1)
gg(i,1)=-miup*ox*E(i,1)/Dp
RR(i,1)=(exp(gg(i,1))-1)/gg(i,1)
gTE(i,1)=miue*ox*E(i,1)/De
RTE(i,1)=(exp(gTE(i,1))-1)/gTE(i,1)
end do
!!!!!!!!!!!!追赶法求密度
do i=1,1000
XL(i,1)=De*ot/(ox*ox*R(i,1))
XLL(i,1)=Dp*ot/(ox*ox*RR(i,1))
XTE(i,1)=5*De*ot /(3*ox*ox*RTE(i,1))
end do
do ii=1,998
   B(ii,1)=1+XL(ii+1,1)+XL(ii,1)*exp(g(ii,1))
   A(ii,1)=-XL(ii+1,1)
   C(ii,1)=-XL(ii+1,1)*exp(g(ii+1,1))
   S(ii,1)=ne(ii+1,1)+4.4*p*exp(-14.0*(p/abs(E(ii+1,1)))**0.5)*miue*abs(E(ii+1,1))*ot*ne(ii+1,1)-bt*ot*ne(ii+1,1)*ne(ii+1,1)*np(ii+1,1)
   BB(ii,1)=1+XLL(ii+1,1)+XLL(ii,1)*exp(gg(ii,1))
   AA(ii,1)=-XLL(ii+1,1)
   CC(ii,1)=-XLL(ii+1,1)*exp(gg(ii+1,1))
   SS(ii,1)=np(ii+1,1)+4.4*p*exp(-14.0*(p/abs(E(ii+1,1)))**0.5)*miue*abs(E(ii+1,1))*ot*ne(ii+1,1)-bt*ot*ne(ii+1,1)*ne(ii+1,1)*np(ii+1,1)
BTE(ii,1)=1+XTE(ii+1,1)+XTE(ii,1)*exp(gTE(ii,1))
ATE(ii,1)=-XTE(ii+1,1)
CTE(ii,1)=-XTE(ii+1,1)*exp(gTE(ii+1,1))
   STE(ii,1)=ne(ii+1,1)-qe*abs(E(ii+1,1))*je(ii+1,1)*ot-bt*ot*ne(ii+1,1)*ne(ii+1,1)*np(ii+1,1)*(3*KB*TE(ii+1,1)/2 )
end do

end do
B(999,1)=1+XL(1000,1)+XL(999,1)*exp(g(999,1))
S(999,1)=ne(1000,1)+4.4*p*exp(-14.0*(p/abs(E(1000,1)))**0.5)*miue*abs(E(1000,1))*ot*ne(1000,1)-bt*ot*ne(1000,1)*ne(1000,1)*np(1000,1)
ne(1,2)=1.0d7
ne(1001,2)=1.0d7
S(1,1)=S(1,1)+XL(1,1)*ne(1,2)
S(999,1)=S(999,1)+XL(1000,1)*exp(g(1000,1))*ne(1001,2)!!!!!!!!!!!!!!!!!设置电子离子密度的边界值
BB(999,1)=1+XLL(1000,1)+XLL(999,1)*exp(gg(999,1))
SS(999,1)=np(1000,1)+4.4*p*exp(-14.0*(p/abs(E(1000,1)))**0.5)*miue*abs(E(1000,1))*ot*ne(1000,1)-bt*ot*ne(1000,1)*ne(1000,1)*np(1000,1)
np(1,2)=np(2,1)
np(1001,2)=np(1000,1)
SS(1,1)=SS(1,1)+XLL(1,1)*np(1,2)
SS(999,1)=SS(999,1)+XLL(1000,1)*exp(gg(1000,1))*np(1001,2)
BTE(999,1)=1+XTE(1000,1)+XTE(999,1)*exp(gTE(999,1))
STE(999,1)=-qe*abs(E(1000,1))*ot*je(1000,1)-bt*ot*ne(1000,1)*ne(1000,1)*np(1000,1)*3*KB*TE(1000,1)/2
TE(1,2)=1.5
TE(1001,2)=1.5
STE(1,1)=STE(1,1)+XTE(1,1)*TE(1,2)
STE(999,1)=STE(999,1)+XTE(1000,1)*exp(gTE(1000,1))*TE(1001,2)!!!!!!!!
!追赶法
S(1,1)=S(1,1)/B(1,1)
T=B(1,1)
SS(1,1)=SS(1,1)/BB(1,1)
TT=BB(1,1)
STE(1,1)=STE(1,1)/BTE(1,1)
TTE=BTE(1,1)

do k=2,999
B(k-1,1)=C(k-1,1)/T(1,1)
T(1,1)=B(k,1)-A(k-1,1)*B(k-1,1)
S(k,1)=(S(k,1)-A(k-1,1)*S(k-1,1))/T(1,1)
BB(k-1,1)=CC(k-1,1)/TT(1,1)
TT(1,1)=BB(k,1)-AA(k-1,1)*BB(k-1,1)
SS(k,1)=(SS(k,1)-AA(k-1,1)*SS(k-1,1))/TT(1,1)
BTE(k-1,1)=CTE(k-1,1)/TTE(1,1)
TTE(1,1)=BTE(k,1)-ATE(k-1,1)*BTE(k-1,1)
STE(k,1)=(STE(k,1)-ATE(k-1,1)*STE(k-1,1))/TTE(1,1)

end do

do k=1,998
S(999-k,1)=S(999-k,1)-B(999-k,1)*S(1000-k,1)
SS(999-k,1)=SS(999-k,1)-BB(999-k,1)*SS(1000-k,1)
STE(999-k,1)=STE(999-k,1)-BTE(999-k,1)*STE(1000-k,1)

end do
do i=2,1000
ne(i,2)=S(i-1,1)
np(i,2)=SS(i-1,1)
TE(i,2)=STE(i-1,1)
end do
!过河拆桥
ne(:,1)=ne(:,2)
np(:,1)=np(:,2)
TE(:,1)=TE(:,2)

if(vg(1,1)>=0)then
   do i=1,1000
      je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))

      jp(i,1)=Dp*(np(i,1)-exp(gg(i,1))*np(i+1,1))/(ox*RR(i,1))
jq(1,1)=3*(De*(ne(i,1)-ne(i+1,1)*exp(gte(i,1))))/5*(ox*Rte(i,1))
   end do
if(jp(1000,1)>0)then
   gama=0.01d0
else
   gama=0.0d0
end if
   je(1000,1)=je(1000,1)-gama*jp(1000,1)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!je jp的边界条件

else
   do i=1,1000
      je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))
jq(i,1)= 3*(De*(ne(i,1)-ne(i+1,1)*exp(gte(i,1))))/5*(ox*Rte(i,1))

      jp(i,1)=Dp*(np(i,1)-exp(gg(i,1))*np(i+1,1))/(ox*RR(i,1))
   end do
if(jp(1000,1)<0)then
   gama=0.01d0
else
   gama=0.0d0
end if
   je(1,1)=je(1,1)-gama*jp(1,1)

end if

   iC=qe*(jp-je)
   icjifen=ox*sum(iC)
   vol=s0*(V0*sin(2.0d0*3.14d0*f*real(j)*ot)-V0*sin(2.0d0*3.14d0*f*real(j-1)*ot))/ot
   ITz(1,1)=(vol+icjifen)/L
   do i=1,1000
      E(i,1)=ot*(ITz(1,1)-iC(i,1))/s0+E(i,1)
   end do

   vg(1,1)=ox*sum(E(:,1))!!!!!2时刻的气体电压

   bdh1(2,1)=(bdh1(1,1)+ic(1,1)*ot)

bdh2(2,1)=(bdh2(1,1)+ic(1000,1)*ot)

bdh1(1,1)=bdh1(2,1)
bdh11(1,1)=-bdh1(1,1)
bdh2(1,1)=bdh2(2,1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

if(mod(j,100)==0)then
open(1,file="ITz.dat")
open(2,file="va.dat")
open(3,file="vg.dat")
open(4,file="bdh11.dat")
open(5,file="bdh2.dat")
open(6,file="Te.dat")

   write(1,*) j,ITz(1,1)
   write(2,*) j,va(1,1)
   write(3,*) j,vg(1,1)
   write(4,*) j,bdh11(1,1)
   write(5,*) j,bdh2(1,1)
write(6,*) j,Te(1,1)

end if

j=j+1
write(*,*)j

close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
stop
end program

Nimo 发表于 2017-10-25 10:35:22

求助，要加一个电子能量守恒方程式求电子温度
重粒子守恒方程求气体温度

li913 发表于 2017-10-26 13:12:39

绝大多数人不了解你的专业，对这个算法更是不懂。你需要提供详细的算法说明，比如公式、文本介绍。

页: [1]

Fortran Coder's Archiver

forrt1:severe(161)