program modelcreater
implicit none
integer :: nx,nz
integer :: i,j,k,n_rec
integer :: isx,isz
real :: sx,sz,dx,dz,offset0,drec
real,allocatable,dimension(:) :: recx,recz
character, allocatable,dimension(:,:) :: icode
real,parameter :: pi=3.1415926
real :: kk,h
integer :: mm,ang
nx=1201 !grid number in x direction
nz=601 !grid number in z direction
sx=50.0 !horizontal position of the source (unit: m)
offset0=1.0 !the nearest offset in horizontal direction
drec=1.0 !the trace interval in horizontal direction
n_rec=100 !the total trace numbers
dx=1.0 !the horizontal spatial step length (unit: m)
dz=1.0 !the vertical spatial step length (unit: m)
isx=nint(sx/dx) !the horizontal grid position of the source
ang=10.0 !Dipping angle of free surface (just for dipping-free-surface model)
h=10.0 !Real thickness of surface layer
allocate(recx(1:n_rec),recz(1:n_rec))
allocate(icode(1:nx,1:nz))
!***********************************************************************************************
!---call the subrotine to configure the matrix of icode
!call Right_dipping_twolayer(icode,nx,nz,ang,h,pi,dx,dz) !for twolayer model with right dipping free surface and layer interface
!call Left_dipping_twolayer(icode,nx,nz,ang,h,pi,dx,dz) !for twolayer model with left dipping free surface and layer interface
!call Left_dipping_twolayer_fault(icode,nx,nz,ang,h,pi,dx,dz) !for twolayer fault model with left dipping free surface and layer interface
call sintopo_twolayer(icode,nx,nz,ang,h,pi,dx,dz) ! for twolayer model with sin function case topographical free surface
!***********************************************************************************************
!***********************************************************************************************
call threelayer2(icode,nx,nz,ang,h,pi,dx,dz)
! output model the model
open(1, file='threelayer2.txt',status='unknown')
do j=1,nz
write(1,100)(icode(i,j),i=1,nx)
enddo
100 format(<nx>a1)
close(1)
open(2,file='threelayer2_srf.txt',status='unknown')
do i=1,nx
do j=1,nz
write(2,200)(i-1)*dx,-1*(j-1)*dz,icode(i,j)
enddo
enddo
200 format(2e,1x,1a)
close(2)
end program modelcreater
!!*****************************************************************************************
!!*************************************topography twolayer model***************************
subroutine sintopo_twolayer(icode,nx,nz,ang,h,pi,dx,dz)
implicit none
integer :: nx,nz
integer :: i,j
character :: icode(1:nx,1:nz)
real :: pi,kk,h,dx,dz
integer :: ang,mm
do i=1,nx
mm=nint(80*sin(pi*(i+200)/400)+0.5)+100
do j=1,nz
if(j<mm)then
icode(i,j)='1'
elseif(j>mm.and.j<mm+nint(h/dz))then !twolayer halfspace with topographic free surface
icode(i,j)='3'
elseif(j>=mm+nint(h/dz))then !twolayer halfspacer with topographic layer interface
icode(i,j)='4'
else
icode(i,j)='2'
endif
enddo
enddo
return
end subroutine sintopo_twolayer
@1988, please explain the problem you encounter: bad result, compilation error? More people will help you if you clearly state the problem.
2 Likes