📄 wwhmodel.f90
字号:
!!!!!!!!!!!!!输出快照的程序!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!输出快照的程序!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!参数设置
! nx-水平方向网格数,nz-深度方向网格数,dx-网格大小
! nmax-时间采样点数,dt-采样间隔,ixs-激发点水平方向网格号,izs-激发点深度方向网格号
! nsnap-快照输出最大炮数,ifsnap-快照输出时间起点步长,idsnap-快照输出时间增加步长
! isotype-震源类型:1/点震源,2/剪切震源,3/垂直震源,4/水平震源
! iabmax-吸收边界网格数,topbc-地面边界类型:free/自由边界,abbc/吸收边界
parameter(nx=600,nz=600,dx=0.2,nz0=5)
parameter(nmax=2001,dt=0.1e-3,ixs=300,izs=5)
parameter(nsnap=40,ifsnap=100,idsnap=100,isotype=1,iabmax=80)
parameter(topbc='free')
!!!!!!!!!!!!定义数组
dimension vx(0:nx,0:nz),vz(0:nx,0:nz)
dimension txx(0:nx,0:nz),tzz(0:nx,0:nz)
dimension txz(0:nx,0:nz),bx(0:nx,0:nz),bz(0:nx,0:nz)
dimension rlamu(0:nx,0:nz),rla(0:nx,0:nz),rmuxz(0:nx,0:nz)
dimension rho(0:nx,0:nz),source(nmax),eponge(0:nx)
dimension recx(0:nx,nmax),recz(0:nx,nmax),recxz(0:nx,nmax)
character*12 ffo11,ffo12,ffo21,ffo22,ffo23
data ffo21,ffo22,ffo23/'recx.grd','recz.grd','recxz.grd'/
!!!!!!!!!!!!!g1,g2时间2阶-空间4阶常数,soufac震源因子
g1=-1.0/24.0
g2=9.0/8.0
dtdx=dt/dx
soufac=dt/(dx*dx)
isnap=0
!!!!!!!!!!!!!生成雷克子波
call wavelet(nmax,dt,source)
print*,'wavelet'
amax=100.0*amax/(0.25*soufac)
!!!!!!!!!!!!!生成地质模型
call model(rlamu,rla,rmuxz,rho,bx,bz,dtdx,nx,nz,bx0,bz0,bz1,rmuxz1,rla0,rla1,rlamu0,rlamu1)
print*,'model'
!!!!!!!!!!!!!生成吸收边界系数
xx=0.305/float(iabmax)
do ib=0,iabmax
aaa=(xx*float(iabmax-ib))**2
eponge(ib)=exp(-aaa)
enddo
!!!!!!!!!!!!!数组初始化
do i=0,nx
do k=0,nz
vx(i,k)=0.0
vz(i,k)=0.0
txx(i,k)=0.0
tzz(i,k)=0.0
txz(i,k)=0.0
enddo
enddo
!!!!!!!!!!!!!获得质点速度(4阶)
do n=1,nmax
do k=2,nz-2
do i=2,nx-2
vx(i,k)=vx(i,k)+bx(i,k)*(g1*(txx(i+1,k)-txx(i-2,k))+g2*(txx(i,k)-txx(i-1,k))+g1*(txz(i,k+2)-txz(i,k-1))+g2*(txz(i,k+1)-txz(i,k)))
vz(i,k)=vz(i,k)+bz(i,k)*(g1*(txz(i+2,k)-txz(i-1,k))+g2*(txz(i+1,k)-txz(i,k))+g1*(tzz(i,k+1)-tzz(i,k-2))+g2*(tzz(i,k)-tzz(i,k-1)))
enddo
enddo
!!!!!!!!!!!!单炮记录赋值
do ix=0,nx
recx(ix,n)=vx(ix,nz0)*1000000000
recz(ix,n)=vz(ix,nz0)*1000000000
recxz(ix,n)=(vx(ix,nz0)+vz(ix,nz0))*10000000/2.0
enddo
!!!!!!!!!!!!!零碎东西
do k=1,1
do i=2,nx-2
vx(i,k)=vx(i,k)+bx(i,k)*(g1*(txx(i+1,k)-txx(i-2,k))+g2*(txx(i,k)-txx(i-1,k))+g1*(txz(i,k+2)-txz(i,k-1))+g2*(txz(i,k+1)-txz(i,k)))
enddo
enddo
do k=2,nz-2
do i=1,1
vz(i,k)=vz(i,k)+bz(i,k)*(g1*(txz(i+2,k)-txz(i-1,k))+g2*(txz(i+1,k)-txz(i,k))+g1*(tzz(i,k+1)-tzz(i,k-2))+g2*(tzz(i,k)-tzz(i,k-1)))
enddo
enddo
if(mod(n,100).eq.0)print*,'timestep=',n
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!速度边界条件
!!!!!!!!!!!!!左右边界吸收条件
do ka=0,nz
do ia=0,iabmax
vx(ia,ka)=vx(ia,ka)*eponge(ia)
vz(ia,ka)=vz(ia,ka)*eponge(ia)
vx(nx-ia,ka)=vx(nx-ia,ka)*eponge(ia)
vz(nx-1-ia,ka)=vz(nx-1-ia,ka)*eponge(ia)
enddo
enddo
!!!!!!!!!!!!!上下边界吸收条件
do ka=0,iabmax
do ia=0,nx
if(topbc.eq.'abbc')then
vx(ia,ka)=vx(ia,ka)*eponge(ka)
vz(ia,ka)=vz(ia,ka)*eponge(ka)
endif
vx(ia,nz-1-ka)=vx(ia,nz-1-ka)*eponge(ka)
vz(ia,nz-ka)=vz(ia,nz-ka)*eponge(ka)
enddo
enddo
!!!!!!!!!!!!!自由表面边界条件
if(topbc.eq.'free')then
do i=2,nx-2
vx(i,0)=vx(i,0)+bx0*(g1*(txx(i+1,0)-txx(i-2,0))+g2*(txx(i,0)-txx(i-1,0))+g1*(txz(i,2)+txz(i,1))+g2*(txz(i,1)-txz(i,0)))
enddo
do i=1,nx-2
vz(i,0)=vz(i,0)+2*bz0*(g1*tzz(i,1)+g2*tzz(i,0))
vz(i,1)=vz(i,1)+bz1*(g1*(tzz(i,2)+tzz(i,0))+g2*(tzz(i,1)-tzz(i,0))+g1*(txz(i+2,1)-txz(i-1,1))+g2*(txz(i+1,1)-txz(i,1)))
enddo
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!均匀弹性介质中应力-应变之间的关系
do k=2,nz-2
do i=2,nx-2
exx=g1*(vx(i+2,k)-vx(i-1,k))+g2*(vx(i+1,k)-vx(i,k))
ezz=g1*(vz(i,k+2)-vz(i,k-1))+g2*(vz(i,k+1)-vz(i,k))
if(k.eq.(nz-2))then
exz=vx(i,k)-vx(i,k-1)+g1*(vz(i+1,k)-vz(i-2,k))+g2*(vz(i,k)-vz(i-1,k))
if(i.eq.1)exz=vx(1,k)-vx(1,k-1)+vz(1,k)-vz(0,k)
endif
if(k.ge.2.and.k.le.(nz-3))then
exz=g1*(vx(i,k+1)-vx(i,k-2))+g2*(vx(i,k)-vx(i,k-1))+g1*(vz(i+1,k)-vz(i-2,k))+g2*(vz(i,k)-vz(i-1,k))
if(i.eq.1)then
exz=vz(1,k)-vz(0,k)+g1*(vx(1,k+1)-vx(i,k-2))+g2*(vz(1,k)-vz(i,k-1))
endif
endif
txx(i,k)=txx(i,k)+rlamu(i,k)*exx+rla(i,k)*ezz
tzz(i,k)=tzz(i,k)+rla(i,k)*exx+rlamu(i,k)*ezz
txz(i,k)=txz(i,k)+rmuxz(i,k)*exz
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!应力应变边界条件
!!!!!!!!!!!!!自由边界条件
if(topbc.eq.'free')then
do i=2,nx-2
txz(i,0)=0.0
txz(i,1)=txz(i,1)+g1*(vz(i+1,1)-vz(i-2,1))+g2*(vz(i,1)-vz(i-1,1))+(g1*(vx(i,2)+vz(i-1,0)-vz(i,0)-vx(i,0))+g2*(vx(i,1)-vx(i,0)))*rmuxz1
enddo
do i=1,nx-2
tzz(i,0)=tzz(i,0)+rlamu0*(vz(i,1)-vz(i,0))+rla0*(g1*(vx(i+2,0)-vx(i-1,0))+g2*(vx(i+1,0)-vx(i,0)))
tzz(i,1)=tzz(i,1)+rlamu1*(g1*(vz(i,3)-vz(i,0))+g2*(vz(i,2)-vz(i,1)))+rla1*(g1*(vx(i+2,1)-vx(i-1,1))+g2*(vx(i+1,1)-vx(i,1)))
enddo
do i=1,nx-2
txx(i,0)=txx(i,0)+rla0*(vz(i,1)-vz(i,0))+rlamu0*(g1*(vx(i+2,0)-vx(i-1,0))+g2*(vx(i+1,0)-vx(i,0)))
txx(i,1)=txx(i,1)+rla1*(g1*(vz(i,3)-vz(i,0))+g2*(vz(i,2)-vz(i,1)))+rlamu1*(g1*(vx(i+2,1)-vx(i-1,1))+g2*(vx(i+1,1)-vx(i,1)))
enddo
endif
!!!!!!!!!!!!!左,右吸收边界条件
do ka=0,nz
do ia=0,iabmax
txz(ia,ka)=txz(ia,ka)*eponge(ia)
txx(ia,ka)=txx(ia,ka)*eponge(ia)
tzz(ia,ka)=tzz(ia,ka)*eponge(ia)
txz(nx-ia,ka)=txz(nx-ia,ka)*eponge(ia)
txx(nx-1-ia,ka)=txx(nx-1-ia,ka)*eponge(ia)
tzz(nx-1-ia,ka)=tzz(nx-1-ia,ka)*eponge(ia)
enddo
enddo
!!!!!!!!!!!!!上,下吸收边界条件
do ka=0,iabmax
do ia=0,nx
if(topbc.eq.'abbc')then
txz(ia,ka)=txz(ia,ka)*eponge(ka)
txx(ia,ka)=txx(ia,ka)*eponge(ka)
tzz(ia,ka)=tzz(ia,ka)*eponge(ka)
endif
txz(ia,nz-ka)=txz(ia,nz-ka)*eponge(ka)
txx(ia,nz-1-ka)=txx(ia,nz-1-ka)*eponge(ka)
tzz(ia,nz-1-ka)=tzz(ia,nz-1-ka)*eponge(ka)
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!震源函数
!!!!!!!!!点震源
if(isotype.eq.1)then
addsou=soufac*source(n)*0.25
txx(ixs-1,izs)=txx(ixs-1,izs)+addsou
txx(ixs,izs)=txx(ixs,izs)+addsou
txx(ixs-1,izs+1)=txx(ixs-1,izs+1)+addsou
txx(ixs,izs+1)=txx(ixs,izs+1)+addsou
tzz(ixs-1,izs)=tzz(ixs-1,izs)+addsou
tzz(ixs,izs)=tzz(ixs,izs)+addsou
tzz(ixs-1,izs+1)=tzz(ixs-1,izs+1)+addsou
tzz(ixs,izs+1)=tzz(ixs,izs+1)+addsou
endif
!!!!!!!!!!!!!剪切震源
if(isotype.eq.2)then
addsou=soufac*source(n)*0.25
vx(ixs-1,izs)=vx(ixs-1,izs)+addsou
vx(ixs,izs)=vx(ixs,izs)+addsou
vx(ixs-1,izs+1)=vx(ixs-1,izs+1)+addsou
vx(ixs,izs+1)=vx(ixs,izs+1)+addsou
vz(ixs-1,izs)=vz(ixs-1,izs)-addsou
vz(ixs,izs)=vz(ixs,izs)-addsou
vz(ixs-1,izs+1)=vz(ixs-1,izs+1)-addsou
vz(ixs,izs+1)=vz(ixs,izs+1)-addsou
endif
!!!!!!!!!!!!!垂直震源
if(isotype==3)then
addsou=0.5*soufac*source(n)/rho(ixs-1,izs+1)
vz(ixs-1,izs+1)=vz(ixs-1,izs+1)+addsou
vz(ixs,izs+1)=vz(ixs,izs+1)+addsou
endif
!!!!!!!!!!!!!水平震源
if(isotype.eq.4)then
addsou=0.5*soufac*source(n)/rho(ixs,izs)
vx(ixs,izs)=vx(ixs,izs)+addsou
vx(ixs,izs+1)=vx(ixs,izs+1)+addsou
endif
!!!!!!!!!!!!!保存快照数据
! if(mod(n,idsnap).eq.0.and.isnap.lt.nsnap)then
! ffo11=' '
! ffo12=' '
! ffo22=' '
! ffo11='kz'//char(isnap/10+48)//char(mod(isnap,10)+48)//'x.grd'
! ffo12='kz'//char(isnap/10+48)//char(mod(isnap,10)+48)//'z.grd'
! ffo22='kz'//char(isnap/10+48)//char(mod(isnap,10)+48)//'x+z.grd'
! open(910,file=ffo11)
! write(910,'(a4)') 'DSAA'
! write(910,'(2i5)') int((nx)/2+1),int((nz)/2+1)
! write(910,'(2i5)') 0,int(nx*dx)
! write(910,'(2i5)') -int(nz*dx),0
! write(910,'(2i5)') -1000,1000
!
! open(920,file=ffo12)
! write(920,'(a4)') 'DSAA'
! write(920,'(2i5)') int((nx)/2+1),int((nz)/2+1)
! write(920,'(2i5)') 0,int(nx*dx)
! write(920,'(2i5)') -int(nz*dx),0
! write(920,'(2i5)') -1000,1000
!
! open(930,file=ffo22)
! write(930,'(a4)') 'DSAA'
! write(930,'(2i5)') int((nx)/2+1),int((nz)/2+1)
! write(930,'(2i5)') 0,int(nx*dx)
! write(930,'(2i5)') -int(nz*dx),0
! write(930,'(2i5)') -1000,1000
!
! do iz=nz,0,-2
! write(910,'(10g)')(vx(ix,iz)*1000000000,ix=0,nx,2)
! write(920,'(10g)')(vz(ix,iz)*1000000000,ix=0,nx,2)
! write(930,'(10g)')((vz(ix,iz)+vx(ix,iz))*1000000000,ix=0,nx,2)
! enddo
! close(910)
! close(920)
! close(930)
! isnap=isnap+1
! endif
enddo
!!!!!!!!!!!!保存单炮记录数据
open(930,file='x.grd')
write(930,'(a4)') 'DSAA'
write(930,'(2i6)') nx+1,nmax
write(930,'(2i6)') 0,int(nx*dx)
write(930,'(2i6)') -int(1000*(nmax-1)*dt),0
write(930,'(2f5.1)') -0.1,0.1
open(940,file='z.grd')
write(940,'(a4)') 'DSAA'
write(940,'(2i6)') nx+1,nmax
write(940,'(2i6)') 0,int(nx*dx)
write(940,'(2i6)') -int(1000*(nmax-1)*dt),0
write(940,'(2f5.1)') -0.1,0.1
open(950,file='x-z.grd')
write(950,'(a4)') 'DSAA'
write(950,'(2i6)') nx+1,nmax
write(950,'(2i6)') 0,int(nx*dx)
write(950,'(2i6)') -int(1000*(nmax-1)*dt),0
write(950,'(2f5.1)') -0.1,0.1
do n=nmax,1,-1
write(930,'(601g14.3)') (recx(ix,n),ix=0,nx)
write(940,'(601g14.3)') (recz(ix,n),ix=0,nx)
write(950,'(601g14.3)') (recxz(ix,n),ix=0,nx)
enddo
close(930)
close(940)
close(950)
!!输出自己需要的部分记录！！！！！！！！！！！！！！！！！！！！！！
open(960,file='x.dat')
open(970,file='z.dat')
open(980,file='x-z.dat')
do n=1,2001
write(960,1000) (recx(ix,n),IX=int(nx/2)+50,int(nx/2)+150,2)
write(970,1000) (recz(ix,n),IX=int(nx/2),int(nx/2)+300)
write(980,1000) (recxz(ix,n),IX=int(nx/2)+50,int(nx/2)+150,2)
ENDDO
1000 FORMAT(1X,301(F9.1,X))
close(960)
close(970)
close(980)
stop
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!模型参数设置
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine model(rlamu,rmu,rmuxz,rho,bx,bz,dtdx,nx,nz,bx0,bz0,bz1,rmuxz1,rmu0,rmu1,rlamu0,rlamu1)
parameter(nmat=3,eps=1.e-6)
dimension rlamu(0:nx,0:nz),rmu(0:nx,0:nz)
dimension rho(0:nx,0:nz),rmuxz(0:nx,0:nz)
dimension bx(0:nx,0:nz),bz(0:nx,0:nz)
dimension alpha(nmat),beta(nmat),dens(nmat)
data alpha(1),beta(1),dens(1)/1000.0,500.0,1800.0/
data alpha(2),beta(2),dens(2)/2200.0,1100.0,2650.0/
data alpha(3),beta(3),dens(3)/1500.0,0.0,1000.0/
Real Ra,Rx,Ry
Rx=300
Ry=260
write(*,*) nx,nz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!第一层
do ix=0,nx
do iz=0,nz
rlamu(ix,iz)=alpha(1)*alpha(1)*dens(1)
rmu(ix,iz)=beta(1)*beta(1)*dens(1)
rho(ix,iz)=dens(1)
enddo
enddo
! do ix=0,nx
! do iz=101,nz
! rlamu(ix,iz)=alpha(2)*alpha(2)*dens(2)
! rmu(ix,iz)=beta(2)*beta(2)*dens(2)
! rho(ix,iz)=dens(2)
! enddo
! enddo
!!!!!!!!!!!!!结束模型构建
do i=0,nx
do k=0,nz
im=max0(0,i-1)
km=max0(0,k-1)
ii=min0(i,nx-1)
kk=min0(k,nz-1)
!!!!!!!!bx
bx(i,k)=0.5*dtdx*((1.0/rho(im,kk))+(1.0/rho(ii,kk)))
!!!!!!!!bz
bz(i,k)=0.5*dtdx*((1.0/rho(ii,km))+(1.0/rho(ii,kk)))
!!!!!!!!muxz
akan=1.0/(rmu(ii,kk)+eps)
akan=akan+(1.0/(rmu(im,kk)+eps))
akan=akan+(1.0/(rmu(ii,km)+eps))
akan=akan+(1.0/(rmu(im,km)+eps))
rmuxz(i,k)=(dtdx*4.0)/akan
if(rmuxz(i,k).lt.1.0)rmuxz(i,k)=0.0
enddo
enddo
do i=0,nx-1
do k=0,nz-1
!!!!!!!!lambda
rmu(i,k)=dtdx*(rlamu(i,k)-(2.0*rmu(i,k)))
!!!!!!!!lambda+2mu
rlamu(i,k)=dtdx*rlamu(i,k)
enddo
enddo
!!!!!!!!自由表面fd系数
bx0=bx(0,0)
bz0=bz(0,0)
bz1=bz(1,0)
rlamu0=rlamu(0,0)
rlamu1=rlamu(1,0)
rmu0=rmu(0,0)
rmu1=rmu(1,0)
rmuxz0=rmuxz(0,0)
!!!!!!!!!!!增加一些条件用来减少“IF”代码
!!!!!!!!!!!vx
do i=0,nx
bx(i,0)=0.0
bx(i,nz-1)=0.0
bx(i,nz)=0.0
enddo
do k=0,nz
bx(0,k)=0.0
bx(1,k)=0.0
bx(nx-1,k)=0.0
bx(nx,k)=0.0
enddo
!!!!!!!!!!!vz
do i=0,nx
bz(i,0)=0.0
bz(i,nz-1)=0.0
bz(i,nz)=0.0
enddo
do k=0,nz
bz(0,k)=0.0
bz(1,k)=0.0
bz(nx-1,k)=0.0
bz(nx,k)=0.0
enddo
!!!!!!!!!!!txx 和 tzz
do i=0,nx
rlamu(i,0)=0.0
rlamu(i,1)=0.0
rlamu(i,nz-2)=0.0
rlamu(i,nz-1)=0.0
rlamu(i,nz)=0.0
rmu(i,0)=0.0
rmu(i,1)=0.0
rmu(i,nz-2)=0.0
rmu(i,nz-1)=0.0
rmu(i,nz)=0.0
enddo
do k=0,nz
rlamu(0,k)=0.0
rlamu(nx-1,k)=0.0
rlamu(nx,k)=0.0
rmu(0,k)=0.0
rmu(nx-1,k)=0.0
rmu(nx,k)=0.0
enddo
!!!!!!!!!!!txz
do i=0,nx
rmuxz(i,0)=0.0
rmuxz(i,1)=0.0
rmuxz(i,nz-1)=0.0
rmuxz(i,nz)=0.0
enddo
do k=0,nz
rmuxz(0,k)=0.0
rmuxz(nx,k)=0.0
enddo
return
end
!!!!!!!!!!!!!sfreq-子波主频,tl-周期数
subroutine wavelet(nt,dt,wav)
dimension wav(*)
do i=1,nt
wav(i)=0.0
enddo
sfreq=80
lwav=1+int(1.55/(sfreq*dt))
do i=1,lwav
wav(i)=ricker(float(i-1)*dt,sfreq)
enddo
return
end
function ricker(time,sfreq)
pix=4.05367
g0=0.5
g1=0.565
g2=0.105
t0=time*sfreq*pix
ricker=0.0
ricker=g0*cos(t0)-g1*cos(2*t0)+g2*cos(3*t0)
return
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -