!!!!!!!!!!!!!輸出快照的程序!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!輸出快照的程序!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!參數(shù)設(shè)置
!	nx-水平方向網(wǎng)格數(shù),nz-深度方向網(wǎng)格數(shù),dx-網(wǎng)格大小
!	nmax-時(shí)間采樣點(diǎn)數(shù),dt-采樣間隔,ixs-激發(fā)點(diǎn)水平方向網(wǎng)格號(hào),izs-激發(fā)點(diǎn)深度方向網(wǎng)格號(hào)
!	nsnap-快照輸出最大炮數(shù),ifsnap-快照輸出時(shí)間起點(diǎn)步長(zhǎng),idsnap-快照輸出時(shí)間增加步長(zhǎng)
!	isotype-震源類型:1/點(diǎn)震源,2/剪切震源,3/垂直震源,4/水平震源
!	iabmax-吸收邊界網(wǎng)格數(shù),topbc-地面邊界類型:free/自由邊界,abbc/吸收邊界
  	  parameter(nx=600,nz=600,dx=0.2,nz0=5)
      parameter(nmax=2001,dt=0.1e-3,ixs=300,izs=5)
      parameter(nsnap=40,ifsnap=100,idsnap=100,isotype=1,iabmax=80)
      parameter(topbc='free')
!!!!!!!!!!!!定義數(shù)組
      dimension vx(0:nx,0:nz),vz(0:nx,0:nz)
      dimension txx(0:nx,0:nz),tzz(0:nx,0:nz)
      dimension txz(0:nx,0:nz),bx(0:nx,0:nz),bz(0:nx,0:nz)
      dimension rlamu(0:nx,0:nz),rla(0:nx,0:nz),rmuxz(0:nx,0:nz)
      dimension rho(0:nx,0:nz),source(nmax),eponge(0:nx)
      dimension recx(0:nx,nmax),recz(0:nx,nmax),recxz(0:nx,nmax)
	  character*12 ffo11,ffo12,ffo21,ffo22,ffo23
	
	  data ffo21,ffo22,ffo23/'recx.grd','recz.grd','recxz.grd'/
!!!!!!!!!!!!!g1,g2時(shí)間2階-空間4階常數(shù),soufac震源因子
	  g1=-1.0/24.0
	  g2=9.0/8.0
	  dtdx=dt/dx
	  soufac=dt/(dx*dx)
	  isnap=0
!!!!!!!!!!!!!生成雷克子波
	  call wavelet(nmax,dt,source)
	  print*,'wavelet'
	  amax=100.0*amax/(0.25*soufac)
!!!!!!!!!!!!!生成地質(zhì)模型

	  call model(rlamu,rla,rmuxz,rho,bx,bz,dtdx,nx,nz,bx0,bz0,bz1,rmuxz1,rla0,rla1,rlamu0,rlamu1)
	  print*,'model'
!!!!!!!!!!!!!生成吸收邊界系數(shù)
	  xx=0.305/float(iabmax)
	  do ib=0,iabmax
	     aaa=(xx*float(iabmax-ib))**2
		 eponge(ib)=exp(-aaa)
      enddo
!!!!!!!!!!!!!數(shù)組初始化
	  do i=0,nx
	     do k=0,nz
		    vx(i,k)=0.0
			vz(i,k)=0.0
			txx(i,k)=0.0
			tzz(i,k)=0.0
			txz(i,k)=0.0
         enddo
	  enddo
!!!!!!!!!!!!!獲得質(zhì)點(diǎn)速度(4階)	  
	  do n=1,nmax
	    do k=2,nz-2
		  do i=2,nx-2
		    vx(i,k)=vx(i,k)+bx(i,k)*(g1*(txx(i+1,k)-txx(i-2,k))+g2*(txx(i,k)-txx(i-1,k))+g1*(txz(i,k+2)-txz(i,k-1))+g2*(txz(i,k+1)-txz(i,k)))
			vz(i,k)=vz(i,k)+bz(i,k)*(g1*(txz(i+2,k)-txz(i-1,k))+g2*(txz(i+1,k)-txz(i,k))+g1*(tzz(i,k+1)-tzz(i,k-2))+g2*(tzz(i,k)-tzz(i,k-1))) 
		  enddo
		enddo
