Write a code that shows the behavior of the seismic wave in the case of cavitation

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 = regtangular

@1988 is this homework? If so, you should mark it as such. Also you should state what kind of help you seek, i.e., what you tried, what worked what didn’t work etc.

1 Like

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')

! 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

do i=1,nrec
ixrec(i)=ixs+offset+(i-1)*interval
izrec(i)=1
enddo

!  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
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
!**********************************************************************
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
!***************************************************************************
! 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
``````

Not a lot of references in the Fortran code ?

I work in vibration propagation through structural objects, using Finite element analysis, so am not familiar with your x/z field model approach.
Anyway I will start off the comments !!
Basically, you need to set the boundary conditions in your mesh x/z where the cavity occurs to S=P=0, then progress through your time loop.
It could be as simple as that, although you also need to understand how the vibration field around a cavity is incompatible with your x/z mesh assumption.
Does it need to become an x/y/z mesh to model non uniformity in the y direction ?
Assuming x is longitudinal, z is vertical, does your mesh refinement need to change near the edges of the rectangle to model a non-uniform field around the edge of the cavity ?
Do your field iteration equations cope with non-uniform grid spacing ?

Again I do not work with a field modelling approach, but you have to understand how the vibration field is distorted around the cavity and how your x/z field assumptions may no longer be valid for this case.
How can you adapt, as x/y/z field modelling will come with a huge increase in memory storage and time/calculation cost.
Do the search

Great opportunity to learn about Fortran coarrays, preferably from those who are already using them.

1 Like