[求助] 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