!!!!!!!!!!!!單炮記錄賦值	
      do ix=0,nx
		recx(ix,n)=vx(ix,nz0)*1000000000
		recz(ix,n)=vz(ix,nz0)*1000000000
		recxz(ix,n)=(vx(ix,nz0)+vz(ix,nz0))*10000000/2.0
 	  enddo
!!!!!!!!!!!!!零碎東西
	    do k=1,1
		   do i=2,nx-2
		    vx(i,k)=vx(i,k)+bx(i,k)*(g1*(txx(i+1,k)-txx(i-2,k))+g2*(txx(i,k)-txx(i-1,k))+g1*(txz(i,k+2)-txz(i,k-1))+g2*(txz(i,k+1)-txz(i,k)))   
			enddo
         enddo

		do k=2,nz-2
		   do i=1,1
		      vz(i,k)=vz(i,k)+bz(i,k)*(g1*(txz(i+2,k)-txz(i-1,k))+g2*(txz(i+1,k)-txz(i,k))+g1*(tzz(i,k+1)-tzz(i,k-2))+g2*(tzz(i,k)-tzz(i,k-1))) 
		  enddo
		enddo

      if(mod(n,100).eq.0)print*,'timestep=',n
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!速度邊界條件

!!!!!!!!!!!!!左右邊界吸收條件
	    do ka=0,nz
		   do ia=0,iabmax
		      vx(ia,ka)=vx(ia,ka)*eponge(ia)
			  vz(ia,ka)=vz(ia,ka)*eponge(ia)
			  vx(nx-ia,ka)=vx(nx-ia,ka)*eponge(ia)
			  vz(nx-1-ia,ka)=vz(nx-1-ia,ka)*eponge(ia)
           enddo
		enddo
!!!!!!!!!!!!!上下邊界吸收條件
        do ka=0,iabmax
		   do ia=0,nx
		      if(topbc.eq.'abbc')then
			     vx(ia,ka)=vx(ia,ka)*eponge(ka)
				 vz(ia,ka)=vz(ia,ka)*eponge(ka)
			  endif
			  vx(ia,nz-1-ka)=vx(ia,nz-1-ka)*eponge(ka)
			  vz(ia,nz-ka)=vz(ia,nz-ka)*eponge(ka)
			enddo
		enddo
!!!!!!!!!!!!!自由表面邊界條件
		if(topbc.eq.'free')then
		  do i=2,nx-2
             vx(i,0)=vx(i,0)+bx0*(g1*(txx(i+1,0)-txx(i-2,0))+g2*(txx(i,0)-txx(i-1,0))+g1*(txz(i,2)+txz(i,1))+g2*(txz(i,1)-txz(i,0)))
		  enddo
		  do i=1,nx-2
		     vz(i,0)=vz(i,0)+2*bz0*(g1*tzz(i,1)+g2*tzz(i,0))
			 vz(i,1)=vz(i,1)+bz1*(g1*(tzz(i,2)+tzz(i,0))+g2*(tzz(i,1)-tzz(i,0))+g1*(txz(i+2,1)-txz(i-1,1))+g2*(txz(i+1,1)-txz(i,1)))
		  enddo
		endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!均勻彈性介質(zhì)中應(yīng)力-應(yīng)變之間的關(guān)系
        do k=2,nz-2
		   do i=2,nx-2
		      exx=g1*(vx(i+2,k)-vx(i-1,k))+g2*(vx(i+1,k)-vx(i,k))
			  ezz=g1*(vz(i,k+2)-vz(i,k-1))+g2*(vz(i,k+1)-vz(i,k))

			  if(k.eq.(nz-2))then
			    exz=vx(i,k)-vx(i,k-1)+g1*(vz(i+1,k)-vz(i-2,k))+g2*(vz(i,k)-vz(i-1,k))
				if(i.eq.1)exz=vx(1,k)-vx(1,k-1)+vz(1,k)-vz(0,k)
		      endif
			  if(k.ge.2.and.k.le.(nz-3))then
			    exz=g1*(vx(i,k+1)-vx(i,k-2))+g2*(vx(i,k)-vx(i,k-1))+g1*(vz(i+1,k)-vz(i-2,k))+g2*(vz(i,k)-vz(i-1,k))
				if(i.eq.1)then
				  exz=vz(1,k)-vz(0,k)+g1*(vx(1,k+1)-vx(i,k-2))+g2*(vz(1,k)-vz(i,k-1))
			    endif
			 endif

		   txx(i,k)=txx(i,k)+rlamu(i,k)*exx+rla(i,k)*ezz
		   tzz(i,k)=tzz(i,k)+rla(i,k)*exx+rlamu(i,k)*ezz
		   txz(i,k)=txz(i,k)+rmuxz(i,k)*exz
		enddo
	  enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!應(yīng)力應(yīng)變邊界條件
