Hi All,
I want to write a code that shows the numerical modeling of surface waves over shallow cavities (underground cavity structure ) at a depth of two meters by using Fortran programming code. But I have a problem how writing the cavity part in the code. The code below is for (Free-surface boundary conditions for elastic staggered-grid modeling schemes) when there is a cavity in the subsurface layer how to write code for the cavity part.
Cavity
• P-wave velocity (Vp ) = 0.0 m/s
• S-wave velocity (Vs) = 0.0 m/s
• Dnsity ( ρ) = 10-6 kg/m3
Cavity depth = 2m
Cavity length = 4m
Cavity width = 2.5m
cavity shape = rectangular
! v_x(x,z,t+dt/2) = v_x(x,z,t-dt/2)
! + dt*rho(x,z)^{-1}*[ dx(-) t_{xx}(x+dx/2,z,t) + dz(-) t_{xz}(x,z+dz/2,t) ] (1)
!
! v_z(x+dx/2,z+dz/2,t+dt/2) = v_z(x+dx/2,z+dz/2,t-dt/2)
! + dt*rho(x+dx/2,z+dz/2)^{-1}*[ dx(+) t_{xz}(x,z+dz/2,t) + dz(+) t_{zz}(x+dx/2,z,t) ] (2)
!
! t_{xx}(x+dx/2,z,t+dt) = t_{xx}(x+dx/2,z,t)
! + dt*lambda((x+dx/2,z)*[dx(+) v_x(x,z,t+dt/2) + dz(-) v_z(x+dx/2,z+dz/2,t+dt/2)] (3)
! + dt*2*mu((x+dx/2,z)*[dx(+) v_x(x,z,t+dt/2)]
!
! t_{zz}(x+dx/2,z,t+dt) = t_{zz}(x+dx/2,z,t)
! + dt*lambda((x+dx/2,z)*[dx(+) v_x(x,z,t+dt/2) + dz(-) v_z(x+dx/2,z+dz/2,t+dt/2)] (4)
! + dt*2*mu((x+dx/2,z)*[dz(-) v_z(x+dx/2,z+dz/2,t+dt/2)]
!
! t_{xz}(x,z+dz/2,t+dt) = t_{xz}(x,z+dz/2,t)
! + dt*mu((x,z+dz/2)*[dz(+) v_x(x,z,t+dt/2) + dx(-) v_z(x+dx/2,z+dz/2,t+dt/2)] (5)
!
!
!
!
!
program fdrayleighwave
implicit none
integer :: nx,nz,nt,ix,iz,lx,lz,i
integer :: itrec
integer :: k0x,klx,k0z,klz
integer :: ixs,izs
integer :: xdir,zdir,forwrd,bckwrd
parameter(xdir=1,zdir=3,forwrd=0,bckwrd=1) !Parameters used when calculating wave field components
real,parameter :: pi = 3.141592654
real,parameter :: twopi = 6.283185307
integer :: present,next
real :: vp1,vs1,dens1,vp2,vs2,dens2,width,height,source_x,source_z
real :: dx,dz,dt,dtrec !Space and time step，
real :: layer
real :: delay,a1,a2,tmax !delayDelay time, tmax maximum cycle time
real :: fp !Dominant frequency of source function
integer :: it !Time loop variable
integer :: nmedia,b_x,b_z
real :: snapshot_time,offset,interval
! the derivative operator half length.
integer :: lop !Differential operator half length
parameter(lop=3) !lop=3Then the accuracy is 6th order
real alpx(lop)
data alpx /1.20282,-0.08276,0.00950/ !Staggered Grid Finite Difference Coefficient
real alpz(lop)
data alpz /1.20282,-0.08276,0.00950/ !Staggered Grid Finite Difference Coefficient
real, allocatable ,dimension(:,:,:) :: txx,tzz,txz,txx_x,txx_z,tzz_x,tzz_z,txz_x,txz_z !Define stress array, variable length
real, allocatable ,dimension(:,:,:) :: vx,vz,vx_x,vx_z,vz_x,vz_z !Define the particle velocity component array
real, allocatable ,dimension(:,:) :: w1,w2,dxx,dzz ! Buffer field value
real, allocatable ,dimension(:,:) :: lam,mu,twomu !Elastic parameters
real, allocatable ,dimension(:,:,:) :: rhoinv ! density
real, allocatable ,dimension(:) :: src !Array of the function value of the hypocenter with respect to time
real, allocatable ,dimension(:) :: vxrec,vzrec !Array of particle velocity component values received by the geophone
integer, allocatable ,dimension(:) :: ixrec,izrec !Array of geophone positions
real, allocatable ,dimension(:,:) :: vp,vs,dens !Density and P-wave Velocity
real :: vzt1,vzt2
real :: vx_snapshot,vz_snapshot
integer nrec
! parameter(nrec=2)
OPEN(111,file="twolayer.txt",STATUS='OLD')
READ(111,*)tmax,dtrec
READ(111,*)snapshot_time
READ(111,*)delay
READ(111,*)fp
READ(111,*)nmedia
READ(111,*)vp1,vs1,dens1
READ(111,*)vp2,vs2,dens2
READ(111,*)width,height
READ(111,*)source_x,source_z
READ(111,*)dt,dx,dz
READ(111,*)nrec,offset,interval
READ(111,*)b_x,b_z
! recording of particle velocities at nrec stations
allocate(vxrec(nrec),vzrec(nrec))
allocate(ixrec(nrec),izrec(nrec))
! grid size
nx=width/dx+1
nz=height/dz+1
allocate(vp(nx,nz),vs(nx,nz),dens(nx,nz),lam(nx,nz),twomu(nx,nz),mu(nx,nz),rhoinv(nx,nz,2))
! grids
do ix=1,nx
do iz=1,12/dz
vp(ix,iz)=vp1
vs(ix,iz)=vs1
dens(ix,iz)=dens1
lam(ix,iz)=dens(ix,iz)*(vp(ix,iz)**2-2.0*vs(ix,iz)**2)*dt
mu(ix,iz)=dens(ix,iz)*(vs(ix,iz)**2)*dt
twomu(ix,iz)=2.0*mu(ix,iz)
rhoinv(ix,iz,1)=dt/dens(ix,iz)
enddo
do iz=12/dz+1,nz
vp(ix,iz)=vp2
vs(ix,iz)=vs2
dens(ix,iz)=dens2
lam(ix,iz)=dens(ix,iz)*(vp(ix,iz)**2-2.0*vs(ix,iz)**2)*dt
mu(ix,iz)=dens(ix,iz)*(vs(ix,iz)**2)*dt
twomu(ix,iz)=2.0*mu(ix,iz)
rhoinv(ix,iz,1)=dt/dens(ix,iz)
enddo
enddo
! p-wave velocity, s-wave velocity and density
! vp=1000.0
! vs=577.0
! dens=2000
! output files
! open(unit=11,file="vx_r1.txt")
! open(unit=12,file="vx_r2.txt")
! open(unit=13,file="vz_r1.txt")
! open(unit=14,file="vz_r2.txt")
open(unit=11,file="vx_kgs.txt")
open(unit=12,file="vz_kgs.txt")
open(unit=15,file="vx_snapshot.txt")
open(unit=16,file="vz_snapshot.txt")
! source position
! ixs=201
! izs=1
ixs=source_x/dx+1
izs=source_z/dz+1
! receiver position 1
do i=1,nrec
ixrec(i)=ixs+offset+(i-1)*interval
izrec(i)=1
enddo
! receiver position 2
! ixrec(2)=321
! izrec(2)=1
! time stepping and sampling interval for recording
! dtrec=0.001 !sampling interval
! dt=0.0002 !Time Step
itrec=nint(dtrec/dt)
! max time (Maximum cycle time)
! tmax=0.8
nt=nint(tmax/dt)+1 !Cycle time steps nint is the rounding function
! source time function
allocate(src(nt)) !Variable length
! The function of the seismic source with respect to time is a first-order Gaussian differential operator function
! parameterization is the same as for analytical calculation.
! delay=0.05
! fp=25.0 !Main frequency
a1=-twopi*fp*exp(0.5)
a2=-0.5*(twopi*fp)**2
call gauss(src,nt,dt,delay,a1,a2) !Call Gaussian function
! derivative operator half lengths
lx=lop
lz=lop
! the dynamic (time dependent) fields have buffer areas for "ix <1", "ix > nx", "iz < 1" and "iz > nz".
! the fields in these buffer areas are set to 0 initially and never updated. all
! dynamic fields are 0 in the buffers areas at all times. the width of the buffer areas are
! the same as the operator half length.
! Calculation area processing, the array starts from the negative
k0x=-(lx-1) !k0x=-2
klx=nx+lx !klx=404
k0z=-(lz-1) !k0z=-2
klz=nz+lz !klz=204
! dynamical fields , wave field value
allocate(txx(k0x:klx,k0z:klz,2),tzz(k0x:klx,k0z:klz,2),txz(k0x:klx,k0z:klz,2))
allocate(vx(k0x:klx,k0z:klz,2),vz(k0x:klx,k0z:klz,2),dxx(1:nx,1:nz),dzz(1:nx,1:nz))
allocate(vx_x(k0x:klx,k0z:klz,2),vx_z(k0x:klx,k0z:klz,2),vz_x(k0x:klx,k0z:klz,2),vz_z(k0x:klx,k0z:klz,2))
allocate(txx_x(k0x:klx,k0z:klz,2),tzz_x(k0x:klx,k0z:klz,2),txz_x(k0x:klx,k0z:klz,2))
allocate(txx_z(k0x:klx,k0z:klz,2),tzz_z(k0x:klx,k0z:klz,2),txz_z(k0x:klx,k0z:klz,2))
! work arrays for dynamical fields
allocate(w1(k0x:klx,k0z:klz),w2(k0x:klx,k0z:klz))
! the model parameters are not defined with buffers since this is not required. the vacuum
! is implemented by the dynamical fields. zero density and zero
! elastic parameters in the buffer areas and outside the buffer areas
! implicitly lead to zero dynamical fields in the buffers and outside the buffers.
! the free surface is implemented at depth "iz=1".(Free interface position)
! two inverse density, "rhoinv", grids are used, one for each component of particle
! velocity.
! the "lam" and "twomu" grids are used for the normal components of stress.
! the "mu" grid is used for the shear stress component.
! initial conditions
txx=0.0
txx_x=0.0
txx_z=0.0
tzz=0.0
tzz_x=0.0
tzz_z=0.0
txz=0.0
txz_x=0.0
txz_z=0.0
vx=0.0
vz=0.0
vx_x=0.0
vx_z=0.0
vz_x=0.0
vz_z=0.0
w1=0.0
w2=0.0
vxrec=0.0
vzrec=0.0
! spatial step lengths
! dx=1.0
! dz=1.0
layer=b_x*dx
! make derivative operators dependent of step length
alpx=alpx/dx
alpz=alpz/dz
! free surface: (Free interface parameter setting)
lam(1:nx,1)=0.0
twomu(1:nx,1)=0.5*twomu(1:nx,1)
rhoinv(1:nx,1,1)=2.0*rhoinv(1:nx,1,1)
! rhoinv(1:nx,1,1) is top nodes for inverse density grid for vx.
! the inverse density grid for vx has last index equal to 1.
! the inverse density grid for vz has last index equal to 2.
! time integration of elastic fields
do it=1,nt !Time loop
present=mod(it+1,2)+1 !Take the remainder
next=mod(it,2)+1 !Take the remainder
do iz=1,nz
do ix=1,nx
dxx(ix,iz)=0.0
dzz(ix,iz)=0.0
if(ix.le.(1+layer/dx))then
dxx(ix,iz)=(2*vp(ix,iz)/layer)*log(1/0.0001)*((1+layer/dx-ix)*dx/layer)**4
endif
if(ix.ge.(nx-layer/dx))then
dxx(ix,iz)=(2*vp(ix,iz)/layer)*log(1/0.0001)*((ix-(nx-layer/dx))*dx/layer)**4
endif
if(iz.ge.(nz-layer/dz))then
dzz(ix,iz)=(2*vp(ix,iz)/layer)*log(1/0.0001)*((iz-(nz-layer/dz))*dz/layer)**4
endif
enddo
enddo
! calculate vx (equation 1 above)
call der(xdir,bckwrd,txx(k0x,k0z,present),w1,nx,nz,alpx,lx,k0x,klx,k0z,klz)
vx_x(1:nx,1:nz,next)=((1-0.5*dt*dxx(1:nx,1:nz))/(1+0.5*dt*dxx(1:nx,1:nz)))*vx_x(1:nx,1:nz,present)+(1/(1+0.5*dt*dxx(1:nx,1:nz)))*rhoinv(1:nx,1:nz,1)/dx*w1(1:nx,1:nz)
call der(zdir,bckwrd,txz(k0x,k0z,present),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
vx_z(1:nx,1:nz,next)=((1-0.5*dt*dzz(1:nx,1:nz))/(1+0.5*dt*dzz(1:nx,1:nz)))*vx_z(1:nx,1:nz,present)+(1/(1+0.5*dt*dzz(1:nx,1:nz)))*rhoinv(1:nx,1:nz,1)/dx*w2(1:nx,1:nz)
! vx(1:nx,1:nz,next)=vx(1:nx,1:nz,present) + rhoinv(1:nx,1:nz,1)*(w1(1:nx,1:nz)+w2(1:nx,1:nz))
vx(1:nx,1:nz,present)=vx_x(1:nx,1:nz,present) + vx_z(1:nx,1:nz,present)
vx(1:nx,1:nz,next)=vx_x(1:nx,1:nz,next) + vx_z(1:nx,1:nz,next)
! calculate vz (equation 2 above)
call der(xdir,forwrd,txz(k0x,k0z,present),w1,nx,nz,alpx,lx,k0x,klx,k0z,klz)
vz_x(1:nx,1:nz,next)=((1-0.5*dt*dxx(1:nx,1:nz))/(1+0.5*dt*dxx(1:nx,1:nz)))*vz_x(1:nx,1:nz,present)+(1/(1+0.5*dt*dxx(1:nx,1:nz)))*rhoinv(1:nx,1:nz,1)/dx*w1(1:nx,1:nz)
call der(zdir,forwrd,tzz(k0x,k0z,present),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
! add directional force in z-direction to w2
call addfz(w2,src(it),ixs,izs,nx,nz,dx,dz,k0x,klx,k0z,klz)
vz_z(1:nx,1:nz,next)=((1-0.5*dt*dzz(1:nx,1:nz))/(1+0.5*dt*dzz(1:nx,1:nz)))*vz_z(1:nx,1:nz,present)+(1/(1+0.5*dt*dzz(1:nx,1:nz)))*rhoinv(1:nx,1:nz,1)/dx*w2(1:nx,1:nz)
vz(1:nx,1:nz,present)=vz_x(1:nx,1:nz,present) + vz_z(1:nx,1:nz,present)
vz(1:nx,1:nz,next)=vz_x(1:nx,1:nz,next) + vz_z(1:nx,1:nz,next)
! vz(1:nx,1:nz,next)=vz(1:nx,1:nz,present) + rhoinv(1:nx,1:nz,2)*(w1(1:nx,1:nz)+w2(1:nx,1:nz))
! recording: (Record the particle velocity component, call the recording sub-function, and output the recorded value)
call recordata(vx,vz,vxrec,vzrec,ixrec,izrec,nrec,nx,nz,k0x,klx,k0z,klz)
if(mod(it,itrec).eq.1) then
do i=1,nrec
write(11,*) vxrec(i)
write(12,*) vxrec(i)
enddo
! write(13,*) vzrec(1)
! write(14,*) vzrec(2)
endif
! change by yxz
if(it.eq.nint(snapshot_time/dt))then
do ix=1,nx
do iz=1,nz
vx_snapshot=0.5*(vx(ix,iz,1)+vx(ix,iz,2))
if(iz.eq.1)then
vzt1=0.25*(vz(ix-1,iz,1) + vz(ix,iz,1))
vzt2=0.25*(vz(ix-1,iz,2) + vz(ix,iz,2))
vz_snapshot=0.5*(vzt1+vzt2)
else
vzt1=0.25*(vz(ix-1,iz-1,1) + vz(ix,iz-1,1) &
+ vz(ix-1,iz,1) + vz(ix,iz,1))
vzt2=0.25*(vz(ix-1,iz-1,2) + vz(ix,iz-1,2) &
+ vz(ix-1,iz,2) + vz(ix,iz,2))
vz_snapshot=0.5*(vzt1+vzt2)
endif
write(15,*) ix,iz,vx_snapshot
write(16,*) ix,iz,vz_snapshot
enddo
enddo
endif
! divergence of particle velocity (Calculate the particle velocity component)
call der(xdir,forwrd,vx(k0x,k0z,next),w1,nx,nz,alpx,lx,k0x,klx,k0z,klz)
txx_x(1:nx,1:nz,next)=((1-0.5*dt*dxx(1:nx,1:nz))/(1+0.5*dt*dxx(1:nx,1:nz)))*txx_x(1:nx,1:nz,present)+(1/(1+0.5*dt*dxx(1:nx,1:nz)))*(lam(1:nx,1:nz)+twomu(1:nx,1:nz))/dx*w1(1:nx,1:nz)
call der(zdir,bckwrd,vz(k0x,k0z,next),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
txx_z(1:nx,1:nz,next)=((1-0.5*dt*dzz(1:nx,1:nz))/(1+0.5*dt*dzz(1:nx,1:nz)))*txx_z(1:nx,1:nz,present)+(1/(1+0.5*dt*dzz(1:nx,1:nz)))*lam(1:nx,1:nz)/dz*w2(1:nx,1:nz)
! time integration of txx (equation 3 above)
! txx(1:nx,1:nz,next)= txx(1:nx,1:nz,present) &
! + twomu(1:nx,1:nz)*w1(1:nx,1:nz)&
! + lam(1:nx,1:nz)*(w1(1:nx,1:nz)+w2(1:nx,1:nz))
txx(1:nx,1:nz,present)= txx_x(1:nx,1:nz,present) + txx_z(1:nx,1:nz,present)
txx(1:nx,1:nz,next)= txx_x(1:nx,1:nz,next) + txx_z(1:nx,1:nz,next)
! time integration of tzz (equation 4 above)
tzz_x(1:nx,1:nz,next)=((1-0.5*dt*dxx(1:nx,1:nz))/(1+0.5*dt*dxx(1:nx,1:nz)))*tzz_x(1:nx,1:nz,present)+(1/(1+0.5*dt*dxx(1:nx,1:nz)))*lam(1:nx,1:nz)/dx*w1(1:nx,1:nz)
tzz_z(1:nx,1:nz,next)=((1-0.5*dt*dzz(1:nx,1:nz))/(1+0.5*dt*dzz(1:nx,1:nz)))*tzz_z(1:nx,1:nz,present)+(1/(1+0.5*dt*dzz(1:nx,1:nz)))*(lam(1:nx,1:nz)+twomu(1:nx,1:nz))/dz*w2(1:nx,1:nz)
! tzz(1:nx,1:nz,next)= tzz(1:nx,1:nz,present) &
! + twomu(1:nx,1:nz)*w2(1:nx,1:nz)&
! + lam(1:nx,1:nz)*(w1(1:nx,1:nz)+w2(1:nx,1:nz))
tzz(1:nx,1:nz,present)= tzz_x(1:nx,1:nz,present) + tzz_z(1:nx,1:nz,present)
tzz(1:nx,1:nz,next)= tzz_x(1:nx,1:nz,next) + tzz_z(1:nx,1:nz,next)
! free surface **************Free interface processing*************************
tzz(1:nx,1,next)=0.0
call der(xdir,bckwrd,vz(k0x,k0z,next),w1,nx,nz,alpx,lx,k0x,klx,k0z,klz)
txz_x(1:nx,1:nz,next)=((1-0.5*dt*dxx(1:nx,1:nz))/(1+0.5*dt*dxx(1:nx,1:nz)))*txz_x(1:nx,1:nz,present)+(1/(1+0.5*dt*dxx(1:nx,1:nz)))*mu(1:nx,1:nz)/dx*w1(1:nx,1:nz)
call der(zdir,forwrd,vx(k0x,k0z,next),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
txz_z(1:nx,1:nz,next)=((1-0.5*dt*dzz(1:nx,1:nz))/(1+0.5*dt*dzz(1:nx,1:nz)))*txz_z(1:nx,1:nz,present)+(1/(1+0.5*dt*dzz(1:nx,1:nz)))*mu(1:nx,1:nz)/dz*w2(1:nx,1:nz)
! time integration of txz (equation 5 above)
! txz(1:nx,1:nz,next)=txz(1:nx,1:nz,present) &
! + mu(1:nx,1:nz)*(w1(1:nx,1:nz)+w2(1:nx,1:nz))
txz(1:nx,1:nz,present)= txz_x(1:nx,1:nz,present) + txz_z(1:nx,1:nz,present)
txz(1:nx,1:nz,next)= txz_x(1:nx,1:nz,next) + txz_z(1:nx,1:nz,next)
enddo !End of time loop
end program fdrayleighwave !End of main function program
!******************************************************************
! Calculate the external process of wave field components
!******************************************************************
subroutine der(axis,stagger,f,df,nx,nz,alp,l,k0x,klx,k0z,klz)
implicit none
integer :: axis,stagger
integer :: k0x,klx,k0z,klz,nx,nz,l
real :: f(k0x:klx,k0z:klz),df(k0x:klx,k0z:klz),alp(l)
integer :: ix,iz,ip,im
if(axis.eq.1) then
do iz=1,nz
do ix=1,nx
ip=ix-stagger !stagger Corresponds to forward or backward
im=ix+1-stagger
!6-order precision staggered grid finite difference calculation
df(ix,iz) = (f(ip+1,iz)-f(im-1,iz))*alp(1) &
+ (f(ip+2,iz)-f(im-2,iz))*alp(2) &
+ (f(ip+3,iz)-f(im-3,iz))*alp(3)
enddo
enddo
elseif(axis.eq.3) then
do iz=1,nz
ip=iz-stagger
im=iz+1-stagger
do ix=1,nx
df(ix,iz) = (f(ix,ip+1)-f(ix,im-1))*alp(1) &
+ (f(ix,ip+2)-f(ix,im-2))*alp(2) &
+ (f(ix,ip+3)-f(ix,im-3))*alp(3)
enddo
enddo
endif
return
end subroutine der
!**********************************************************************
!The external function of the wavefield value processing after adding the source term...
!that is, the processing of the wavefield value in the area near the source
!**********************************************************************
subroutine addfz(f,src,ixs,izs,nx,nz,dx,dz,k0x,klx,k0z,klz)
implicit none
integer :: ixs,izs
integer :: k0x,klx,k0z,klz,nx,nz
real :: f(k0x:klx,k0z:klz),src,dx,dz
! add vertical force in reference grid position. assume field f is on same sub grid as vz.
! use second order approximation to delta distribution to achieve this
if(ixs<=1.or.ixs>=nx) then
print*,'source cannot be on edge of (or outside) grid in x direction!'
stop
endif
if(izs>=nz) then
print*,'source cannot be on bottom edge of (or outside) grid in z direction!'
stop
endif
if(izs<1) then
print*,'source cannot be above top edge of grid in z direction!'
stop
endif
if(izs==1) then ! the source is set on the free-surface
f(ixs-1,izs)=f(ixs-1,izs)+0.25*src/(dx*dz)
f(ixs,izs)=f(ixs,izs)+0.25*src/(dx*dz)
else
f(ixs-1,izs-1)=f(ixs-1,izs-1)+0.25*src/(dx*dz)
f(ixs,izs-1)=f(ixs,izs-1)+0.25*src/(dx*dz)
f(ixs-1,izs)=f(ixs-1,izs)+0.25*src/(dx*dz)
f(ixs,izs)=f(ixs,izs)+0.25*src/(dx*dz)
endif
return
end subroutine addfz
!***************************************************************************
! Calculate the required output, the records at positions r1 and r2
! The return value is the particle velocity component at the receiving point vxrec(nrec),vzrec(nrec)
!***************************************************************************
subroutine recordata(vx,vz,vxrec,vzrec,ixrec,izrec,nrec,nx,nz,k0x,klx,k0z,klz)
implicit none
integer nrec !Number of detectors
integer :: ixrec(nrec),izrec(nrec) !Location of grid node where the geophone is located
integer :: k0x,klx,k0z,klz,nx,nz ! k0x,klx,k0z,klzIndicates the boundary position of the finite difference calculation area, nx, nz indicate the effective calculation area size
real :: vxrec(nrec),vzrec(nrec)
real :: vx(k0x:klx,k0z:klz,2),vz(k0x:klx,k0z:klz,2)
integer :: irec,ixr,izr ! Detector number cycle variable, geophone position variable
real :: vzt1,vzt2 !Used to find the comprehensive average of time and space when vz
! record vx and vz at position on reference grid and at "present" time step.
! thus, both fields require temporal interpolation. vx is on reference grid.
! vz require spatial interpolation also. interpolation operators are second
! order only, in both time and space.
do irec=1,nrec
ixr=ixrec(irec)
izr=izrec(irec)
!The following conditional statements are that the detector is outside the effective calculation area
if(ixr.le.1.or.ixr.ge.nx) then
print*,'receiver cannot be on edge of (or outside) grid in x direction!'
stop
endif
if(izr.ge.nz) then
print*,'receiver cannot be on bottom edge of (or outside) grid in z direction!'
stop
endif
if(izr.lt.1) then
print*,'receiver cannot be above top edge of grid in z direction!'
stop
endif
! vx is on reference grid. need only time interpolation.
! For vx, it only needs to average the value of the time before and after
vxrec(irec)=0.5*(vx(ixr,izr,1) + vx(ixr,izr,2))
! vz need spatial and temporal interpolation.
! For vz, space and time integration is required
if(izr.eq.1) then !The processing of the average wave field value at the free interface
vzt1=0.25*(vz(ixr-1,izr,1) + vz(ixr,izr,1)) ! 1 means the previous moment
vzt2=0.25*(vz(ixr-1,izr,2) + vz(ixr,izr,2)) ! 2 means the next moment
vzrec(irec)=0.5*(vzt1+vzt2)
else !Average processing of wave field values at other spatial positions
vzt1=0.25*(vz(ixr-1,izr-1,1) + vz(ixr,izr-1,1) &
+ vz(ixr-1,izr,1) + vz(ixr,izr,1))
vzt2=0.25*(vz(ixr-1,izr-1,2) + vz(ixr,izr-1,2) &
+ vz(ixr-1,izr,2) + vz(ixr,izr,2))
vzrec(irec)=0.5*(vzt1+vzt2)
endif
enddo
return
end subroutine recordata
!***********************************************************************
!Source function, the return value is a calculated array src(nt)
!***********************************************************************
subroutine gauss(src,nt,dt,t0,a1,a2)
implicit none
integer :: nt
real :: src(nt),dt,t0,a1,a2 !Source function value, time step, wavelet delay time, coefficients used in the calculation of Gaussian time function
integer :: it !Loop variable
real :: t,tau
do it=1,nt
t=(it-1)*dt
tau=t-t0
src(it)=2.0*a1*tau*exp(a2*tau*tau)
enddo
return
end subroutine gauss