9-模拟理想海上勘探条件下合成野外磁带-p100.txt

地球物理大师Claebout的几个有用的程序

	implicit none

	!# synthetic marine data tape movie generation

	integer:: it,nt,ih,nh,is,ns,iz,nz,it0,iy,kbyte
	real p(512),b(512),refl(25,16),z(25),geol(25),rand

	open(3,file="plot.txt") 
	nt=512;nh=48;ns=10;nz=25;kbyte=1
	
	do iz=1,nz		!# reflector depth (共(nz+1)层,z(iz)表示第iz层的深度)
		call random_seed(iz)
		call random_number(rand)
		z(iz)=nt*rand	!# random() is on the interval (0.,1.)
		write(*,*) rand,z(iz)
	enddo
	
	do iz=1,nz		!# reflector strength with depth
		call random_seed(iz)
		call random_number(rand)
		geol(iz)=2*rand-1.
	enddo
	
	do is=1,ns		!# give texture to the geology
		do iz=1,nz
			call random_seed(iz)
			call random_number(rand)
			refl(iz,is)=(1.+rand)*geol(iz)
		enddo
	enddo
	
	do it=1,nt		!# prepare a wavelet
		b(it)=exp(-it*0.08)*sin(0.5*it-0.5)
	enddo

	do is=ns,1,-1		!# shots run backwards
		do ih=1,nh	!# down cable h=(g-s)/2 
			iy=(is-1)+(ih-1)	!# y=(g+s)/2=midpoint
			iy=1+(iy-ns*(iy/ns))	!# periodic with midpoint
			do it=1,nt
				p(it)=0.
			enddo
			do iz=1,nz		!# add in a hyperbola for each layer
				it0=sqrt(z(iz)**2+100.*(ih-1)**2)
				do it=1,nt-it0	!# add in the wavelet
					p(it+it0)=p(it+it0)+refl(iz,iy)*b(it)
				enddo
			enddo
			write(3,rec=kbyte)(p(it),it=1,nt) 
			kbyte=kbyte+nt*4
		enddo
	enddo
	close(3)
	stop
	!end