!!!!!!!!!!!!!自由邊界條件

	  if(topbc.eq.'free')then
	    do i=2,nx-2
		  txz(i,0)=0.0
		  txz(i,1)=txz(i,1)+g1*(vz(i+1,1)-vz(i-2,1))+g2*(vz(i,1)-vz(i-1,1))+(g1*(vx(i,2)+vz(i-1,0)-vz(i,0)-vx(i,0))+g2*(vx(i,1)-vx(i,0)))*rmuxz1
		enddo
		do i=1,nx-2
		   tzz(i,0)=tzz(i,0)+rlamu0*(vz(i,1)-vz(i,0))+rla0*(g1*(vx(i+2,0)-vx(i-1,0))+g2*(vx(i+1,0)-vx(i,0)))
		   tzz(i,1)=tzz(i,1)+rlamu1*(g1*(vz(i,3)-vz(i,0))+g2*(vz(i,2)-vz(i,1)))+rla1*(g1*(vx(i+2,1)-vx(i-1,1))+g2*(vx(i+1,1)-vx(i,1)))
		enddo
		do i=1,nx-2
		   txx(i,0)=txx(i,0)+rla0*(vz(i,1)-vz(i,0))+rlamu0*(g1*(vx(i+2,0)-vx(i-1,0))+g2*(vx(i+1,0)-vx(i,0)))
		   txx(i,1)=txx(i,1)+rla1*(g1*(vz(i,3)-vz(i,0))+g2*(vz(i,2)-vz(i,1)))+rlamu1*(g1*(vx(i+2,1)-vx(i-1,1))+g2*(vx(i+1,1)-vx(i,1)))
		enddo
	endif
!!!!!!!!!!!!!左,右吸收邊界條件
	 do ka=0,nz
	   do ia=0,iabmax
	      txz(ia,ka)=txz(ia,ka)*eponge(ia)
		  txx(ia,ka)=txx(ia,ka)*eponge(ia)
		  tzz(ia,ka)=tzz(ia,ka)*eponge(ia)
		  txz(nx-ia,ka)=txz(nx-ia,ka)*eponge(ia)
		  txx(nx-1-ia,ka)=txx(nx-1-ia,ka)*eponge(ia)
		  tzz(nx-1-ia,ka)=tzz(nx-1-ia,ka)*eponge(ia)
		enddo
	enddo
!!!!!!!!!!!!!上,下吸收邊界條件
	do ka=0,iabmax
	   do ia=0,nx
	      if(topbc.eq.'abbc')then
		  txz(ia,ka)=txz(ia,ka)*eponge(ka)
		  txx(ia,ka)=txx(ia,ka)*eponge(ka)
		  tzz(ia,ka)=tzz(ia,ka)*eponge(ka)
		  endif
		  txz(ia,nz-ka)=txz(ia,nz-ka)*eponge(ka)
		  txx(ia,nz-1-ka)=txx(ia,nz-1-ka)*eponge(ka)
		  tzz(ia,nz-1-ka)=tzz(ia,nz-1-ka)*eponge(ka)
		enddo
	enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!震源函數(shù)

