seismic avo.txt

这是运用Fortran语言编写的地震波AVO正演的程序源码

program main
implicit none
real,dimension(1:200)::ra1,ta1,angleh1,angle1,y,x1,x2,r1,rangle1,r2,r,b
real,dimension(1:200)::ra2,ta2,angleh2,angle2
real,dimension(1:200)::wav
real,dimension(200,500):: seis,ref
real,dimension(500)::velp,vels,den,buff
real vp1,vp2,vp3,vs1,vs2,vs3,lo1,lo2,lo3,vpb1,vpb2,h1,h2
real vph1,vph2,vsh1,vsh2,loh1,loh2,vpc1,vpc2,vsc1,vsc2,loc1,loc2
real  lja1,lja2,x,f,f1,l,w0,scaler,c,sss
real f0,t0,samp,wl,samptt,cc,ff,t,lb2,bsb1,bsb2,bsb3
real m1,n1,w1,q1,z1,j1,s1,s2,s3,s4,m2,n2,w2,q2,z2,j2,cs1,cs2,res,rrm
integer i,n,j,k,d,s,dd,lsyn,lw,lb1
real ,parameter ::pi=3.1415926,a=12.5
n=0
x=0.01
!参数的值
lb1=125
vp1=2000.0
vs1=sqrt((0.5-bsb1)/(1.0-bsb1))*vp1
bsb1=0.3
lo1=1.9

vp2=2500.0
vs2=sqrt((0.5-bsb2)/(1.0-bsb2))*vp1
bsb2=0.4
lo2=2.2

vp3=3000.0
vs3=sqrt((0.5-bsb3)/(1.0-bsb3))*vp1
bsb3=0.3
lo3=2.4

h1=160.0
h2=210.0
vpb1=vp2/vp1
vpb2=vp3/vp2

vph1=(vp1+vp2)/2.0
vsh1=(vs1+vs2)/2.0
loh1=(lo1+lo2)/2.0
vpc1=vp2-vp1
vsc1=vs2-vs1
loc1=lo2-lo1

vph2=(vp3+vp2)/2.0
vsh2=(vs3+vs2)/2.0
loh2=(lo3+lo2)/2.0
vpc2=vp3-vp2
vsc2=vs3-vs2
loc1=lo3-lo2




!......................................第一层的angle1

lja1=asin(1/vpb1)

lja2=asin(1/vpb2)

print*,"critical angle is："
print*,lja1
print*,"angle1 of the first floor:"
      
   do i=1,50
       l=a*(i-1)
       ra1(i)=atan(l/(2.0*h1))
rangle1(i)=ra1(i)*180/pi
if (ra1(i).lt.lja1)  then
   ta1(i)=asin(vpb1*sin(ra1(i)))  
else
  ta1(i)=0.0
	end if

angleh1(i)=(ra1(i)+ta1(i))/2.0
angle1(i)=angleh1(i)*180/pi
write(*,*)"i=",i,ra1(i),angle1(i)
enddo

!.............................第2层的angle2

do i=1,50
  do
  y(i)=a*(i-1)
  x1(i)=x
  F=vp2**2*(x-y(i))**2*(x**2+h2**2)-vp1**2*x**2*((x-y(i))**2+h1**2)
  F1=2*vp2**2*(x-y(i))*(x**2+h2**2)+2*vp2**2*x*(x-y(i))**2-2*vp1**2*x*((x-y(i))**2+h1**2)-2*vp1**2*x**2*(x-y(i))
  x=x1(i)-F/F1
  x2(i)=x
  n=n+1
  
  if(abs(x2(i)-x1(i))<1e-5)exit
  end do
   !print*,"changdu"
   !print*,"i=",i,x2(i)
  end do


print*,"critical angle is："
print*,lja2
print*,"angle2 of the second floor"
do i=1,50
   ra2(i)=atan(x2(i)/(2.0*h2))
if (ra2(i).lt.lja2)  then
    ta2(i)=asin(vpb2*sin(ra2(i)))   
	else
  ta2(i)=90.0
	end if
	angleh2(i)=(ra2(i)+ta2(i))/2.0
    angle2(i)=angleh2(i)*180/pi
	
!write(*,*)"i=",i,ra2(i),angle2(i)
enddo



!.......................calculate  reflectance 反射系数的求取
m1=vpc1/vph1
n1=loc1/loh1
w1=vsc1/vsh1
q1=(vsh1/vph1)**2
z1=(m1+n1)/2.0
j1=m1/2.0-4*q1*w1-2*q1*n1    !negative
write(*,*) "j1=",j1

open(unit=500,file='r1floor.dat',status='old')
!first  floor
do i=1,50
s1=(sin(angleh1(i)))**2
s2=(tan(angleh1(i)))**2
r1(i)= z1+j1*s1+m1*(s2-s1)/2.0

write(*,*) angle1(i),sin(angle1(i))**2,r1(i)
!write(500,*) angle1(i),r1(i)
enddo


