自己编写的交错网格有限差分的fortran程序，显示没有语言错误但是运行没有结果且界面闪退
自己编写的交错网格有限差分的fortran程序，显示没有语言错误但是运行没有结果，出现如图情况，请求各位帮忙解决，急，谢谢！
program main
real dx,dz !网格间距
integer pml,modx,modz,n !吸收层厚度，模型网格大小（包含PML层），差分阶数
integer mx,mz !计算模型网格大小（包含PML层）
integer sx,sz,ns,ds,rd,off,dj !震源坐标,ns炮数，ds炮间隔,rd接收道数，最小偏移距,道距
integer tt,nt,out_t !采样时长,采样数
real mf,dt !子波主频,间隔（ms）
integer,allocatable::nsx(:),nsz(:)
real,allocatable::v(:,:),st(:,:)
real,allocatable::jl(:,:,:),jlw(:,:),out_w(:,:)
real start, finish
character*40 outfile,input_filename,outfile2
integer::i,j,k
!输入参数
input_filename='cmd.txt'
call input_cmd(input_filename,dx,dz,pml,n,modz,modx,tt,dt,ns,ds,rd,off,dj,sz,sx,mf,out_t)
!差分网格大小
mx=n*2+modx
mz=n*2+modz
!采样点数
nt=int(tt/dt)
allocate(nsx(ns),nsz(ns))
do i=1,ns
nsx(i)=sx+(i-1)*ds
nsz(i)=sz+n
end do
allocate(v(mz,mx))
allocate(jl(ns,nt,rd),jlw(nt,rd),out_w(mz,mx)) !jl(炮数，采样数，接收道数)
allocate(st(mz,mx))
!输入速度模型并标记地层
v=0
call vr(v,st,mz,mx,n)
do i=1,ns
call calculate_jcwg(i,pml,mx,mz,n,nsx(i),nsz(i),rd,off,dj,dt,nt,mf,dx,dz,v,st,j1w,out_w,out_t)
forall(j=1:nt,k=1:rd)
jl(i,j,k)=jlw(j,k)
end forall
end do
outfile='jlw.sgy'
outfile2='w.sgy'
!sgy格式输出数据
call out_sgy(nt,rd,ns,jl,outfile,dt)
call out_sgy(mz,mx,1,out_w,outfile2,dz)
end program main
subroutine input_cmd(input_filename,dx,dz,pml,n,modz,modx,tt,dt,ns,ds,rd,off,dj,sz,sx,mf,out_t)
character*40 input_filename
real dx,dz
integer pml,modx,modz,n
integer sx,sz,ns,ds,rd,off,dj
integer tt,nt,out_t
real mf,dt
input_filename = 'cmd.txt'
open(111,file=input_filename)
read(111,*)dx,dz
read(111,*)pml
read(111,*)n
read(111,*)modz,modx
read(111,*)tt,dt
read(111,*)ns,ds
read(111,*)rd,off,dj
read(111,*)sz,sx
read(111,*)mf
read(111,*)out_t
out_t=out_t/0.2
close(111)
end subroutine
subroutine calculate_jcwg(ins,pml,mx,mz,n,sx,sz,rd,off,dj,dt,nt,mf,dx,dz,v,st,j1w,out_w,out_t)
real dx,dz !节点长度
integer pml,n !吸收层厚度，空间阶数n*2
integer mx,mz !立方单元数
integer sx,sz !震源节点号
integer nt !总步数
integer rd,off,dj,out_t
real mf,dt !子波主频，采样间隔
real wavelet(nt) !地震子波值
real cn(n)
real vmax
integer ins,m
real v(mz,mx),st(mz,mx),dtdx
real U(mz,mx),W(mz,mx)
real p(mz,mx)
real decz(mz),decx(mx)
real pzw,pxu
real ppzw(mz,mx),ppxu(mz,mx)
real ppxp(mz,mx),ppzp(mz,mx)
real pxp,pzp
character*20 outfile,fft
character*4 filenumb
real d1,d2,a
real jlw(nt,rd),out_w(mz,mx) !切片输出
U=0.
W=0.
p=0.
wavelet=0.
ppzw=0
ppxu=0
ppxp=0
ppzp=0
call calcn(n,cn) !计算空间差分系数
call calwave(nt,wavelet,mf,dt) !计算地震子波
open(11,file='wavelet.dat')
do i=1,nt
write(11,*)wavelet(i)
end do
close(11)
call calpmldec(decz,decx,mx,mz,dx,n,pml,v) !计算pml衰减因子
if(maxval(v)*dt/1000/dx>sqrt(2.0)/2)then
write(*,*)'模型不满足稳定性条件'
pause
end if
a=sqrt(-1.0*log(0.01)/((10.0*dx)**2.0)); !空间衰减因子，其值越大则子波沿空间衰减越迅速
dtdx=dt/1000/dx
!开始计算
do it=1,nt
forall(i=n+1:mz-n,j=n+1:mx-n)
p(i,j) = p(i,j) + wavelet(it) * exp(-a**2 * (((i-sz)*dz)**2 + ((j-sx)*dx)**2))
end forall
do m=1,mx
p(n+1,m) = 0 ! 自由表面应力赋0
end do
do m=1,n
do k=1,mx
p(m,k) = -p(2*(n+1)-m,k) ! 应力关于自由表面上下反对称
end do
end do
do j=n+1,mx-n
do i=n+1,mz-n
if(st(i,j)==1)then !进行声波波动方程
!普通交错网格
pxp=0
do m=1,n
pxp=pxp+cn(m)*(p(i,j+m)-p(i,j-m+1)) !对x方向压力进行中心差分
end do
pzp=0
do m=1,n
pzp=pzp+cn(m)*(p(i+m,j)-p(i-m+1,j)) !对z方向压力进行中心差分
end do
ppxp(i,j)=ppxp(i,j)+(pxp-ppxp(i,j))*decx(j)*dt/1000
ppzp(i,j)=ppzp(i,j)+(pzp-ppzp(i,j))*decz(i)*dt/1000
u(i,j)=u(i,j)+(pxp-ppxp(i,j))*dtdx
w(i,j)=w(i,j)+(pzp-ppzp(i,j))*dtdx
!普通交错网格
!u
pxu=0
do m=1,n
pxu=pxu+cn(m)*(u(i,j+m-1)-u(i,j-m)) !对x方向速度进行中心差分
end do
!w
pzw=0
do m=1,n
pzw=pzw+cn(m)*(w(i+m-1,j)-w(i-m,j)) !对z方向速度进行中心差分
end do
ppxu(i,j)=ppxu(i,j)+(pxu-ppxu(i,j))*decx(j)*dt/1000
ppzw(i,j)=ppzw(i,j)+(pzw-ppzw(i,j))*decz(i)*dt/1000
p(i,j)=p(i,j)+(v(i,j)**2*(pxu-ppxu(i,j))+v(i,j)**2*(pzw-ppzw(i,j)))*dtdx !根据速度反演压力
else !进行弹性波波动方程
!普通交错网格
!p
pxp=0
do m=1,n
pxp=pxp+cn(m)*(p(i,j+m-1)-p(i,j-m)) !对x方向应力进行中心差分
end do
pzp=0
do m=1,n
pzp=pzp+cn(m)*(p(i+m,j)-p(i-m+1,j)) !对z方向应力进行中心差分
end do
ppxp(i,j)=ppxp(i,j)+(pxp-ppxp(i,j))*decx(j)*dt/1000
ppzp(i,j)=ppzp(i,j)+(pzp-ppzp(i,j))*decz(i)*dt/1000
u(i,j)=u(i,j)+(pxp-ppxp(i,j))*dtdx
w(i,j)=w(i,j)+(pzp-ppzp(i,j))*dtdx
!求取速度分量结束
!普通交错网格
!u
pxu=0
do m=1,n
pxu=pxu+cn(m)*(u(i,j+m)-u(i,j-m+1))
end do
!w
pzw=0
do m=1,n
pzw=pzw+cn(m)*(w(i+m-1,j)-w(i-m,j))
end do
ppzw(i,j)=ppzw(i,j)+(pzw-ppzw(i,j))*decz(i)*dt/1000
ppxu(i,j)=ppxu(i,j)+(pxu-ppxu(i,j))*decx(j)*dt/1000
p(i,j)=p(i,j)+(v(i,j)*v(i,j)*(pxu-ppxu(i,j))+v(i,j)*v(i,j)*(pzw-ppzw(i,j)))*dtdx
!求取应力分量结束
end if
end do
end do
do i=1,rd
!jlw(it,i)=w(sz,sx+(i-1)*dj-off)
jlw(it,i)=w(25,sx+(i-1)*dj-off) !从某一平面x方向输出
!起伏海底输出
!m=(i-1)*6+35
!if(m<335) then
! jlw(it,i)=w((1265-m)/6,m)
!else
! m=(i-1)*6+35
! jlw(it,i)=w((m+595)/6,m)
!end if
end do
if(it==out_t) then
out_w=w
!!输出某一偏移距的波形图
!open(333,file='250道.txt')
!do i=1,mz
!write(333,*) w(i,250)
!end do
!close(333)
else
end if
end do
open(555,file='20米深-50道.txt')
do i=1,nt
write(555,*) jlw(i,50) !输出某一道检波器单道记录
end do
close(555)
end subroutine
subroutine calcn(n,cn)
integer n,l
real cn(n)
real A(n,n),B(n)
do jc=1,n
do ic=1,n
A(ic,jc)=(2*jc-1)**(2*ic-1)
end do
end do
B=0.
B(1)=1.0
call AGAUS(A,B,n,cn,l)
if(l==0) then
write(*,*)'矩阵奇异，计算cn错误'
stop
end if
end subroutine
subroutine calwave(wl,wavelet,mf,dt) !计算地震子波
integer wl
real mf,pi,dt
real wavelet(wl)
real delay
pi=3.141592653
delay=1.5*sqrt(6.0)/(pi*mf)
do iw=1,wl
wavelet(iw)=(1-2.0*pi**2*mf**2*(1.0*(iw*dt/1000-delay))**2)* &
exp(-pi**2*mf**2*(1.0*(iw*dt/1000-delay))**2)
end do
end subroutine
subroutine out_sgy(m,n,ns,a,filename,dt)
integer m,n,ns
real dt
real a(ns,m,n)
character*20 filename
integer*2 seg_y_head(120)
integer*4 seg_y_head_4(60)
equivalence(seg_y_head,seg_y_head_4)
integer rec_long,i,j
seg_y_head(58)=m !采样点数
seg_y_head(59)=int(dt*1000) !采样间隔
rec_long=240+m*4 !计算道长
open(unit=13,file=filename,access='direct',status='replace',form='binary',recl=rec_long)
do is=1,ns
do i=1,n,1 !读记录到矩阵形式
write(13,rec=i+(is-1)*n)seg_y_head,a(is,:,i)
end do
end do
close(13)
end subroutine
subroutine calpmldec(decz,decx,mx,mz,dx,n,pml,v) !计算pml衰减系数
integer mx,mz,n,pml,x
real decz(mz),decx(mx),dx
real v(mz,mx)
real r,pi,c
pi=3.141592653
r=0.001
decz=0.
decx=0.
!do iz=n+1,pml+n+1
! x=n+pml+1-iz
! decz(iz)=v*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
!end do
do iz=mz-pml-n,mz-n
x=iz-(mz-pml-n)
c=v(iz,n+1)!选择z方向的速度作为衰减系数的计算
decz(iz)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=n+1,pml+n+1
x=n+pml+1-jx
c=v(n+1,jx)!选择x方向的速度作为衰减系数的计算
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=mx-pml-n,mx-n
x=jx-(mx-pml-n)
c=v(n+1,jx)!选择x方向的速度作为衰减系数的计算
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
end subroutine
SUBROUTINE AGAUS(A,B,N,X,L) !全选主元高斯消去法
real A(N,N),X(N),B(N),JS(N)
!DIMENSION A(N,N),X(N),B(N),JS(N)
!DOUBLE PRECISION A,B,X,T
L=1
DO 50 K=1,N-1
D=0.0
DO 210 I=K,N
DO 210 J=K,N
IF (ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
JS(K)=J
IS=I
END IF
210 CONTINUE
IF (D+1.0.EQ.1.0) THEN
L=0
ELSE
IF (JS(K).NE.K) THEN
DO 220 I=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
220 CONTINUE
END IF
IF (IS.NE.K) THEN
DO 230 J=K,N
T=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=T
230 CONTINUE
T=B(K)
B(K)=B(IS)
B(IS)=T
END IF
END IF
IF (L.EQ.0) THEN
WRITE(*,100)
RETURN
END IF
DO 10 J=K+1,N
A(K,J)=A(K,J)/A(K,K)
10 CONTINUE
B(K)=B(K)/A(K,K)
DO 30 I=K+1,N
DO 20 J=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
20 CONTINUE
B(I)=B(I)-A(I,K)*B(K)
30 CONTINUE
50 CONTINUE
IF (ABS(A(N,N))+1.0.EQ.1.0) THEN
L=0
WRITE(*,100)
RETURN
END IF
X(N)=B(N)/A(N,N)
DO 70 I=N-1,1,-1
T=0.0
DO 60 J=I+1,N
T=T+A(I,J)*X(J)
60 CONTINUE
X(I)=B(I)-T
70 CONTINUE
100 FORMAT(1X,' FAIL ')
JS(N)=N
DO 150 K=N,1,-1
IF (JS(K).NE.K) THEN
T=X(K)
X(K)=X(JS(K))
X(JS(K))=T
END IF
150 CONTINUE
RETURN
END SUBROUTINE
subroutine vr(v,st,mz,mx,n)
integer mz,mx,i,n,j
real v(mz,mx),st(mz,mx)
! 设定地层模型
do j=1,mx
do i=1,mz
if (i <= 200) then
v(i,j) = 1500
st(i,j) = 1
elseif (i > 200 .and. i <= 300) then
if (j >= 500) then
v(i,j) = 1500
st(i,j) = 1
elseif(j<=400) then
v(i,j)=2000
st(i,j)=2
elseif(j > (400+100.*(i-200)/100.) .and. j <= 500) then
v(i,j)=1500
st(i,j)=1
else
v(i,j)=2000
st(i,j)=2
end if
elseif(i > 300.and.i <= 400) then
v(i,j)=2000
st(i,j)=2
else
v(i,j)=2500
st(i,j)=2
end if
end do
end do
return
end subroutine
```c
从错误描述来看,问题可能出在:
数组下标越界。要检查所有数组访问是否在合理范围内,特别是循环中。
除数为0。在矩阵求逆等运算中要检查除数是否为0,以免出错。
内存不足。要检查程序是否申请了过多内存,导致内存溢出。
其他异常情况。要考虑所有可能的异常情况,做好异常处理。
从代码来看,有以下几点需要检查:
数组jlw, out_w的大小是否设置正确,是否会出现下标越界。
subroutine AGAUS中的除数是否会为0,如果是需要做判断。
申请的内存是否过大,导致内存溢出。可以检查内存使用情况。
是否考虑到了所有的异常情况,比如文件打开失败等。
总的来说,需要检查:
所有数组访问是否在合理范围内,没有下标越界。
所有的除法运算是否会出现除数为0的情况,如果有需要加判断。
程序申请的内存总量是否太大,内存使用情况。
所有的异常情况是否都有异常处理,文件操作等。
以上都是Fortran程序常见的错误,希望能对您有所帮助。如果问题 persists,欢迎提供更详细的错误信息,我们可以进一步分析。
提示winsig.c里面的418行报错，在418行的地方打一个断点，看看运行到此处的数据情况是否符合预期，判定原因
该回答引用ChatGPT
您提供的代码和说明来看,可能存在以下几个问题:
- 数组越界。在对数组进行索引访问时,需要确保索引在数组的有效范围内。如果访问了数组范围外的元素,会导致错误。
- 内存溢出。如果申请的数组太大,超出内存限制,也会导致错误。
- 死循环。如果存在某些循环条件导致循环无法正常终止,也会产生界面无响应的情况。
- 其他逻辑错误。除上述几点外,任何其他逻辑错误也可能导致程序运行异常。
针对上述可能的问题,我有以下几点建议: - 仔细检查所有数组的索引访问,确保索引在有效范围内。特别是在循环中,需要确保开始索引和结束索引正确。
- 检查程序申请内存的总量,确保不超过系统内存限制。如果申请的数组过大,可以适当调整数组的维度使内存消耗减小。
- 检查所有循环结构,确保循环能够正常终止。需要考虑所有可能的循环终止条件。
- 逐步调试程序,检查关键逻辑和算法是否正确。需要考虑各种边界条件和异常情况。
- 可以尝试减小问题的规模和复杂度,先实现一个简单的示例进行测试。等示例运行正常后,再逐步扩大规模。
- 可以在 subroutine 和 function 中插入 write 语句,打印关键变量的值,观察程序的运行过程,从而发现异常。
- 可以使用工具如 gdb 对程序进行调试,设置断点并检查变量的值。
代码
`fortran
subroutine fdtd2d()
implicit none
integer n !空间阶数
integer nt !总步数
integer mx,mz !立方单元数
integer sx,sz !震源节点号
real mf,dt !子波主频,采样间隔
real wavelet(nt) !地震子波
real cn(n)
real vmax
integer ins,m
real, allocatable :: v(:,:),st(:,:)
real U(:,:),W(:,:)
real p(:,:)
real decz(:),decx(:)
real dtdx
if (maxval(v)*dt/1000/dx>sqrt(2.0)/2)then
write(,)'模型不满足稳定性条件'
stop
end if
call calcn(n,cn) !计算空间差分系数
call calwave(nt,wavelet,mf,dt) !计算地震子波
allocate(v(mz,mx),st(mz,mx))
call vr(v,st,mz,mx,n) !设定模型
allocate(U(mz,mx),W(mz,mx),p(mz,mx))
call calpmldec(decz,decx,mx,mz,dx,n)
dtdx=dt/1000/dx
!开始计算
do it=1,nt
do i = 1, mz
do j = 1, mx
if (i <= mx .and. i > 0 .and. j <= mz .and. j > 0) then
p(i,j) = p(i,j) + wavelet(it) * exp(-((i-sz)*dz)**2/((ins*dx)**2) &
-((j-sx)*dx)**2/((ins*dx)**2))
else
write(,) '数组索引越界!'
end if
end do
end do
!其他语句...
!插入write语句,打印变量
if(it==nt/2) then
write(,) 'it=', it, 'p(25,25)=', p(25,25)
end if
!设置循环终止条件
do i=1,n
if (条件) exit
end do
end do
end subroutine
subroutine calcn(n,cn)
integer n,l
real cn(n)
real A(n,n),B(n)
do jc=1,n
do ic=1,n
A(ic,jc)=(2*jc-1)**(2*ic-1)
end do
end do
B=0.
B(1)=1.0
call AGAUS(A,B,n,cn,l)
if(l==0) then
write(,)'矩阵奇异,计算cn错误'
stop
end if
end subroutine
subroutine calwave(wl,wavelet,mf,dt) !计算地震子波
integer wl
real mf,pi,dt
real wavelet(wl)
real delay
pi=3.141592653
delay=1.5*sqrt(6.0)/(pi*mf)
do iw=1,wl
wavelet(iw)=(1-2.0*pi**2*mf**2*(1.0*(iw*dt/1000-delay))**2)* &
exp(-pi**2*mf**2*(1.0*(iw*dt/1000-delay))**2)
end do
end subroutine
subroutine out_sgy(m,n,ns,a,filename,dt)
integer m,n,ns
real dt
real a(ns,m,n)
character*20 filename
integer*2 seg_y_head(120)
integer*4 seg_y_head_4(60)
equivalence(seg_y_head,seg_y_head_4)
integer rec_long,i,j
seg_y_head(58)=m !采样点数
seg_y_head(59)=int(dt*1000) !采样间隔
rec_long=240+m*4 !计算道长
open(unit=13,file=filename,access='direct',status='replace', recl=rec_long)
do j=1,n
do i=1,ns
write(13,rec=i)seg_y_head,(a(i,k,j),k=1,m)
end do
end do
close(13)
end subroutine
subroutine calpmldec(decz,decx,mx,mz,dx,n,pml,v)
integer mx,mz,n,pml
real dx
real v(mz,mx)
real decz(mz),decx(mx)
real r,x,c
data r/0.00001/
do jz=1,n
x=jz
c=v(jz,1)
decz(jz)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jz=mz-n+1,mz
x=mz-jz+1
c=v(mz-n+1,1)
decz(jz)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=1,n
x=jx
c=v(1,jx)
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=n+1,pml+n+1
x=n+pml+1-jx
c=v(n+1,jx)!选择x方向的速度作为衰减系数的计算
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
do jx=mx-pml-n,mx-n
x=jx-(mx-pml-n)
c=v(n+1,jx)!选择x方向的速度作为衰减系数的计算
decx(jx)=c*3/2/pml/dx*log(1.0/r)*(1.0*x/pml)**2
end do
end subroutine
SUBROUTINE AGAUS(A,B,N,X,L)
real A(N,N),X(N),B(N),JS(N)
L=1
DO 50 K=1,N-1
D=0.0
DO 210 I=K,N
DO 210 J=K,N
IF (ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
JS(K)=J
IS=I
END IF
210 CONTINUE
IF (D+1.0.EQ.1.0) THEN
L=0
ELSE
IF (JS(K).NE.K) THEN
DO 220 I=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
220 CONTINUE
END IF
IF (IS.NE.K) THEN
DO 230 J=K,N
T=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=T
230 CONTINUE
T=B(K)
B(K)=B(IS)
B(IS)=T
END IF
END IF
IF (L.EQ.0) THEN
WRITE(*,100)
RETURN
END IF
DO 10 J=K+1,N
A(K,J)=A(K,J)/A(K,K)
10 CONTINUE
B(K)=B(K)/A(K,K)
DO 30 I=K+1,N
DO 20 J=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
20 CONTINUE
B(I)=B(I)-A(I,K)*B(K)
30 CONTINUE
50 CONTINUE
IF (ABS(A(N,N))+1.0.EQ.1.0) THEN
L=0
WRITE(*,100)
RETURN
END IF
X(N)=B(N)/A(N,N)
DO 70 I=N-1,1,-1 T=0.0
DO 60 J=I+1,N
T=T+A(I,J)*X(J)
60 CONTINUE
X(I)=(B(I)-T)/A(I,I)
70 CONTINUE
100 FORMAT(' 矩阵奇异,不能求逆')
RETURN
END SUBROUTINE
subroutine vr(v,st,mz,mx,n)
integer mz,mx,n,i,j
real v(mz,mx),st(mz,mx)
v=2000.
st=2 !初始条件为全体积的P波速度
do i = 1, mz
do j = 1, mx
if(i <= n) then
if(j > (200+50.*(i-100)/100.) .and. j <= 250) then
v(i,j)=1500
st(i,j)=1
elseif(j > (400+100.*(i-200)/100.) .and. j <= 500) then
v(i,j)=1500
st(i,j)=1
else
v(i,j)=2000
st(i,j)=2
end if
elseif(i > n .and. i <= 2*n) then
if(j > (200+50.*(i-100-n)/100.) .and. j <= 250) then
v(i,j)=1500
st(i,j)=1
elseif(j > (400+100.*(i-200-n)/100.) .and. j <= 500) then
v(i,j)=1500
st(i,j)=1
else
v(i,j)=2000
st(i,j)=2
end if
elseif(i > 2*n .and. i <= 300) then
if(j > (200+50.*((i-100-2*n)/100.)) .and. j <= 250) then
v(i,j)=1000
st(i,j)=1
elseif(j > (400+100.*((i-200-2*n)/100.)) .and. j <= 500) then
v(i,j)=1000
st(i,j)=1
else
v(i,j)=1500
st(i,j)=2
end if
elseif(i > 300.and.i <= 400) then
v(i,j)=2000
st(i,j)=2
else
v(i,j)=2500
st(i,j)=2
end if
end do
end do
end subroutine