!!!!!!!!!點(diǎn)震源
 if(isotype.eq.1)then
	   addsou=soufac*source(n)*0.25
	   txx(ixs-1,izs)=txx(ixs-1,izs)+addsou
	   txx(ixs,izs)=txx(ixs,izs)+addsou
	   txx(ixs-1,izs+1)=txx(ixs-1,izs+1)+addsou
	   txx(ixs,izs+1)=txx(ixs,izs+1)+addsou

	   tzz(ixs-1,izs)=tzz(ixs-1,izs)+addsou
	   tzz(ixs,izs)=tzz(ixs,izs)+addsou
	   tzz(ixs-1,izs+1)=tzz(ixs-1,izs+1)+addsou
	   tzz(ixs,izs+1)=tzz(ixs,izs+1)+addsou
	 endif
!!!!!!!!!!!!!剪切震源

	 if(isotype.eq.2)then
	   addsou=soufac*source(n)*0.25
	  vx(ixs-1,izs)=vx(ixs-1,izs)+addsou
	   vx(ixs,izs)=vx(ixs,izs)+addsou
	   vx(ixs-1,izs+1)=vx(ixs-1,izs+1)+addsou
	   vx(ixs,izs+1)=vx(ixs,izs+1)+addsou

	   vz(ixs-1,izs)=vz(ixs-1,izs)-addsou
	   vz(ixs,izs)=vz(ixs,izs)-addsou
	   vz(ixs-1,izs+1)=vz(ixs-1,izs+1)-addsou
	   vz(ixs,izs+1)=vz(ixs,izs+1)-addsou
	 endif
	 !!!!!!!!!!!!!垂直震源
       if(isotype==3)then
		  addsou=0.5*soufac*source(n)/rho(ixs-1,izs+1)
		  vz(ixs-1,izs+1)=vz(ixs-1,izs+1)+addsou
		  vz(ixs,izs+1)=vz(ixs,izs+1)+addsou
		endif
	 !!!!!!!!!!!!!水平震源
    	if(isotype.eq.4)then
		  addsou=0.5*soufac*source(n)/rho(ixs,izs)
		  vx(ixs,izs)=vx(ixs,izs)+addsou
		  vx(ixs,izs+1)=vx(ixs,izs+1)+addsou
		endif
!!!!!!!!!!!!!保存快照數(shù)據(jù)
!     if(mod(n,idsnap).eq.0.and.isnap.lt.nsnap)then
!        ffo11=' '
!		ffo12=' '
!		ffo22=' '
!		ffo11='kz'//char(isnap/10+48)//char(mod(isnap,10)+48)//'x.grd'	   
! 		ffo12='kz'//char(isnap/10+48)//char(mod(isnap,10)+48)//'z.grd'	
! 		ffo22='kz'//char(isnap/10+48)//char(mod(isnap,10)+48)//'x+z.grd'	 	     
! 		open(910,file=ffo11)
! 		write(910,'(a4)') 'DSAA'
! 		write(910,'(2i5)') int((nx)/2+1),int((nz)/2+1)
! 		write(910,'(2i5)') 0,int(nx*dx)
! 		write(910,'(2i5)') -int(nz*dx),0
! 		write(910,'(2i5)') -1000,1000
! 
! 		open(920,file=ffo12)
! 		write(920,'(a4)') 'DSAA'
! 		write(920,'(2i5)') int((nx)/2+1),int((nz)/2+1)
! 		write(920,'(2i5)') 0,int(nx*dx)
! 		write(920,'(2i5)') -int(nz*dx),0
! 		write(920,'(2i5)') -1000,1000
! 
! 		open(930,file=ffo22)
! 		write(930,'(a4)') 'DSAA'
! 		write(930,'(2i5)') int((nx)/2+1),int((nz)/2+1)
! 		write(930,'(2i5)') 0,int(nx*dx)
! 		write(930,'(2i5)') -int(nz*dx),0
! 		write(930,'(2i5)') -1000,1000
! 
! 	    do iz=nz,0,-2
! 		write(910,'(10g)')(vx(ix,iz)*1000000000,ix=0,nx,2)
! 		write(920,'(10g)')(vz(ix,iz)*1000000000,ix=0,nx,2) 
! 		write(930,'(10g)')((vz(ix,iz)+vx(ix,iz))*1000000000,ix=0,nx,2) 
! 		enddo
! 		close(910)
! 		close(920)
! 		close(930)
! 		isnap=isnap+1
! 	  endif
   enddo
