test5.f90

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

!  test5.f90 
!
!  FUNCTIONS:
!	test5      - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: test5
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

	program test5

	implicit none

	! Variables
	!# 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
		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,*)(p(it),it=1,nt) !kbyte,
!			kbyte=kbyte+nt*4
		enddo
	enddo
	close(3)
	stop

	end program test5