	
	include 'psv_deriv.f'

	subroutine psv_evolution

	include 'common.h'

	use derivative

	real, dimension (nx,nz)  :: rtemp1,rtemp2
	real, dimension (4,nx,nz) :: e,s 

!hpf$ align with tt :: rtemp1,rtemp2
!hpf$ align (*,:,:) with tt(:,:) :: e,s


	if(mod(it,icheck)==0)then
	write(*,*)' Time step   ',it,'   Max(u_r)  =  ',maxval(u(it2,2,:,:))
	write(*,*)' Time step   ',it,'   Min(u_r)  =  ',minval(u(it2,2,:,:))
	write(*,*)' Time step   ',it,'   Max(u_t)  =  ',maxval(u(it2,1,:,:))
	write(*,*)' Time step   ',it,'   Min(u_t)  =  ',minval(u(it2,1,:,:))
c	write(*,*)'                   Maxloc(u)  =  ',maxloc(u)
	if(maxval(u)>1.)then
	write(*,*)' unstable at  it = ',it
	stop
	endif
	endif


c symmetry with respect to source location

	forall(ix=1:isx-1)u(it2,1,isx-ix,:)=-u(it2,1,isx+ix,:)	
	forall(ix=1:isx)u(it2,2,isx+1-ix,:)=u(it2,2,isx+ix,:)	

c free surface


        forall (iz=1:izfree-1) u(it2,2,:,izfree-iz)=
     &                                          u(it2,2,:,izfree+iz)
        forall (iz=1:izfree) u(it2,1,:,izfree+1-iz)=
     &                                          u(it2,1,:,izfree+iz)


c eoshift simulations at sides and at bottom

	u(it2,1,nx-izfree:nx,:)=0.
	u(it2,2,nx-izfree:nx,:)=0.
	u(it2,1,:,nz-izfree:nz)=0.
	u(it2,2,:,nz-izfree:nz)=0.


c end-of-shift simulation to omit eoshift 


c absorbing boundaries


c calculate strain elements :
c		e1 = eps_r r
c		e2 = eps_phi phi
c		e3 = eps_theta theta
c		e4 = eps_r theta

c e1
c
c e1 = d_r u^r 

	e=0.
	rtemp2=u(it2,2,:,:)
	call pder(rtemp1,rtemp2,dr,nx,nz,-1,2,nop)	
	e(1,:,:)=rtemp1

c e2 = u^r/r + 1/r cotg theta u^theta

	call inter(rtemp1,rtemp2,
     &		  nx,nz,-1,2,nop)

	e(2,:,:)=rtemp1*r1inv

c we need this later

	e(3,:,:) = e(2,:,:)

	rtemp2=u(it2,1,:,:)
	call inter(rtemp1,rtemp2,nx,nz,-1,1,nop)

	e(2,:,:) = e(2,:,:)+ cot2*rtemp1*r1inv

c e3 = 1/r d_theta u^theta + 1/r u^r

        call pder(rtemp1,rtemp2,dphi,nx,nz,-1,1,nop)

	e(3,:,:) = e(3,:,:) + rtemp1*r1inv

c e4 = .5 ( 1/r d_theta u^r + d_r u^theta - 1/r u_theta)

	rtemp2=u(it2,2,:,:)	
        call pder(rtemp1,rtemp2,dphi,nx,nz,1,1,nop)
	s(1,:,:)=rtemp1	
	rtemp2=u(it2,1,:,:)
        call pder(rtemp1,rtemp2,dr,nx,nz,1,2,nop)

	e(4,:,:) = 0.5*(r2inv*s(1,:,:) + rtemp1)

        call inter(rtemp1,rtemp2,nx,nz,1,2,nop)

	e(4,:,:)=e(4,:,:)-0.5*rtemp1*r2inv

c set off diag elements zero at source

	e(4,isx,:)=0.

c calculate stress


	s=0.

	s(1,:,:) = (lam+2*mu)*e(1,:,:)+lam*(e(2,:,:)+e(3,:,:))

	s(2,:,:) = (lam+2*mu)*e(2,:,:)+lam*(e(1,:,:)+e(3,:,:))

	s(3,:,:) = (lam+2*mu)*e(3,:,:)+lam*(e(2,:,:)+e(1,:,:))

	s(4,:,:)=2*mu*e(4,:,:)

c moments

	if(source_type==2)then
	s(1,isx:isx+1,isz:isz+1)=s(1,isx:isx+1,isz:isz+1)+so(1,it)
	s(2,isx:isx+1,isz:isz+1)=s(2,isx:isx+1,isz:isz+1)+so(1,it)
	s(3,isx:isx+1,isz:isz+1)=s(3,isx:isx+1,isz:isz+1)+so(1,it)
	endif