!!!!!!!!!!!!保存單炮記錄數(shù)據(jù)
    open(930,file='x.grd')
	write(930,'(a4)') 'DSAA'
	write(930,'(2i6)') nx+1,nmax
	write(930,'(2i6)') 0,int(nx*dx)
	write(930,'(2i6)') -int(1000*(nmax-1)*dt),0
	write(930,'(2f5.1)') -0.1,0.1

	open(940,file='z.grd')
	write(940,'(a4)') 'DSAA'
	write(940,'(2i6)') nx+1,nmax
	write(940,'(2i6)') 0,int(nx*dx)
	write(940,'(2i6)') -int(1000*(nmax-1)*dt),0
	write(940,'(2f5.1)') -0.1,0.1

	open(950,file='x-z.grd')
	write(950,'(a4)') 'DSAA'
	write(950,'(2i6)') nx+1,nmax
	write(950,'(2i6)') 0,int(nx*dx)
	write(950,'(2i6)') -int(1000*(nmax-1)*dt),0
	write(950,'(2f5.1)') -0.1,0.1

	do n=nmax,1,-1
	  write(930,'(601g14.3)') (recx(ix,n),ix=0,nx)
	  write(940,'(601g14.3)') (recz(ix,n),ix=0,nx)
	  write(950,'(601g14.3)') (recxz(ix,n),ix=0,nx)
	enddo
	
	close(930)
	close(940)
	close(950)
!!輸出自己需要的部分記錄！！！！！！！！！！！！！！！！！！！！！！
	open(960,file='x.dat')
	open(970,file='z.dat')
	open(980,file='x-z.dat')
    do  n=1,2001
	  write(960,1000) (recx(ix,n),IX=int(nx/2)+50,int(nx/2)+150,2)
	  write(970,1000) (recz(ix,n),IX=int(nx/2),int(nx/2)+300)
	  write(980,1000) (recxz(ix,n),IX=int(nx/2)+50,int(nx/2)+150,2)
	ENDDO
1000 FORMAT(1X,301(F9.1,X))
	close(960)
	close(970)
	close(980)
   stop
    end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!模型參數(shù)設(shè)置
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     subroutine model(rlamu,rmu,rmuxz,rho,bx,bz,dtdx,nx,nz,bx0,bz0,bz1,rmuxz1,rmu0,rmu1,rlamu0,rlamu1)

     parameter(nmat=3,eps=1.e-6)
	 dimension rlamu(0:nx,0:nz),rmu(0:nx,0:nz)
     dimension rho(0:nx,0:nz),rmuxz(0:nx,0:nz)
     dimension bx(0:nx,0:nz),bz(0:nx,0:nz)
     dimension alpha(nmat),beta(nmat),dens(nmat)
 data alpha(1),beta(1),dens(1)/1000.0,500.0,1800.0/
	 data alpha(2),beta(2),dens(2)/2200.0,1100.0,2650.0/
	 data alpha(3),beta(3),dens(3)/1500.0,0.0,1000.0/
	 Real Ra,Rx,Ry
     	  Rx=300
    	  Ry=260
	 write(*,*) nx,nz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!第一層     
      do ix=0,nx
	  do iz=0,nz
		rlamu(ix,iz)=alpha(1)*alpha(1)*dens(1)
		rmu(ix,iz)=beta(1)*beta(1)*dens(1)
		rho(ix,iz)=dens(1)
	  enddo
	enddo
!	 do ix=0,nx
!	  do iz=101,nz
!		rlamu(ix,iz)=alpha(2)*alpha(2)*dens(2)
!		rmu(ix,iz)=beta(2)*beta(2)*dens(2)
!		rho(ix,iz)=dens(2)
!	  enddo
!	enddo

!!!!!!!!!!!!!結(jié)束模型構(gòu)建
    do i=0,nx
	   do k=0,nz
	      im=max0(0,i-1)
		  km=max0(0,k-1)
		  ii=min0(i,nx-1)
		  kk=min0(k,nz-1)
!!!!!!!!bx
          bx(i,k)=0.5*dtdx*((1.0/rho(im,kk))+(1.0/rho(ii,kk)))
