您所在位置：主题：【不怕笑话】自编小程序求助，我想毕业。。。
_xiaoxi_
[专家分：0] 发布于 2010-05-27 21:13:00
我是做水动力数学模拟的，自己编了个小程序来算离散方程组，主要是 u k s vt 这四个变量，我调试了很久，程序算出的结果一直发散。由于本人是菜鸟中的菜鸟，格式调试完之后，不知道是不是编程上的不妥，导致在计算过程中会出现 nan。特在此求助，有闲心和善心的的行内人帮忙看看，多谢啦！
前面定义常数部分和罗里吧嗦的各数值没有贴上来。主要迭代计算部分如下：
PS：我这个程序需要的内存空间貌似比较大，nt的值也就能设成100000，有没有办法减小内存占用值，这样我可以把时间步长设小一些，这样说不定可以收敛。
program main
use canshu
implicit none
!水流部分数组
real(kind=8),allocatable::u(:,:),k(:,:),s(:,:),u1(:,:),k1(:,:),s1(:,:),vt(:,:),vt1(:,:) !s()系耗散系数
real(kind=8),allocatable::utt(:),dz(:),dz1(:),z(:)!摩阻流速，网格长
real(kind=8),allocatable::au(:,:),bu(:,:),cucu(:,:),ak(:,:),bk(:,:),ckck(:,:),as(:,:),bs(:,:),cscs(:,:) !系数
real(kind=8),allocatable::du(:,:),dk(:,:),ds(:,:) !方程常数
!real(kind=8),allocatable::f1(:,:),f2(:,:),fu(:,:) !新模型 系数
real(kind=8),allocatable::fufu(:),fk(:),fs(:),gu(:),gk(:),gs(:) !追赶法过渡数组
real,allocatable::cha(:) !判断
real::dt,z0,r,xu1,xu2,xk1,xk2,xs1,xs2 !离散格式重复项
integer::i,j,ji,nz,nt
real::maxcha
real::c1,c2,cu,ck,cs,a !调节参数
call shurucanshu()
call bosu()
call bianjiecenghou()
call yundongnianzhi()
!输入参数值
c1=1.44
c2=1.92
ck=1.0
cs=1.3
cu=0.09
open(unit=1,file="xianshicanshu.txt")
write(1,*)"d=",d
write(1,*) "f=",fw
write(1,*)"c1=",c1,"c2=",c2,"cu=",cu,"ck=",ck,"cs=",cs
!write(*,*)"input 网格总数nz"
!read(*,*)nz
nz=100
write(*,*)"nz=",nz
!write(*,*)"input 时间段数nt"
!read(*,*)nt
nt=10000
dt=t/nt !时间步长
a=0.5 !加权值
allocate(u(0:nz,0:nt),k(1:nz,0:nt),s(1:nz,0:nt),u1(0:nz,0:nt),k1(0:nz,0:nt),s1(0:nz,0:nt),vt(0:nz+1,0:nt+1),vt1(0:nz+1,0:nt+1))
allocate(utt(0:nt),dz(0:nz+1),dz1(0:nz+1),z(0:nz))
allocate(au(1:nz-1,1:nt),bu(1:nz-1,1:nt),cucu(1:nz-1,1:nt),ak(1:nz-1,1:nt),bk(1:nz,1:nt),ckck(1:nz-1,1:nt),as(1:nz-1,1:nt),bs(1:nz-1,1:nt),cscs(1:nz-1,1:nt))
allocate(du(1:nz-1,1:nt),dk(1:nz-1,1:nt),ds(1:nz-1,1:nt))
allocate(fufu(0:nz),fk(0:nz),fs(0:nz),gu(0:nz+1),gk(0:nz+1),gs(0:nz+1))
!allocate(fu(1:nz-1,1:nt),f1(1:nz-1,0:nt),f2(1:nz-1,0:nt))
allocate(cha(0:nz))
!摩阻流速
do j=0,nt
utt(j)=sqrt(abs(0.5*fw*u1m**2*(u01*cos(pi*j/(0.5*nt)-pi/4)+u02*cos(pi*j/(0.25*nt)-pi/2))/u01)+0.5*fw*uc)
!write(*,*) "utt=",utt(j)
end do
!定网格长
z0=ks/30
r=1.1
do i=1,nz
z(i)=d*(r**i-1)/(r**nz-1)
dz(i)=d*(r-1)/(r**nz-1)*r**(i-1) !u
end do
do i=1,nz
dz1(i)=0.5*(dz(i)+dz(i+1)) !k,ε,vt
end do
!初始条件
do i=1,nz
u1(i,0)=0
k1(i,0)=0.001
s1(i,0)=0.002
vt1(i,0)=0.0001
end do
!边界条件
do j=0,nt
u1(0,j)=0
u1(nz,j)=(u01*cos(pi*j/(0.5*nt))+u02*cos(pi*j/(0.25*nt))+uc)*log(z(nz)/z0)/log((z0+d)/z0)
k1(0,j)=utt(j)**2/sqrt(cu)
!write(*,*) "k1 =",k1(0,j)
s1(0,j)=utt(j)**3/(kapa*(0.5*dz(1))+z0)
!write(*,*) s1(0,j)
vt1(0,j)=utt(j)*kapa*(0.5*dz(1)+z0)
!write(*,*)vt1(0,j)
end do
!计算（流速部分）
maxcha=0.00001
do while(maxcha>0.000001)
ji=ji+1
do j=1,2
do i=1,nz-1
!u方程
xu1=(vm+vt1(i+1,j-1))/(dz(i+1)*dz1(i))
xu2=(vm+vt1(i,j-1))/(dz(i)*dz1(i))
au(i,j)=a*dt*xu1
cucu(i,j)=a*dt*xu2
bu(i,j)=1+au(i,j)+cucu(i,j)
du(i,j)=u1(i,j-1)+a*dt*(utt(j)**2/h-2*pi/nt*u01*sin(2*pi*j/nt)-4*pi/nt*u02*sin(4*pi*j/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))&
+(1-a)*dt*(xu1*u1(i+1,j-1)+xu2*u1(i-1,j-1)-xu1*u1(i,j-1)-xu2*u1(i,j-1)+utt(j-1)**2/h-2*pi/nt*u01*sin(2*pi*(j-1)/nt)-4*pi/nt*u02*sin(4*pi*(j-1)/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))
!特殊值末一行
du(nz-1,j)=u1(nz-1,j-1)+a*dt*(utt(j)**2/h-2*pi/nt*u01*sin(2*pi*j/nt)-4*pi/nt*u02*sin(4*pi*j/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))&
+(1-a)*dt*((vm+vt1(nz-1,j-1))/(dz(nz-1)*dz1(nz-1))*u1(i-1,j-1)-(vm+vt1(nz-1,j-1))/(dz(nz-1)*dz1(nz-1))*u1(i,j-1)+utt(j-1)**2/h-2*pi/nt*u01*sin(2*pi*(j-1)/nt)-4*pi/nt*u02*sin(4*pi*(j-1)/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))
fufu(0)=0
fufu(i)=au(i,j)/(bu(i,j)-cucu(i,j)*fufu(i-1))
fufu(nz-1)=0
gu(0)=0
gu(i)=(du(i,j)+cucu(i,j)*gu(i-1))/(bu(i,j)-cucu(i,j)*fufu(i-1))
end do
do i=nz-1,1,-1
u1(i,j)=fufu(i)*u1(i+1,j)+gu(i)
end do
!write(*,*)du(:,1)
continue
do i=1,nz-1
!k方程
xk1=(vt1(i,j-1)+vt1(i+1,j-1)+2*ck*vm)/ck/dz1(i)/dz(i)*0.5
xk2=(vt1(i-1,j-1)+vt1(i,j-1)+2*ck*vm)/ck/dz1(i+1)/dz(i)*0.5
ak(i,j)=a*dt*xk1
ckck(i,j)=a*dt*xk2
bk(i,j)=1+ak(i,j)+ckck(i,j)+s1(i,j-1)/k1(i,j-1)*a*dt
dk(i,j)=k1(i,j-1)+a*dt*vt1(i,j-1)*(u1(i,j)-u1(i-1,j))**2/(dz(i)**2)&
+(1-a)*dt*(xk1*k1(i+1,j-1)+xk2*k1(i-1,j-1)-xk1*k1(i,j-1)-xk2*k1(i,j-1)+vt1(i,j-1)*(u1(i,j-1)-u1(i-1,j-1))**2/dz(i)**2-s1(i,j-1))
!特殊值
dk(nz-1,j)=k1(1,j-1)+a*dt*vt1(1,j-1)*(u1(2,j)-u1(1,j))**2/(dz(1)**2)&
+(1-a)*dt*((vt1(nz-1,j-1)+vt1(nz-2,j-1)+2*ck*vm)/(ck*dz1(nz-2)+2*dz(nz-1))*k1(nz-2,j-1)&
-(vt1(nz-1,j-1)+vt1(nz-2,j-1)+2*ck*vm)/(ck*dz1(nz-2)*2*dz(nz-1))*k1(nz-1,j-1)+vt1(nz-1,j-1)*(u1(nz-1,j-1)-u1(nz-2,j-1))**2/dz(nz-1)**2-s1(nz-1,j-1))
bk(nz-1,j)=1+ckck(nz-1,j)+s1(nz-1,j-1)/k1(nz-1,j-1)*a*dt
!ε方程
!f1(i,j)=1
!f2(i,j)=1
!fu(i,j)=1
!fu(i,j)=(1+4.1/k1(i,j)**1.5*(s1(i,j)*vm)**0.75)*(1-exp(-(vm*s1(i,j))**0.25*z(i)/(vm*15.75)))**2
!f2(i,j)=(1-0.3*exp(-(k1(i,j)**4/(s1(i,j)*vm*6.5)**2)*(1-exp(-vm*s1(i,j)**0.25*z(i)/(vm*3.64)))**2
xs1=(vt1(i,j-1)+vt1(i+1,j-1)+2*cs*vm)/ck/dz1(i)/dz(i)*0.5
xs2=(vt1(i-1,j-1)+vt1(i,j-1)+2*cs*vm)/ck/dz1(i-1)/dz(i)*0.5
as(i,j)=a*dt*xs1
cscs(i,j)=a*dt*xs2
bs(i,j)=1+as(i,j)+cscs(i,j)+a*dt*c2*s1(i,j-1)/k1(i,j-1)
ds(i,j)=s1(i,j-1)+a*dt**c1*s1(i,j-1)/k1(i,j-1)*vt1(i,j-1)*(u1(i,j)-u1(i-1,j))**2/dz(i)**2&
+(1-a)*dt*(xs1*s1(i+1,j-1)+xs2*s1(i-1,j-1)-xs1*s1(i,j-1)-xs2*s1(i,j-1)&
+c1*s1(i,j-1)/k1(i,j-1)*vt1(i,j-1)*(u1(i,j-1)-u1(i-1,j-1))**2/dz(i)**2-c2*s1(i,j-1)**2/k1(i,j-1)**2)
!write(*,*)ds(:,1)
!特殊值
ds(nz-1,j)=s1(nz-1,j-1)+a*dt**c1*s1(nz-1,j-1)/k1(nz-1,j-1)*vt1(nz-1,j-1)*(u1(nz-1,j)-u1(nz-2,j))**2/dz(nz-1)**2&
+(1-a)*dt*((vt1(nz-2,j-1)+vt1(nz-1,j-1)+2*cs*vm)/ck/dz1(nz-2)/dz(i)*0.5*s1(nz-1,j-1)-(vt1(nz-2,j-1)+vt1(nz-1,j-1)+2*cs*vm)/ck/dz1(nz-2)/dz(i)*0.5*s1(nz-1,j-1)&
+c1*s1(nz-1,j-1)/k1(nz-1,j-1)*vt1(nz-1,j-1)*(u1(nz-1,j-1)-u1(nz-2,j-1))**2/dz(nz-1)**2-c2*s1(nz-1,j-1)**2/k1(nz-1,j-1)**2)
bs(nz-1,j)=1+cscs(nz-1,j)+a*dt*c2*s1(nz-1,j-1)/k1(nz-1,j-1)
!追赶法求解
fk(1)=0
fk(i)=ak(i,j)/(bk(i,j)-ckck(i,j)*fk(i-1))
fs(nz)=1
fs(i)=as(i,j)/(bs(i,j)-cscs(i,j)*fs(i-1))
!fs(nz-1)=0
gk(1)=0
gk(i)=(dk(i,j)+ak(i,j)*gk(i-1))/(bk(i,j)-ckck(i,j)*fk(i-1))
gs(1)=0
gs(i)=(ds(i,j)+as(i,j)*gs(i-1))/(bs(i,j)-cscs(i,j)*fs(i-1))
end do
continue
do i=nz-1,1,-1
k1(i,j)=fk(i)*k1(i+1,j)+gk(i)
s1(i,j)=fs(i)*s1(i+1,j)+gs(i)
end do
k1(nz,j)=k1(nz-1,j)
s1(nz,j)=s1(nz-1,j)
continue
do i=1,nz
vt1(i,j)=abs((cu*k1(i,j)**2)/s1(i,j))
write(*,*)vt1(i,j)
end do
continue
end do
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cha(:)=abs(u1(1,nt)-u1(1,0))
!取最大值
maxcha=cha(1)
do i=2,nz-1
if(cha(i)>maxcha)then
maxcha=cha(i)
else
continue
end if
end do
write(*,*)ji,"maxcha",maxcha
!赋值
do i=1,nz
u1(i,0)=u1(i,nt)
k1(i,0)=k1(i,nt)
s1(i,0)=s1(i,nt)
vt1(i,0)=vt1(i,nt)
end do
end do
前面定义常数部分和罗里吧嗦的各数值没有贴上来。主要迭代计算部分如下：
PS：我这个程序需要的内存空间貌似比较大，nt的值也就能设成100000，有没有办法减小内存占用值，这样我可以把时间步长设小一些，这样说不定可以收敛。
program main
use canshu
implicit none
!水流部分数组
real(kind=8),allocatable::u(:,:),k(:,:),s(:,:),u1(:,:),k1(:,:),s1(:,:),vt(:,:),vt1(:,:) !s()系耗散系数
real(kind=8),allocatable::utt(:),dz(:),dz1(:),z(:)!摩阻流速，网格长
real(kind=8),allocatable::au(:,:),bu(:,:),cucu(:,:),ak(:,:),bk(:,:),ckck(:,:),as(:,:),bs(:,:),cscs(:,:) !系数
real(kind=8),allocatable::du(:,:),dk(:,:),ds(:,:) !方程常数
!real(kind=8),allocatable::f1(:,:),f2(:,:),fu(:,:) !新模型 系数
real(kind=8),allocatable::fufu(:),fk(:),fs(:),gu(:),gk(:),gs(:) !追赶法过渡数组
real,allocatable::cha(:) !判断
real::dt,z0,r,xu1,xu2,xk1,xk2,xs1,xs2 !离散格式重复项
integer::i,j,ji,nz,nt
real::maxcha
real::c1,c2,cu,ck,cs,a !调节参数
call shurucanshu()
call bosu()
call bianjiecenghou()
call yundongnianzhi()
!输入参数值
c1=1.44
c2=1.92
ck=1.0
cs=1.3
cu=0.09
open(unit=1,file="xianshicanshu.txt")
write(1,*)"d=",d
write(1,*) "f=",fw
write(1,*)"c1=",c1,"c2=",c2,"cu=",cu,"ck=",ck,"cs=",cs
!write(*,*)"input 网格总数nz"
!read(*,*)nz
nz=100
write(*,*)"nz=",nz
!write(*,*)"input 时间段数nt"
!read(*,*)nt
nt=10000
dt=t/nt !时间步长
a=0.5 !加权值
allocate(u(0:nz,0:nt),k(1:nz,0:nt),s(1:nz,0:nt),u1(0:nz,0:nt),k1(0:nz,0:nt),s1(0:nz,0:nt),vt(0:nz+1,0:nt+1),vt1(0:nz+1,0:nt+1))
allocate(utt(0:nt),dz(0:nz+1),dz1(0:nz+1),z(0:nz))
allocate(au(1:nz-1,1:nt),bu(1:nz-1,1:nt),cucu(1:nz-1,1:nt),ak(1:nz-1,1:nt),bk(1:nz,1:nt),ckck(1:nz-1,1:nt),as(1:nz-1,1:nt),bs(1:nz-1,1:nt),cscs(1:nz-1,1:nt))
allocate(du(1:nz-1,1:nt),dk(1:nz-1,1:nt),ds(1:nz-1,1:nt))
allocate(fufu(0:nz),fk(0:nz),fs(0:nz),gu(0:nz+1),gk(0:nz+1),gs(0:nz+1))
!allocate(fu(1:nz-1,1:nt),f1(1:nz-1,0:nt),f2(1:nz-1,0:nt))
allocate(cha(0:nz))
!摩阻流速
do j=0,nt
utt(j)=sqrt(abs(0.5*fw*u1m**2*(u01*cos(pi*j/(0.5*nt)-pi/4)+u02*cos(pi*j/(0.25*nt)-pi/2))/u01)+0.5*fw*uc)
!write(*,*) "utt=",utt(j)
end do
!定网格长
z0=ks/30
r=1.1
do i=1,nz
z(i)=d*(r**i-1)/(r**nz-1)
dz(i)=d*(r-1)/(r**nz-1)*r**(i-1) !u
end do
do i=1,nz
dz1(i)=0.5*(dz(i)+dz(i+1)) !k,ε,vt
end do
!初始条件
do i=1,nz
u1(i,0)=0
k1(i,0)=0.001
s1(i,0)=0.002
vt1(i,0)=0.0001
end do
!边界条件
do j=0,nt
u1(0,j)=0
u1(nz,j)=(u01*cos(pi*j/(0.5*nt))+u02*cos(pi*j/(0.25*nt))+uc)*log(z(nz)/z0)/log((z0+d)/z0)
k1(0,j)=utt(j)**2/sqrt(cu)
!write(*,*) "k1 =",k1(0,j)
s1(0,j)=utt(j)**3/(kapa*(0.5*dz(1))+z0)
!write(*,*) s1(0,j)
vt1(0,j)=utt(j)*kapa*(0.5*dz(1)+z0)
!write(*,*)vt1(0,j)
end do
!计算（流速部分）
maxcha=0.00001
do while(maxcha>0.000001)
ji=ji+1
do j=1,2
do i=1,nz-1
!u方程
xu1=(vm+vt1(i+1,j-1))/(dz(i+1)*dz1(i))
xu2=(vm+vt1(i,j-1))/(dz(i)*dz1(i))
au(i,j)=a*dt*xu1
cucu(i,j)=a*dt*xu2
bu(i,j)=1+au(i,j)+cucu(i,j)
du(i,j)=u1(i,j-1)+a*dt*(utt(j)**2/h-2*pi/nt*u01*sin(2*pi*j/nt)-4*pi/nt*u02*sin(4*pi*j/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))&
+(1-a)*dt*(xu1*u1(i+1,j-1)+xu2*u1(i-1,j-1)-xu1*u1(i,j-1)-xu2*u1(i,j-1)+utt(j-1)**2/h-2*pi/nt*u01*sin(2*pi*(j-1)/nt)-4*pi/nt*u02*sin(4*pi*(j-1)/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))
!特殊值末一行
du(nz-1,j)=u1(nz-1,j-1)+a*dt*(utt(j)**2/h-2*pi/nt*u01*sin(2*pi*j/nt)-4*pi/nt*u02*sin(4*pi*j/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))&
+(1-a)*dt*((vm+vt1(nz-1,j-1))/(dz(nz-1)*dz1(nz-1))*u1(i-1,j-1)-(vm+vt1(nz-1,j-1))/(dz(nz-1)*dz1(nz-1))*u1(i,j-1)+utt(j-1)**2/h-2*pi/nt*u01*sin(2*pi*(j-1)/nt)-4*pi/nt*u02*sin(4*pi*(j-1)/nt)+u01*cos(pi*(j-1)*2/nt)+u02*cos(pi*(j-1)*4/nt))
fufu(0)=0
fufu(i)=au(i,j)/(bu(i,j)-cucu(i,j)*fufu(i-1))
fufu(nz-1)=0
gu(0)=0
gu(i)=(du(i,j)+cucu(i,j)*gu(i-1))/(bu(i,j)-cucu(i,j)*fufu(i-1))
end do
do i=nz-1,1,-1
u1(i,j)=fufu(i)*u1(i+1,j)+gu(i)
end do
!write(*,*)du(:,1)
continue
do i=1,nz-1
!k方程
xk1=(vt1(i,j-1)+vt1(i+1,j-1)+2*ck*vm)/ck/dz1(i)/dz(i)*0.5
xk2=(vt1(i-1,j-1)+vt1(i,j-1)+2*ck*vm)/ck/dz1(i+1)/dz(i)*0.5
ak(i,j)=a*dt*xk1
ckck(i,j)=a*dt*xk2
bk(i,j)=1+ak(i,j)+ckck(i,j)+s1(i,j-1)/k1(i,j-1)*a*dt
dk(i,j)=k1(i,j-1)+a*dt*vt1(i,j-1)*(u1(i,j)-u1(i-1,j))**2/(dz(i)**2)&
+(1-a)*dt*(xk1*k1(i+1,j-1)+xk2*k1(i-1,j-1)-xk1*k1(i,j-1)-xk2*k1(i,j-1)+vt1(i,j-1)*(u1(i,j-1)-u1(i-1,j-1))**2/dz(i)**2-s1(i,j-1))
!特殊值
dk(nz-1,j)=k1(1,j-1)+a*dt*vt1(1,j-1)*(u1(2,j)-u1(1,j))**2/(dz(1)**2)&
+(1-a)*dt*((vt1(nz-1,j-1)+vt1(nz-2,j-1)+2*ck*vm)/(ck*dz1(nz-2)+2*dz(nz-1))*k1(nz-2,j-1)&
-(vt1(nz-1,j-1)+vt1(nz-2,j-1)+2*ck*vm)/(ck*dz1(nz-2)*2*dz(nz-1))*k1(nz-1,j-1)+vt1(nz-1,j-1)*(u1(nz-1,j-1)-u1(nz-2,j-1))**2/dz(nz-1)**2-s1(nz-1,j-1))
bk(nz-1,j)=1+ckck(nz-1,j)+s1(nz-1,j-1)/k1(nz-1,j-1)*a*dt
!ε方程
!f1(i,j)=1
!f2(i,j)=1
!fu(i,j)=1
!fu(i,j)=(1+4.1/k1(i,j)**1.5*(s1(i,j)*vm)**0.75)*(1-exp(-(vm*s1(i,j))**0.25*z(i)/(vm*15.75)))**2
!f2(i,j)=(1-0.3*exp(-(k1(i,j)**4/(s1(i,j)*vm*6.5)**2)*(1-exp(-vm*s1(i,j)**0.25*z(i)/(vm*3.64)))**2
xs1=(vt1(i,j-1)+vt1(i+1,j-1)+2*cs*vm)/ck/dz1(i)/dz(i)*0.5
xs2=(vt1(i-1,j-1)+vt1(i,j-1)+2*cs*vm)/ck/dz1(i-1)/dz(i)*0.5
as(i,j)=a*dt*xs1
cscs(i,j)=a*dt*xs2
bs(i,j)=1+as(i,j)+cscs(i,j)+a*dt*c2*s1(i,j-1)/k1(i,j-1)
ds(i,j)=s1(i,j-1)+a*dt**c1*s1(i,j-1)/k1(i,j-1)*vt1(i,j-1)*(u1(i,j)-u1(i-1,j))**2/dz(i)**2&
+(1-a)*dt*(xs1*s1(i+1,j-1)+xs2*s1(i-1,j-1)-xs1*s1(i,j-1)-xs2*s1(i,j-1)&
+c1*s1(i,j-1)/k1(i,j-1)*vt1(i,j-1)*(u1(i,j-1)-u1(i-1,j-1))**2/dz(i)**2-c2*s1(i,j-1)**2/k1(i,j-1)**2)
!write(*,*)ds(:,1)
!特殊值
ds(nz-1,j)=s1(nz-1,j-1)+a*dt**c1*s1(nz-1,j-1)/k1(nz-1,j-1)*vt1(nz-1,j-1)*(u1(nz-1,j)-u1(nz-2,j))**2/dz(nz-1)**2&
+(1-a)*dt*((vt1(nz-2,j-1)+vt1(nz-1,j-1)+2*cs*vm)/ck/dz1(nz-2)/dz(i)*0.5*s1(nz-1,j-1)-(vt1(nz-2,j-1)+vt1(nz-1,j-1)+2*cs*vm)/ck/dz1(nz-2)/dz(i)*0.5*s1(nz-1,j-1)&
+c1*s1(nz-1,j-1)/k1(nz-1,j-1)*vt1(nz-1,j-1)*(u1(nz-1,j-1)-u1(nz-2,j-1))**2/dz(nz-1)**2-c2*s1(nz-1,j-1)**2/k1(nz-1,j-1)**2)
bs(nz-1,j)=1+cscs(nz-1,j)+a*dt*c2*s1(nz-1,j-1)/k1(nz-1,j-1)
!追赶法求解
fk(1)=0
fk(i)=ak(i,j)/(bk(i,j)-ckck(i,j)*fk(i-1))
fs(nz)=1
fs(i)=as(i,j)/(bs(i,j)-cscs(i,j)*fs(i-1))
!fs(nz-1)=0
gk(1)=0
gk(i)=(dk(i,j)+ak(i,j)*gk(i-1))/(bk(i,j)-ckck(i,j)*fk(i-1))
gs(1)=0
gs(i)=(ds(i,j)+as(i,j)*gs(i-1))/(bs(i,j)-cscs(i,j)*fs(i-1))
end do
continue
do i=nz-1,1,-1
k1(i,j)=fk(i)*k1(i+1,j)+gk(i)
s1(i,j)=fs(i)*s1(i+1,j)+gs(i)
end do
k1(nz,j)=k1(nz-1,j)
s1(nz,j)=s1(nz-1,j)
continue
do i=1,nz
vt1(i,j)=abs((cu*k1(i,j)**2)/s1(i,j))
write(*,*)vt1(i,j)
end do
continue
end do
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cha(:)=abs(u1(1,nt)-u1(1,0))
!取最大值
maxcha=cha(1)
do i=2,nz-1
if(cha(i)>maxcha)then
maxcha=cha(i)
else
continue
end if
end do
write(*,*)ji,"maxcha",maxcha
!赋值
do i=1,nz
u1(i,0)=u1(i,nt)
k1(i,0)=k1(i,nt)
s1(i,0)=s1(i,nt)
vt1(i,0)=vt1(i,nt)
end do
end do