!...............................second floor
print*,"second floor"
m2=vpc2/vph2
n2=loc2/loh2
w2=vsc2/vsh2
q2=(vsh2/vph2)**2
z2=(m2+n2)/2.0
j2=m2/2.0-4*q2*w2-2*q2*n2    !negative
!open(unit=600,file='r2floor.dat',status='old')
   do i=1,50
    s3=(sin(angleh2(i)))**2
    s4=(tan(angleh2(i)))**2
    r2(i)= z2+j2*s2+m2*(s4-s3)/2.0

!write(600,*) angle2(i),r2(i)
!write(*,*) angle2(i),r2(i)
enddo



            print*,'子波的值为：'
!..........................Calculate seismic wavelet 计算地震子波
     !open(unit=100,file='zb4.dat',status='old')
                f0 = 35.0
	            t0 = 0.0
	            samp  = 2.0
	            wl = 126
				w0 = 0.0
	            samptt= samp/1000.0    
		        lw=wl/samp +0.5
			    c=1.0
			    sss=5.5                 	         
			    if( mod(lw*1.,2.).eq.0 ) then
                  lw = lw+1 
		   		  t0 = t0+samp*2.0   
                endif


print*,'lw====',lw

       call wavelet(f0,t0,w0,c,sss,samptt,wav,lw)
	   lsyn = lb1+lw-1
	 
	  print*,'地震记录为'

	   scaler=100.0
	   open(unit=100,file='dzjl001.dat',status='old')
       open(unit=200,file='dzjl002.dat',status='old')
	   open(unit=300,file='dzjl003.dat',status='old')
       do i=1,50
          do j=1,lb1
	   	     if (j.eq.41) then
		        cs1=angleh1(i)
	            r(j)=cs1
		     elseif (j.eq.76) then
 		        cs2=angleh2(i)
		        r(j)=cs2
	         elseif(j.le.1.or.j.ge.127)then
                r(j)=0.0
             else
                r(j)=0.0                  
             endif
	         res=0.0
	         if(r(j).lt.res) r(j)=0.0
		  enddo
		  call fold(lb1,r,lw,wav,lsyn,buff,scaler)
		  do j=1,lsyn
		!  if(i.eq.19) print*,'i===',i,j,buff(j),r(j)
		     seis(i,j)=buff(j)
	      enddo 
	   enddo
	   !print*,"第",i,"道时：",' 各点的反射系数为：'                 
!	  write(*,*) r
      
			
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
!					do   d=1,lb1		
!					 do  s=1,lw				  
!					  k=d+s-1
!					  c(k)=r(d)*b(s)+c(k)					 
!                      enddo					  
!					  enddo
!					   write(*,*) c

!------Write Seismic data
             do j=1,lsyn
                write(100,301) j*samp,(seis(i,j),i=1,18)
			   write(200,301) j*samp,(seis(i,j),i=19,36)
			   write(300,301) j*samp,(seis(i,j),i=37,50)

			 enddo
301          format(1x,f6.1,18f10.3)

      stop
	  end
              !     do  i=1,k
	              
                !      enddo
              !  format(1x,i5,f15.7)
	                 ! close(20)

			   
   !synthet seismic

    !  do   400   d=1,lb1
    !         rrm=0.0
    !  do   500   s=1,lb2
    !         rrm = r(d+s-1)*b(lb1-s+1)+rrm
    !   500 continue
	!	        k= lb1
    !           c(k)=rrm
	!		    write(100,*)   d,c(k)
	                  
	                
   !    400 continue
                       
	    
              
!				
!	     open(20,file='dzjl.execle',status='new')
!
 !        do 200 i=1,sum
!	       write(20,301) i,b(i)
!200        continue
!301        format(1x,i5,f15.7)
!	       close(20)
      Subroutine   fold(la,a,lb,b,lc,c,scaler)    
      dimension    a(la),b(lb),c(lc)                                
      call   zero(lc,c)
	  do   10   i=1,la
	  do   10   j=1,lb
	  k    = i+j-1
10	  c(k) = a(i)*b(j)+c(k)
      do   20   i=1,la+lb-1
20    c(i) = scaler*c(i)                                           
!          K    = MIN0( LA, LB )/2                          
!          LEN  = MAX0( LA, LB )                      
!          DO   30   I = 1, LEN                                           
!          J    = K + I - 1                                               
!30        C(I) = C(J)                                                    
      return                                                          
      end

	Subroutine   zero( l0, o )                                      
      dimension    o(l0)                                              
      do   10   i = 1, l0                                             
10    o(i)  = 0.                                                      
      return                                                          
      end

	Subroutine wavelet(f,t0,w0,c,s,sampt,b,lb)
	dimension  b(lb)
	w   = w0/180.*3.1415927
	pi2 = 2.*3.1415927*f
	a   = -f*f*alog(s)
	t   = -t0
	print*,'lb===',lb
      do  10  i = 1, lb
        b(i)= c * exp( a*t*t ) * cos( pi2*t + w )
        ccc=b(i)
!       print*,'i =',i,' b= ',ccc
10    t   = t + sampt
      return
      end