📄 2d_elastic.f90
字号:
program Forword_Elastic
!定义变量
! implicit none
integer*2 nx,nz,nt,np,I,J,K,T,f,krec,XS,ZS,line,&
rightbdr,leftbdr,lowbdr,topbdr,nLayer,zbdt,nrec,kstep,snapt,&
I0,J0,J00,bdrx,bdrz,kdot,kall,tall,kspace,kfirst,&
Dep(10),vmax,vmin,temp
integer :: time(0:1)
real Xmax,Zmax,dx,dz,pi,fq,dt,dtt,minute,second,C,timediff,Lo(10),V1(10),&
pa3,wmax,umax,tx2,tz2,dtime,dxspace,const1,const2
character*4 fname
real ,allocatable :: U0(:,:),U1(:,:),U2(:,:),W(:,:),WW(:,:)
real ,allocatable :: Vel(:,:),fv(:),Abc(:,:),Den(:,:),Imp(:,:)
call system_clock(time(0))
write(*,*)'计算开始，请耐心等待！'
write(*,*)time(0)
pi=3.1415927
!读取模型及其网格化数据
OPEN (1,FILE='E:\2D_Elastic_Modeling\Model.par')
read (1,*)Xmax !模型x方向的大小,单位m
read (1,*)Zmax !模型z方向的大小,单位m
read (1,*)dx !x方向采样间隔,单位m
read (1,*)dz !z方向采样间隔,单位m
read (1,*)dtt !时间采样间隔,单位ms
read (1,*)nt !采样点数
read (1,*)XS,ZS !震源的坐标x,z,单位是m
read (1,*)line !测线位置(深度)
read (1,*)f !主频
read (1,*)nLayer !层数
!模型参数
read (1,*)(Dep(I),I=1,nLayer) !各层的网格数,即各层厚度
read (1,*)(V1(I),I=1,nlayer) !各层速度
read (1,*)(Lo(I),I=1,nlayer) !各层密度
!数据保存参数
read (1,*)nrec !保存文件道数
read (1,*)krec,kstep,kdot !记录道数的起始点,跳道数,跳点数--用于剖面输出
read (1,*)snapt !快照间隔
close (1)
dt=dtt/1000 !单位化为s
fq=f
zbdt=int(1/(4*f*dt)) !"zbdt"时间延迟=延迟1/4个周期
nx=Xmax/dx+1 !x方向网格点数
nz=Zmax/dz+1 !z方向网格点数
I0 = int(XS/dx+1) !震源在x,z方向的网格点号
J0 = int(ZS/dz+1) !I0,J0是否要+1??
J00= int(line/dz+1) !测线在z方向的网格点号
bdrx=50 !x,z方向吸收边界宽度
bdrz=50
rightbdr=nx+bdrx
leftbdr=-bdrx
topbdr=-bdrz
lowbdr=nz+bdrz
tx2=dt*dt/(dx*dx) !tx,tz:临时变量,用于后面时间和空间差分计算
tz2=dt*dt/(dz*dz)
kspace=dx*kstep !kstep=1,相邻道的间隔(单位为m)
kfirst=I0+(krec-2) !krec=1,震源前的第一个点??
! if ((krec-I0)+(nrec-1)*kstep>160) then
! kall=(160-krec+I0)/kstep
! else
kall=nrec !地震剖面的记录道数
! endif
tall=nt/kdot !nt=4000,kdot=1,地震剖面的采样点数
!屏幕显示data.dat中的内容
write(*,'( "xmodel,zmodel=",2F8.1)')Xmax,Zmax
write(*,'( "dx,dz=",2F6.2)')dx,dz
write(*,'( "dt=",F9.5)')dt
write(*,'( "nt=",I6)')nt
write(*,'( "xsource,zsource=",2I3)')XS,ZS
write(*,'( "krce,kstep,kdot=",5x,3I3)')krec,kstep,kdot
write(*,'( "f=",I2)')f
write(*,'( "I0,J0=",I4,5x,I4)')I0,J0
write(*,'( "J00=",I3)')J00
write(*,'( "snapt",I4)')snapt
write(*,'( "nlayer=",I1)')nLayer
write(*,*)(Dep(I),I=1,nLayer)
write(*,*)(V1(I),I=1,nLayer)
write(*,*)(Lo(I),I=1,nLayer)
!数组定界
allocate (U0(leftbdr:rightbdr,topbdr:lowbdr),U1(leftbdr:rightbdr,topbdr:lowbdr),U2(leftbdr:rightbdr,topbdr:lowbdr))
allocate (Vel(leftbdr:rightbdr,topbdr:lowbdr),Den(leftbdr: rightbdr,topbdr: lowbdr),&
W(0:nx,0:nt),WW(0:nx,0:nt),fv(0:nt),Abc(leftbdr:rightbdr,topbdr:lowbdr),&
Imp(leftbdr:rightbdr,topbdr:lowbdr))
!U0[-30:231,0:121]-,U1[,]-,U2[,]-,W-单炮剖面,WW-抽道后的单炮剖面,
!Den-density,Imp-impedance,Vel-velocity,Dep-depth
!Abc[,]-Absorbing Boundary Condition,fv[,]-子波,
!模型参数初始化
Do J=-bdrz,Dep(1) !第一层
Do I=-bdrx,nx+bdrx
Den(I,J)=Lo(1)
Vel(I,J)=V1(1)
Imp(I,J)=0
end do
end do
const1=Dep(1)
const2=const1+Dep(2)
DO K=2,nLayer-1 !中间各层
do J=const1+1,const2
do I=-bdrx,nx+bdrx
Den(I,J)=Lo(K)
Vel(I,J)=V1(K)
Imp(I,J)=0
end do
end do
const1=const2
const2=const1+Dep(K+1)
end do
do J=const1+1,nz+bdrz !最后一层
do I=-bdrx,nx+bdrx
Den(I,J)=Lo(nLayer)
Vel(I,J)=V1(nLayer)
Imp(I,J)=0
end do
end do
do J=0,nt
fv(J)=0 !"fv"为子波
do I=0,nx
W(I,J)=0 !"W"为单炮记录剖面
end do
end do
!波场初始化
do J=-bdrz,(nz+bdrz)
do I=-bdrx,nx+bdrx
U0(I,J)=0
U1(I,J)=0
U2(I,J)=0
end do
end do
!波阻抗计算
Do I=0,nx
Do J=0,nz
Imp(I,J)= (Vel(I,J)*Den(I,J)-Vel(I,J+1)*Den(I,J+1))/&
(Vel(I,J+1)*Den(I,J+1)+Vel(I,J)*Den(I,J))
END DO
END DO
!稳定性判定
write(*,*)'网格稳定性判断开始'
! do I=1,nlayer
pa3=V1(nlayer)*dtt/dx
dtime=0.5*(dx/V1(nlayer)) !一个网格传播的时间的1/2(按最大速度算)
dxspace=(V1(1)/(30*f)) !一个波长的(1/30--经验值)(按最小速度算)
! write(*,*) pa3,dt,dtime,dx,dxspace
if(pa3>sqrt(3.0/8.0)) then !.or. (pa3>=1.and.dt>=dtime).or.(dx>=dxspace)
write(*,*)'网格划分不稳定??'
! pause
! stop
endif
! end do
write(*,*)'网格划分稳定'
!吸收边界条件系数
C=0.00025
! c1=0.015
! AA=(0.015*0.8)**2/(bdrz**2)
DO J=-bdrz,nz+bdrz
DO I=-bdrx,nx+bdrx
if (I<=0.or.J<=0) then
! Abc(I,J)=exp(-AA*((bdrx+I)**4))
! Abc(I,J)=exp(c1*(nz-J))
! Abc(I,J)=1-exp(-c1*(bdrx+I)) !吸收小
Abc(I,J)=exp(C*(I)) !吸收大
else if(I>=nx.or.J>=nz) then
Abc(I,J)=Abc(nx-I,J)
else
Abc(I,J)=1.0 !无吸收
end if
end do
end do
!子波函数
! if(fv(k)>1.0E-5) then
! tt=(k-3*zbdt)*dt
! fv(k)=tt*exp(-a*tt*tt) !徐－高斯函数
! tt=pi*f*(k*dt-1/(1.5*fq)); !带频率的高斯函数
! fv(k)=exp(-2*tt**2);
! a2=pi*f*(k*dt-1/f) !雷克子波
! fv(k)=(1-2*a2*a2)*exp(-(a2*a2))
Do I=0,8*zbdt-1 !"k<8*zbdt"两个周期
tt=pi*f*(I+zbdt)*dt !向后延迟1/4个周期
fv(I)=(Vel(I0,J0)*Vel(I0,J0)*dt*dt/(dx*dz))*tt*exp(-2*tt**2) !带频率的高斯一阶导数!
end do
! write (*,*)(fv(I),I=0,8*zbdt-1)
T=0
!波动方程正演
! do k=-3*zbdt,nt
do K=0,nt
! IF (K<8*zbdt) THEN
! DO I=-bdrx+2,nx+bdrx-2 !界面同时爆炸
! DO J=-bdrz+2,nz+bdrz-2
! U1(I,J)=U1(I,J)+Imp(I,J)*fv(K)
! END DO
! END DO
! END IF
DO I=-bdrx+2,nx+bdrx-2
DO J=-bdrz+2,nz+bdrz-2
const1= Vel(I,J)**2*tx2
const2= Vel(I,J)**2*tz2
! U2(I,J)= (2-2.5*const1-2.5*const2)*U1(I,J)-U0(I,J)+&
! (const1/12)*(16*(U1(I+1,J)+U1(I-1,J)+U1(I,J+1)+U1(I,J-1))-&
! (U1(I+2,J)+U1(I-2,J)+U1(I,J+2)+U1(I,J-2)))
U2(I,J)= (2-2.5*const1-2.5*const2)*U1(I,J)-U0(I,J)+&
0.5*const1*(U1(I+1,J)+U1(I-1,J))+0.5*const2*(U1(I,J+1)+U1(I,J-1))
IF (I==I0.and.J==J0) THEN !加点震源力
U2(I,J)= U2(I,J)+fv(K)
endif
end do
End do
!检波器接收和边界吸收
Do I=1,nx-1
W(I,K)=U2(I,J00) !单炮记录输出W(1:nx-1,0:nt)
End do
DO J=-bdrz,nz+bdrz
DO I=-bdrx,nx+bdrx
U2(I,J)=U2(I,J)*Abc(I,J)
End do
End do
!快照输出
IF (K==T) then
if (T<10) then
WRITE(fname,'(I1)')T
else if (T>=10.and.T<100) then
WRITE(fname,'(I2)')T
else if (T>=100.and.T<1000) then
WRITE(fname,'(I3)')T
else if (T>=1000.and.T<10000) then
WRITE(fname,'(I4)')T
end if
OPEN (2,FILE='E:\2D_Elastic_Modeling\result\Usnapt'//fname//'.grd') !
WRITE(2,'(A4)') 'DSAA'
WRITE(2,'(2I6)') nx+1,nz+1
WRITE(2,'(2I6)') 0.0,nx
WRITE(2,'(2I6)') 0.0,nz
umax =0.0
do J=0,nz
do I=0,nx !数据归一化
if (abs(U2(I,J))>umax ) then
umax=abs(U2(I,J))
end if
end do
end do
WRITE(2,'(2E12.4)') -umax/10,umax/10
DO J=0,nz
write(2,'(9E12.4)') (U2(I,J),I=0,nx) !快照输出U2(0:nx,0,nz)
end do
close (2)
! if(T==1400) T=3100
T=T+snapt
write(*,*)K
end if
DO J=-bdrz,nz+bdrz
DO I=-bdrx,nx+bdrx
U0(I,J)=U1(I,J)
U1(I,J)=U2(I,J)
end do
end do
end do
!结果输出
OPEN (3,FILE='E:\2D_Elastic_Modeling\result\outz.grd')
WRITE(3,'(A4)') 'DSAA'
WRITE(3,'(2I6)') nx-1,nt+1 !纵横向坐标点数
WRITE(3,'(2I6)') 0.0,nx-2 !剖面横向长度
WRITE(3,*) 0.0,-nt*dt*1000 !剖面纵向长度(此处dt可当成每格的长度)
wmax=0.0
do J=0,nt
do I=1,nx-1
if (abs(W(I,J))>wmax ) then
wmax=abs(W(I,J))
end if
end do
end do
WRITE(3,*) -wmax/100,wmax/100
do J=0,nt
write(3,'(10E12.4)')(W(I,J),I=1,nx-1) !完整输出－以彩色影像显示
end do
close (3)
do J=1,tall !数据抽道
do I=1,kall
WW(I,J)=W((I-1)*kstep+krec,(J-1)*kdot) !???
end do
end do
! do J=zbdt,nt+zbdt,kdot
! write(15,'(10E12.4)')(Wx(I,J),I=krec,krec+kall,kstep) !抽道后以彩色影像显示
! write(16,'(10E12.4)')(Wz(I,J),I=krec,krec+kall*kstep,kstep)
! end do
OPEN (4,FILE='E:\2D_Elastic_Modeling\result\profig.grd')
WRITE(4,'(A4)') 'DSAA'
WRITE(4,'(2I6)') kall,tall
WRITE(4,'(2I6)') 0.0,kall-1
WRITE(4,'(2I6)') 0.0,tall-1
WRITE(4,*) -wmax/100,wmax/100
do J=1,tall
write(4,'(10E12.4)')(WW(I,J),I=1,kall)
end do
close (4)
!地震记录数据输出－sgy
! call write_sgy(Wzz,kall,tall,dt*kdot,kspace,kfirst)
call system_clock(time(1))
timediff = time(1)-time(0)
write(*,*)'恭喜你，计算完成！计算所有时间如下：单位：S'
minute=int(timediff/600000)
timed=timediff/600000
second=int((timed-minute)*60)
write(*,*)time(1)
! write(*,'(I4,A1,I2,A1)')minute,'m',second,'s'
write(*,*)minute,'m',second,'s'
write(*,*)int(timediff/10000)
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -