

	program elastic

	parameter(nx=100,ny=nx,nz=ny)
	parameter(maxnt=500)
	parameter(maxnr=60)
	parameter(nop=4)

	real, dimension(nx,ny,nz) :: ux,uy,uz,sxx,syy,szz,sxy,sxz,syz,mxz
	real, dimension(nx,ny,nz) :: lam,mu,rho,vol,rx,ry,rz

	common /cinc / dx

	real s(maxnt),si(maxnt), sis(maxnt), M0
	real, dimension(maxnt*maxnr) ::seisx,seisy,seisz
	real, dimension(maxnt*maxnr) ::rotx,roty,rotz
	
	character*50 seisfile

	write(*,*) ' ' 
	write(*,*)' 3-D Elastic FD  '
	write(*,*)' ' 


	seisfile='Data/rot'
	dt=.005        ! time increment
	dx=100.		! space increment
	vp0=5196.
	vs0=3000.
	rho0=2600.
	f0=10.
	aa=1/f0
	nt=300		! number of time steps
	icheck=1
	isnap=25
	iout=1
	isx=nx/2
	isy=ny/2
	isz=int(1./2*nz)	
	irec=2
	ir0=10
	ir1=90
	izrec=int(1./4.*nz)
	M0=2.e16/(dx**3)
	mxz=0.
	rise_time=.1
	nrise_time=rise_time/dt


! initialise source
	
	rmax=0.
	rimax=0.
	rmin=1.e10
	rimin=1.e10

	t0=3.*aa/dt
	sum=0;
	do i=1,nt
	s(i)=exp(- 1./aa**2*((i-t0)*dt)**2)
	sum=sum+s(i)
	enddo
	s=M0*s/sum
	
!	s=0.
!	s(1)=1.
	open(10,file='source')
	do i=1,nt
	write(10,*)s(i)
	enddo
	close(10)
		

