! Location of geophone 1
ixrec(1)=1500 ! x-direction
izrec(1)=20 ! z-direction
! Location of geophone 2
ixrec(2)=1500 x-direction
izrec(2)=300 ! z-direction
Hallo,
wouldn’t it be a solution to just increase the size of the arays
ixrec and izrec
Of course, you may also have to change other parts (loops ?) of the code to
account for the increased number of geophones.
program fdrayleighwave
implicit none
integer :: nx,nz,nt,ix,iz,lx,lz
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 :: dens,vp,vs !Density and P-wave Velocity
real :: dx,dz,dt,dtrec !Space and time step，
real :: delay,a1,a2,tmax !delay Delay time, tmax maximum cycle time
real :: fp ! frequency Dominant of source function
integer :: it !Time loop variable
! 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 !Define stress array, variable length
real, allocatable, dimension(:,:,:) :: vx,vz !Define the particle velocity component array
real, allocatable ,dimension(:,:) :: w1,w2 ! 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
integer nrec
parameter(nrec=2)
! recording of particle velocities at nrec stations
allocate(vxrec(nrec),vzrec(nrec))
allocate(ixrec(nrec),izrec(nrec))
! grid size Calculation area size
nx=2500
nz=2000
! p-wave velocity, s-wave velocity and density
vp=2000.0
vs=1200.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=17,file="vx_0.15s.txt")
open(unit=18,file="vz_0.15s.txt")
! source position (Source location)
ixs=500
izs=10
! Location of geophone 1
ixrec(1)=1500
izrec(1)=20
! Location of geophone 2
ixrec(2)=1500
izrec(2)=300
! 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))
! Dynamic field value calculation array
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.
allocate(lam(nx,nz),twomu(nx,nz),mu(nx,nz),rhoinv(nx,nz,2)) !Model parameter array
! initial conditions
txx=0.0
tzz=0.0
txz=0.0
vx=0.0
vz=0.0
w1=0.0
w2=0.0
vxrec=0.0
vzrec=0.0
! grids nodes
lam=dens*(vp**2-2.0*vs**2)*dt
mu=dens*(vs**2)*dt
twomu=2.0*mu
rhoinv=dt/dens
! spatial step lengths
dx=1.0
dz=1.0
! 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
! calculate vx (equation 1 above)
call der(xdir,bckwrd,txx(k0x,k0z,present),w1,nx,nz,alpx,lx,k0x,klx,k0z,klz)
call der(zdir,bckwrd,txz(k0x,k0z,present),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
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))
! calculate vz (equation 2 above)
call der(xdir,forwrd,txz(k0x,k0z,present),w1,nx,nz,alpx,lx,k0x,klx,k0z,klz)
call der(zdir,forwrd,tzz(k0x,k0z,present),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
! Add an external force in the z direction to w2
call addfz(w2,src(it),ixs,izs,nx,nz,dx,dz,k0x,klx,k0z,klz)
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
write(11,*) vxrec(1)
write(12,*) vxrec(2)
write(13,*) vzrec(1)
write(14,*) vzrec(2)
endif
if (it.eq.751)then
do iz=1,nz
do ix=1,nx
write(17,*) (vx(ix,iz,next))
write(18,*) (vz(ix,iz,next))
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)
call der(zdir,bckwrd,vz(k0x,k0z,next),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
! 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))
! time integration of tzz (equation 4 above)
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))
! 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)
call der(zdir,forwrd,vx(k0x,k0z,next),w2,nx,nz,alpz,lz,k0x,klx,k0z,klz)
! 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))
enddo !End of time loop
end program fdrayleighwave !End of main function program
Please help
Hallo,
while I am not really able to provide you with a rewritten program,
here are some ideas.
integer nrec
parameter(nrec=2)
! recording of particle velocities at nrec stations
allocate(vxrec(nrec),vzrec(nrec))
allocate(ixrec(nrec),izrec(nrec))
The number of stations is stored in the parameter nrec.
So if the number of stations is increased, you have to change nrec accordingly.
ok I will try , Thank you very much
Note that you very likely have do change other parts of your code.
For example the output files seem to written separately for the geophones.
So if the number of geophones is increased to another, but fixed number, you can just add
new lines for the.
If the number of geophones should be handleed flexibly you may need a loop construct.
Note that I have not really scanned through the entire code, you may also have to change other parts of the code.