!!!!!!!!bz
		  bz(i,k)=0.5*dtdx*((1.0/rho(ii,km))+(1.0/rho(ii,kk)))
!!!!!!!!muxz
         akan=1.0/(rmu(ii,kk)+eps)
		  akan=akan+(1.0/(rmu(im,kk)+eps))
		  akan=akan+(1.0/(rmu(ii,km)+eps))
		  akan=akan+(1.0/(rmu(im,km)+eps))
		  rmuxz(i,k)=(dtdx*4.0)/akan
		  if(rmuxz(i,k).lt.1.0)rmuxz(i,k)=0.0
		enddo
	enddo

    do i=0,nx-1
	   do k=0,nz-1
!!!!!!!!lambda
	      rmu(i,k)=dtdx*(rlamu(i,k)-(2.0*rmu(i,k)))
!!!!!!!!lambda+2mu		 
		  rlamu(i,k)=dtdx*rlamu(i,k)
	   enddo
	enddo
!!!!!!!!自由表面fd系數(shù)

    bx0=bx(0,0)
	bz0=bz(0,0)
	bz1=bz(1,0)
	rlamu0=rlamu(0,0)
	rlamu1=rlamu(1,0)
	rmu0=rmu(0,0)
	rmu1=rmu(1,0)
	rmuxz0=rmuxz(0,0)
!!!!!!!!!!!增加一些條件用來減少“IF”代碼
!!!!!!!!!!!vx

    do i=0,nx
	   bx(i,0)=0.0
	   bx(i,nz-1)=0.0
	   bx(i,nz)=0.0
	enddo
	do k=0,nz
	   bx(0,k)=0.0
	   bx(1,k)=0.0
	   bx(nx-1,k)=0.0
	   bx(nx,k)=0.0
	enddo
!!!!!!!!!!!vz
	do i=0,nx
	   bz(i,0)=0.0
	   bz(i,nz-1)=0.0
	   bz(i,nz)=0.0
	enddo
	do k=0,nz
	   bz(0,k)=0.0
	   bz(1,k)=0.0
	   bz(nx-1,k)=0.0
	   bz(nx,k)=0.0
	enddo
!!!!!!!!!!!txx 和 tzz
    do i=0,nx
	   rlamu(i,0)=0.0
	   rlamu(i,1)=0.0
	   rlamu(i,nz-2)=0.0
	   rlamu(i,nz-1)=0.0
	   rlamu(i,nz)=0.0
	   rmu(i,0)=0.0
	   rmu(i,1)=0.0
	   rmu(i,nz-2)=0.0
	   rmu(i,nz-1)=0.0
	   rmu(i,nz)=0.0
	enddo
	do k=0,nz
	   rlamu(0,k)=0.0
	   rlamu(nx-1,k)=0.0
	   rlamu(nx,k)=0.0
	   rmu(0,k)=0.0
	   rmu(nx-1,k)=0.0
	   rmu(nx,k)=0.0
	enddo
!!!!!!!!!!!txz
    do i=0,nx
	   rmuxz(i,0)=0.0
	   rmuxz(i,1)=0.0
	   rmuxz(i,nz-1)=0.0
	   rmuxz(i,nz)=0.0
	enddo
	do k=0,nz
	   rmuxz(0,k)=0.0
	   rmuxz(nx,k)=0.0
	enddo

    return
	end
!!!!!!!!!!!!!sfreq-子波主頻,tl-周期數(shù)

    subroutine wavelet(nt,dt,wav)
	dimension wav(*)
	do i=1,nt
	wav(i)=0.0
	enddo

    sfreq=80
	lwav=1+int(1.55/(sfreq*dt))

    do i=1,lwav
	   wav(i)=ricker(float(i-1)*dt,sfreq)
	enddo

    return
	end


    function ricker(time,sfreq)
	  pix=4.05367
	  g0=0.5
	  g1=0.565
	  g2=0.105
	  t0=time*sfreq*pix
	  ricker=0.0
	  ricker=g0*cos(t0)-g1*cos(2*t0)+g2*cos(3*t0)
	  return
	end