! initialise elastic constants

	rho=rho0
	mu=rho0*vs0**2
	lam=rho0*vp0**2-2*mu

	write(*,*)' Opening  snapshot files ... '
	open(10,file=seisfile(1:lnblnk(seisfile))//'_xs')
	open(11,file=seisfile(1:lnblnk(seisfile))//'_ys')
	open(12,file=seisfile(1:lnblnk(seisfile))//'_zs')
	open(30,file=seisfile(1:lnblnk(seisfile))//'_rxs')
	open(31,file=seisfile(1:lnblnk(seisfile))//'_rys')
	open(32,file=seisfile(1:lnblnk(seisfile))//'_rzs')

! checking

	write(*,*)' Checking ... '
	write(*,*)' seisfile = ',seisfile
	write(*,*)' lnblnk(seisfile) : ',lnblnk(seisfile)
	write(*,*)' nx = ',nx
	write(*,*)' ny = ',ny
	write(*,*)' nz = ',nz
	write(*,*)' nt = ',nt
	write(*,*)' dx = ',dx
	write(*,*)' isx = ',isx
	write(*,*)' isy = ',isy
	write(*,*)' isz = ',isz
	write(*,*)' vp0 = ',vp0
	write(*,*)' vs0 = ',vs0
	write(*,*)' rho0 = ',rho0
	write(*,*)' max vs = ',maxval(sqrt(mu/rho))
	write(*,*)' max vp = ',maxval(sqrt((lam+2*mu)/rho))
	write(*,*)' icheck = ',icheck 
	write(*,*)' isnap  = ',isnap 
	write(*,*)' M0  = ',M0
	write(*,*)' *********************** '
	write(*,*)' Stability : ',dt/dx*vp0
	write(*,*)' *********************** '	


	start=secnds(0.0)

	do it=1,nt	! time extrapolation

! Strike slip

	mxz(isx,isy,isz)=s(it)

! checking 

	if(mod(it,icheck)==0)then
	write(*,*)' Timestep : ',it
	write(*,*)minval(ux),' < ux < ',maxval(ux)
	write(*,*)minval(uy),' < uy < ',maxval(uy)
	write(*,*)minval(uz),' < uz < ',maxval(uz)
	time=secnds(start)
	total_time=time/icheck*nt
	write(*,*)' Pred. total : ',total_time/60.,'mins. (',
     &float(it)/nt*100,'% done)'
	start=secnds(0.0)
	endif

! acceleration

	ux = ux + (1/rho * (d('x',sxx,'-') + d('y',sxy,'-') + 
     &	d('z',sxz+mxz,'-') ))*dt
	uy = uy + (1/rho * (d('x',sxy,'+') + d('y',syy,'+') + 
     &	d('z',syz,'-') ))*dt
	uz = uz + (1/rho * (d('x',sxz+mxz,'+') + d('y',syz,'-') + 
     &	d('z',szz,'+') ))*dt

! stresses

	vol = (d('x',ux,'+') + d('y',uy,'-') + d('z',uz,'-'))

	sxx = sxx + 
     &		  ( lam *  vol + 2*mu * d('x',ux,'+')  ) * dt

	syy = syy + 
     &		  ( lam *  vol + 2*mu * d('y',uy,'-')  ) * dt

	szz = szz +
     &		  ( lam *  vol + 2*mu * d('z',uz,'-')  ) * dt 


	sxz = sxz + ( mu * (d('z',ux,'+') + d('x',uz,'-') ) ) * dt

	sxy = sxy + ( mu * (d('y',ux,'+') + d('x',uy,'-') ) ) * dt

	syz = syz + ( mu * (d('z',uy,'+') + d('y',uz,'+') ) ) * dt

! rotation 

	rx=d('y',uz,'+') - d('z',uy,'+')
	ry=d('z',ux,'+') - d('x',uz,'-')
	rz=d('x',uy,'-') - d('y',ux,'+')

! snapshot output 

	if(mod(it,isnap)==0)then
	write(*,*)' Outputting  snap ',it/isnap,' ... '
	do i=1,nx,iout
	do j=1,nz,iout
	write(10,*)ux(i,ny/2,j)
	write(11,*)uy(i,ny/2,j)
	write(12,*)uz(i,ny/2,j)
	write(30,*)rx(i,ny/2,j)
	write(31,*)ry(i,ny/2,j)
	write(32,*)rz(i,ny/2,j)
	enddo
	enddo
	endif

! seismogram output

	ii=0
	do i=ir0,ir1,irec
	ii=ii+1
	seisx((ii-1)*nt+it)=ux(i,ny/2,izrec)
	seisy((ii-1)*nt+it)=uy(i,ny/2,izrec)
	seisz((ii-1)*nt+it)=uz(i,ny/2,izrec)
	rotx((ii-1)*nt+it)=rx(i,ny/2,izrec)
	roty((ii-1)*nt+it)=ry(i,ny/2,izrec)
	rotz((ii-1)*nt+it)=rz(i,ny/2,izrec)
	enddo
	nr=ii

	enddo

! sesimogram output

	
	open(unit=20,file=seisfile(1:lnblnk(seisfile))//'_x')
	open(unit=21,file=seisfile(1:lnblnk(seisfile))//'_y')
	open(unit=22,file=seisfile(1:lnblnk(seisfile))//'_z')
	open(unit=40,file=seisfile(1:lnblnk(seisfile))//'_rx')
	open(unit=41,file=seisfile(1:lnblnk(seisfile))//'_ry')
	open(unit=42,file=seisfile(1:lnblnk(seisfile))//'_rz')

	write(*,*)' Outputting',nr,' traces '
	do i=1,nt*nr
	write(20,*)seisx(i)
	write(21,*)seisy(i)
	write(22,*)seisz(i)
	write(40,*)rotx(i)
	write(41,*)roty(i)
	write(42,*)rotz(i)
	enddo
	

	contains

	function d(comp,f,dir)

	real f(:,:,:)
	real d(ubound(f,1),ubound(f,2),ubound(f,3))
	character*1 comp,dir
	common /cinc /dx
      parameter(maxnop=8)
	real g(maxnop)

! derivative operators
!
	if(nop==4) then
	  g(1)=.0416666666
	  g(2)=-1.125
	endif

! the symmetry
!
	do i=nop/2+1,nop
	  g(i) = -g(nop-i+1)
	enddo

	d=0.

	if(comp=='x') then

	if ( dir == '-' ) then

	   if(nop==2) then
		  do i=2,nx
			do j=1,ny
			  do k=1,nz
				d(i,j,k)=(f(i,j,k)-f(i-1,j,k))/dx
			  enddo
			enddo
		  enddo	
	    endif		

	    if(nop==4) then
	  		do i=3,nx-1
			  do j=1,ny
				do k=1,nz
				d(i,j,k)=( g(1)*f(i-2,j,k)+ g(2)*f(i-1,j,k)+ 
     &				 g(3)*f(i  ,j,k)+ g(4)*f(i+1,j,k) )/dx
				enddo
			  enddo
			enddo 
	     endif		

      elseif(dir=='+') then

      if(nop==2) then
		do i=1,nx-1
	      do j=1,ny
	        do k=1,nz
			d(i,j,k)=(f(i+1,j,k)-f(i,j,k))/dx	
		    enddo
		  enddo
		enddo
      endif

      if(nop==4) then
		do i=2,nx-2
	      do j=1,ny
	        do k=1,nz
			d(i,j,k)=(g(1)*f(i-1,j,k)+ g(2)*f(i,  j,k)+ 
     &           g(3)*f(i+1,j,k)+ g(4)*f(i+2,j,k))/dx 
			enddo
		  enddo
		enddo
      endif	
  		 	
	endif

	elseif( comp == 'y' ) then
  
      if ( dir == '-' ) then
   
      if(nop==2) then
	  	do i=1,nx
	      do j=2,ny
	        do k=1,nz
			d(i,j,k)=(f(i,j,k)-f(i,j-1,k))/dx
			enddo
		  enddo
   	    enddo
      endif		

      if(nop==4) then
	  	do i=1,nx
	      do j=3,ny-1
	        do k=1,nz
			d(i,j,k)=( g(1)*f(i,j-2,k)+ g(2)*f(i,j-1,k)+ 
     &            g(3)*f(i,j,  k)+ g(4)*f(i,j+1,k) )/dx 
			enddo
		  enddo
		enddo
      endif		

      elseif( dir == '+' ) then

      if(nop==2) then
	  	do i=1,nx
	      do j=1,ny-1
	        do k=1,nz
			d(i,j,k)=(f(i,j+1,k)-f(i,j,k))/dx	
		    enddo
		  enddo
		enddo
      endif
  
      if(nop==4) then
	  	do i=1,nx
	      do j=2,ny-2
	        do k=1,nz
			d(i,j,k)=( g(1)*f(i,j-1,k)+ g(2)*f(i,j,  k)+ 
     &              g(3)*f(i,j+1,k)+ g(4)*f(i,j+2,k) )/dx 
			enddo
		  enddo
		enddo
      endif

      endif

      elseif( comp == 'z') then

      if ( dir == '-' ) then

      if(nop==2) then
		do i=1,nx
	      do j=1,ny
	        do k=2,nz
			d(i,j,k)=(f(i,j,k)-f(i,j,k-1))/dx
	        enddo
	      enddo
	    enddo
      endif		

      if(nop==4) then
	    do i=1,nx
	      do j=1,ny
	        do k=3,nz-1
			d(i,j,k)=( g(1)*f(i,j,k-2)+ g(2)*f(i,j,k-1)+ 
     &              g(3)*f(i,j,k  )+ g(4)*f(i,j,k+1) )/dx
			enddo
		  enddo
		enddo
      endif		

      elseif(dir =='+') then

      if(nop==2)then
		do i=1,nx
	      do j=1,ny
	        do k=1,nz-1
			d(i,j,k)=(f(i,j,k+1)-f(i,j,k))/dx
		    enddo
		  enddo
		enddo	
      endif

      if(nop==4) then
		do i=1,nx
	      do j=1,ny
	        do k=2,nz-2
			d(i,j,k)=( g(1)*f(i,j,k-1)+ g(2)*f(i,j,k  )+ 
     &               g(3)*f(i,j,k+1)+ g(4)*f(i,j,k+2) )/dx 
			enddo
		  enddo
		enddo
      endif		
	 		
      endif

      endif


	end function d

	function lnblnk(c)
	
	integer lnblnk
	character*50 c

	ii=0
	do i=1,50
	if(c(i:i) ==' ')ii=ii+1
	enddo
	lnblnk=50-ii

	end function lnblnk


	end program elastic