c boundary conditions (stress) -> antisymmetry

	forall(iz=1:izfree)s(1,:,izfree+1-iz)=-s(1,:,izfree+iz)
	forall(iz=1:izfree)s(2,:,izfree+1-iz)=-s(2,:,izfree+iz)
	forall(iz=1:izfree)s(3,:,izfree+1-iz)=-s(3,:,izfree+iz)
	forall(iz=1:izfree-1)s(4,:,izfree-iz)=-s(4,:,izfree+iz)
	s(4,:,izfree)=0.	


c calculate acceleration

c e(1,:,:) and e(2,:,:) will carry acc theta and r, resp.
c e(3-4,:,:) act as temporay arrays

c for u_r :
c     e(3,:,:)=d_r sigma_r^r
c     e(4,:,:)=d_theta sigma_theta^r

	rtemp2=s(1,:,:)
        call pder(rtemp1,rtemp2,dr,nx,nz,1,2,nop)
	e(3,:,:)=rtemp1
	rtemp2=s(4,:,:)
        call pder(rtemp1,rtemp2,dphi,nx,nz,-1,1,nop)
	e(4,:,:)=rtemp1

	e(2,:,:) =   r2inv * e(4,:,:)+e(3,:,:) 

	
c     e(1,:,:) = s(1) -> u_r
c     e(3,:,:) = s(2) -> u_r
c     e(4,:,:) = s(3) -> u_r

	rtemp2=s(1,:,:)
	call inter(rtemp1,rtemp2,nx,nz,1,2,nop) 
	e(1,:,:)=rtemp1
	rtemp2=s(2,:,:)
	call inter(rtemp1,rtemp2,nx,nz,1,2,nop) 
	e(3,:,:)=rtemp1
	rtemp2=s(3,:,:)
	call inter(rtemp1,rtemp2,nx,nz,1,2,nop) 
	e(4,:,:)=rtemp1
	rtemp2=s(4,:,:)
	call inter(rtemp1,rtemp2,nx,nz,-1,1,nop) 

	e(2,:,:) = rhoinv * ( e(2,:,:) + r2inv*(2*e(1,:,:)-e(3,:,:)-e(4,:,:)
     &			      + rtemp1*cot2 )    ) 
	
c for u_theta:

c 	e(3,:,:)=d_r sigma_r^theta
c	e(4,:,:)=d_theta sigma_theta^theta
c	e(5,:,:)=s2 -> u_phi


	rtemp2=s(4,:,:)
        call pder(rtemp1,rtemp2,dr,nx,nz,-1,2,nop)
	e(3,:,:)=rtemp1
	rtemp2=s(3,:,:)
        call pder(rtemp1,rtemp2,dphi,nx,nz,1,1,nop)
	e(4,:,:)=rtemp1

	forall (ix=1:isx-1) e(1,ix,:) = r1inv(ix,:) * e(4,ix,:) + e(3,ix,:) 
	forall (ix=isx+1:nx) e(1,ix,:) = r1inv(ix,:) * e(4,ix,:) + e(3,ix,:) 

	rtemp2=s(4,:,:)
        call inter(rtemp1,rtemp2,nx,nz,-1,2,nop)
	e(3,:,:)=rtemp1
	rtemp2=s(2,:,:)
        call inter(rtemp1,rtemp2,nx,nz,1,1,nop)
	e(4,:,:)=rtemp1
	rtemp2=s(3,:,:)
        call inter(rtemp1,rtemp2,nx,nz,1,1,nop)

	forall (ix=1:isx-1)  e(1,ix,:) = e(1,ix,:) + r1inv(ix,:)* 
     &                       ( (-e(4,ix,:)+rtemp1(ix,:))*cot1(ix,:)+3*e(3,ix,:) )
	forall (ix=isx+1:nx) e(1,ix,:) = e(1,ix,:) + r1inv(ix,:)* 
     &                       ( (-e(4,ix,:)+rtemp1(ix,:))*cot1(ix,:)+3*e(3,ix,:) )
	e(1,isx,:) = 0.
	
	e(1,:,:)=e(1,:,:)*rhoinv


c add source term

	if(source_type==1)then
	e(2,:,:)=e(2,:,:)+gauss*so(1,it)*rhoinv
	endif	

c extrapolate


	u(it3,:,:,:)= 2*u(it2,:,:,:) - u(it1,:,:,:) + e(1:2,:,:)*(dt**2)



	return
	end	
