hechengjilu.f90

此程序用Fortran语言编写

!=======================================================
!     功能：地表模型合成地震记录
!=======================================================
      program synth_surface
	  character *60 title,infile1,infile2,infile3,outfile
	  integer n,m,n1
	  real  dt,amp
	  integer*2 head(120)
	  real,allocatable::wave(:),r(:),y(:)
!======================================================
!     input parameter

      open(12,file="model_par.txt",status="old")
	  do i=1,3
	  read(12,'(a60)')Title
	  end do
!s1:  input wavelet file
	  read(12,'(a30,a30)') Title,infile1 
	  write(*,'(a30,a30)') Title,infile1
!s2:  input re_coefficient	   
      read(12,'(a30,a30)') Title,infile2 
	  write(*,'(a30,a30)') Title,infile2    
   
!s3:  name outfile	 
	  read(12,'(a30,a30)') Title,outfile
	  write(*,'(a30,a30)') Title,outfile 
!s4:  input wavelet length	     
	  read(12,'(a30,i8)') Title,n
	  write(*,'(a30,i8)') Title,n         
!s5:  input the record length	   
      read(12,'(a30,i8)') Title,m
	  write(*,'(a30,i8)') Title,m 
!s6:  input the amplitude	        
      read(12,'(a30,i8)') Title,amp
	  write(*,'(a30,i8)') Title,amp
!s7:  input sampling space(ms)	          
      read(12,'(a30,f8.1)') Title,dt 
	  write(*,'(a30,f8.1)') Title,dt 
!s8   input output trace Sum_number
      read(12,'(a30,i8)') Title,n_tra 
	  write(*,'(a30,i8)') Title,n_tra	       
	  close(12) 
	  allocate(wave(1:n),r(1:m),y(1:m))
!=============================================
!     inputwavelet
	  open(unit=15,file=infile1)
	  do i=1,n
	  read(15,*)wave(i)
	  end do
	  close(15)
!=============================================
!     calculate re_coefficient
      call cal_re_co(m,dt,infile2,r)


!=============================================
!     convolution
	  call convolution(n,m,m,wave,r,y,amp)
!=============================================
!     output data_sgy
      len=4*m+240
	  open(unit=20,file=outfile,access='direct',recl=len, status='replace')
	  head(58)=m
	  head(59)=int(dt*1000)
	  do j=1,n_tra  
      write(20,rec=j) head,(y(i),i=1,m)
	  end do
	  close(20)

	end program


  subroutine convolution(n1,n2,m,f,r,y,amp)
!======================================================
!     Parameters:
!       f     :seismic wavelet
!       r     :reflection conefficient
!       y     :result of convolution
!       n1    :length of the seismic wavelet
!       n2    :length of reflection coefficient series
!       m     :length of seismic record
!======================================================
	  integer  n1,n2,m,i,j,k,kk
	  real  f(n1),r(n2),y(m),amp
	   y=0.0
	   do j=1,n1
	      kk=n2-j+1
		  do i=1,kk
		     k=i+j-n1/2
			 if(k>0) y(k)=y(k)+amp*f(j)*r(i)
		  end do
	   end do
	  
   end subroutine
   
   subroutine cal_re_co(m,dt,infile,r)
	 character *(*)  infile
	 character * 60 title
     integer n,m
     real s1,s2,dt,r(1:m) 
     real,allocatable::v(:),den(:),h(:),Vr(:),re(:),t(:)
     integer,allocatable::num(:)
     open(unit=12,file=infile,status='old')
!======================================
     do k=1,3
	 read(12,*) title
	 end do
! input the number of layers
     read(12,'(a30,i8)') title,n
     write(*,'(a30,i8)') title,n
     allocate(v(1:n),den(1:n),h(1:n-1),Vr(1:n-1),re(1:n-1),t(1:n-1),num(1:n-1))
! input the velocity density and the depth
     do i=1,n-1
     read(12,*) v(i),den(i),h(i)
     end do
     read(12,*) v(n),den(n)
     close(12)
!===================================
! calculate the RMS velocity
     Vr(1)=v(1)
     do i=2,n-1
     s1=h(1)*v(1);s2=h(1)/v(1)
	   do j=2,i
         s1=s1+(h(j)-h(j-1))*v(j)
	     s2=s2+(h(j)-h(j-1))/v(j)   
       end do
	 Vr(i)=sqrt(s1/s2)
    end do
! calculate the re_coefficient
    r=0.0
    do i=1,n-1
    re(i)=(v(i+1)*den(i+1)-v(i)*den(i))/(v(i+1)*den(i+1)+v(i)*den(i))
    t(i)=2*h(i)/Vr(i)
    num(i)=int(t(i)*1e+3/dt)
    r(num(i))=re(i) 
    end do
    open(unit=17,file="re_co.txt",status='replace')
    do i=1,m
    write(17,*) r(i)  
    end do
    close(17)
  end subroutine