	
	subroutine psv_evolution


	include 'common.h'


	real, dimension (nx,nz)  :: rtemp1,rtemp2
	real, dimension (nx,nz) :: e1,e2,e3,e4
!hpf$	distribute (*,block), shadow(0,4) :: rtemp1,rtemp2
!hpf$	distribute (*,block), shadow(0,4) :: e1,e2,e3,e4



	if(mod(it,icheck)==0)then
	write(*,*)' Time step   ',it,':'
	write(*,*)minval(u21),' < v_t < ',maxval(u21)
	write(*,*)minval(u22),' < v_r < ',maxval(u22)
c	write(*,*)minval(window),' < window < ',maxval(window)
	if(maxval(u21)>100..or.maxval(u22)>100.)then
	write(*,*)' unstable at  it = ',it
	stop
	endif
	endif

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Preparation and calculation of accelerations cccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	s1=s1*window
	s2=s2*window
	s3=s3*window
	s4=s4*window

c moments

	if(source_type==2)then
	   if(alpha==0.)then
	   s1(isx:isx+1,isz:isz+1)=s1(isx:isx+1,isz:isz+1)+so(1,it)
	   s2(isx:isx+1,isz:isz+1)=s2(isx:isx+1,isz:isz+1)+so(1,it)
	   s3(isx:isx+1,isz:isz+1)=s3(isx:isx+1,isz:isz+1)+so(1,it)
	   else
	   s1=s1+gauss*so(1,it)
	   s2=s2+gauss*so(1,it)
	   s3=s3+gauss*so(1,it)
	   endif
	endif

c boundary conditions (stress) -> antisymmetry

	forall(iz=1:izfree)s1(:,izfree+1-iz)=-s1(:,izfree+iz)
	forall(iz=1:izfree)s2(:,izfree+1-iz)=-s2(:,izfree+iz)
	forall(iz=1:izfree)s3(:,izfree+1-iz)=-s3(:,izfree+iz)
	forall(iz=1:izfree-1)s4(:,izfree-iz)=-s4(:,izfree+iz)
	s4(:,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

        call pder(e3,s1,dr,nx,nz,1,2,nop)
        call pder(e4,s4,dphi,nx,nz,-1,1,nop)
	e2 =   r2inv * e4+e3 

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

	call inter(e1,s1,nx,nz,1,2,nop) 
	call inter(e3,s2,nx,nz,1,2,nop) 
	call inter(e4,s3,nx,nz,1,2,nop) 
	call inter(rtemp1,s4,nx,nz,-1,1,nop) 
	e2 = window * rhoinv * ( e2 + r2inv*(2*e1-e3-e4+rtemp1*cot2)) 
	
c for u_theta:

        call pder(e3,s4,dr,nx,nz,-1,2,nop)
        call pder(e4,s3,dphi,nx,nz,1,1,nop)

	forall (ix=1:isx-1)  e1(ix,:) = r1inv(ix,:) * e4(ix,:) + e3(ix,:) 
	forall (ix=isx+1:nx) e1(ix,:) = r1inv(ix,:) * e4(ix,:) + e3(ix,:) 

        call inter(e3,s4,nx,nz,-1,2,nop)
        call inter(e4,s2,nx,nz,1,1,nop)
        call inter(rtemp1,s3,nx,nz,1,1,nop)

	forall (ix=1:isx-1)  e1(ix,:) = e1(ix,:) + r1inv(ix,:)* 
     &       ( (-e4(ix,:)+rtemp1(ix,:))*cot1(ix,:)+3*e3(ix,:) )
	forall (ix=isx+1:nx)  e1(ix,:) = e1(ix,:) + r1inv(ix,:)* 
     &       ( (-e4(ix,:)+rtemp1(ix,:))*cot1(ix,:)+3*e3(ix,:) )

	e1(isx,:) = 0.
	
	e1 = e1*rhoinv*window


c add source term

	if(source_type==1)then
	e2=e2+gauss*so(1,it)*rhoinv
	endif	

ccccccccccccccccccccccccccccccccccccccccccccc
c extrapolate velocities cccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccc

	u21 = u21  + e1*dt
	u22 = u22  + e2*dt

ccccccccccccccccccccccccccccccccccccccccccccc
c end of velocity calculation ccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccc





cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccc preparation and calculation of stresses cccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c symmetry with respect to source location

	forall(ix=1:isx-1)u21(isx-ix,:)  =-u21(isx+ix,:)	
	forall(ix=1:isx)  u22(isx+1-ix,:)= u22(isx+ix,:)	

c free surface


        forall (iz=1:izfree-1) u22(:,izfree-iz)=
     &                                          u22(:,izfree+iz)
        forall (iz=1:izfree) u21(:,izfree+1-iz)=
     &                                          u21(:,izfree+iz)


c eoshift simulations at sides and at bottom

	u21(nx-izfree:nx,:)=0.
	u22(nx-izfree:nx,:)=0.
	u21(:,nz-izfree:nz)=0.
	u22(:,nz-izfree:nz)=0.

	u21=window*u21
	u22=window*u22


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 

	call pder(e1,u22,dr,nx,nz,-1,2,nop)	

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

	call inter(rtemp1,u22,nx,nz,-1,2,nop)
	e2=rtemp1*r1inv

c we need this later

	e3 = e2
	call inter(rtemp1,u21,nx,nz,-1,1,nop)
	e2 = e2+ cot2*rtemp1*r1inv

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

        call pder(rtemp1,u21,dphi,nx,nz,-1,1,nop)
	e3 = e3 + rtemp1*r1inv

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

        call pder(rtemp1,u22,dphi,nx,nz,1,1,nop)
	s1=rtemp1	
        call pder(rtemp1,u21,dr,nx,nz,1,2,nop)

	e4 = 0.5*(r2inv*s1 + rtemp1)
        call inter(rtemp1,u21,nx,nz,1,2,nop)
	e4 = e4-0.5*rtemp1*r2inv

c set off diag elements zero at source

	e4(isx,:)=0.

c calculate stress


	s=0.

	s1 = (lam+2*mu)*e1+lam*(e2+e3)
	s2 = (lam+2*mu)*e2+lam*(e1+e3)
	s3 = (lam+2*mu)*e3+lam*(e2+e1)
	s4 = 2*mu*e4

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Extrapolation of stresses cccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	s1= s1o + dt*s1
	s2= s2o + dt*s2
      s3= s3o + dt*s3
      s4= s4o + dt*s4

	s1o=s1
	s2o=s2
	s3o=s3
	s4o=s4


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c End of stress calculation ccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc





	contains

c --------------------------------------------------------

	subroutine pder(h,f,dx,nx,nz,dir,dim,nop)

	real h(nx,nz)
	real f(nx,nz)
!hpf$	distribute *(*,block), shadow(0,4) :: f
!hpf$   align h with f


	parameter(maxnop=8)
        real g(maxnop)
	integer dir,dim
	real dx

c derivative operators

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

	if(nop==6)then
	g(1)=-0.00468750000000
     	g(2)= 0.06510416666667
        g(3)=-1.17187500000000
	endif

        if(nop==8)then
        g(1)=0.1344048E-02
        g(2)=-0.1324718E-01
        g(3)=0.8948020E-01
        g(4)=-0.1211563E+01
        endif

c the symmetry

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

	if (dim == 1) then
	if ( dir == -1 )then

		if(nop==2)then
		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=(f(i,j)-f(i-1,j))/dx	
		end forall
		endif		

		if(nop==4)then
		forall (i=3:nx-3,j=3:nz-3)
		h(i,j)=(g(1)*f(i-2,j)+ 
     &		        g(2)*f(i-1,j)+ 
     &		        g(3)*f(i  ,j)+ 
     &		        g(4)*f(i+1,j))/dx 
		end forall
		endif		

		if(nop==8)then
		forall (i=5:nx-5,j=5:nz-5)
		h(i,j)=(g(1)*f(i-4,j)+ 
     &		        g(2)*f(i-3,j)+ 
     &		        g(3)*f(i-2,j)+ 
     &		        g(4)*f(i-1,j)+ 
     &		        g(5)*f(i  ,j)+ 
     &		        g(6)*f(i+1,j)+ 
     &		        g(7)*f(i+2,j)+ 
     &		        g(8)*f(i+3,j))/dx 
		end forall
		endif		

	else

		if(nop==2)then
		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=(f(i+1,j)-f(i,j))/dx	
		end forall
		endif


		if(nop==4)then
		forall (i=3:nx-3,j=3:nz-3)
		h(i,j)=(g(1)*f(i-1,j)+ 
     &		        g(2)*f(i,j)+ 
     &		        g(3)*f(i+1,j)+ 
     &		        g(4)*f(i+2,j))/dx 
		end forall
		endif		

		if(nop==8)then
		forall (i=5:nx-5,j=5:nz-5)
		h(i,j)=(g(1)*f(i-3,j)+ 
     &		        g(2)*f(i-2,j)+ 
     &		        g(3)*f(i-1,j)+ 
     &		        g(4)*f(i  ,j)+ 
     &		        g(5)*f(i+1,j)+ 
     &		        g(6)*f(i+2,j)+ 
     &		        g(7)*f(i+3,j)+ 
     &		        g(8)*f(i+4,j))/dx 
		end forall
		endif		

	endif

	else

	if ( dir == -1 )then

		if(nop==2)then
		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=(f(i,j)-f(i,j-1))/dx	
		end forall
		endif		

		if(nop==4)then
		forall (i=3:nx-3,j=3:nz-3)
		h(i,j)=(g(1)*f(i,j-2)+ 
     &		        g(2)*f(i,j-1)+ 
     &		        g(3)*f(i,j)+ 
     &		        g(4)*f(i,j+1))/dx 
		end forall
		endif		

		if(nop==8)then
		forall (i=5:nx-5,j=5:nz-5)
		h(i,j)=(g(1)*f(i,j-4)+ 
     &		        g(2)*f(i,j-3)+ 
     &		        g(3)*f(i,j-2)+ 
     &		        g(4)*f(i,j-1)+ 
     &		        g(5)*f(i,j)+ 
     &		        g(6)*f(i,j+1)+ 
     &		        g(7)*f(i,j+2)+ 
     &		        g(8)*f(i,j+3))/dx 
		end forall
		endif		

	else

		if(nop==2)then
		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=(f(i,j+1)-f(i,j))/dx	
		end forall
		endif


		if(nop==4)then
		forall (i=3:nx-3,j=3:nz-3)
		h(i,j)=(g(1)*f(i,j-1)+ 
     &		        g(2)*f(i,j)+ 
     &		        g(3)*f(i,j+1)+ 
     &		        g(4)*f(i,j+2))/dx 
		end forall
		endif		

		if(nop==8)then
		forall (i=5:nx-5,j=5:nz-5)
		h(i,j)=(g(1)*f(i,j-3)+ 
     &		        g(2)*f(i,j-2)+ 
     &		        g(3)*f(i,j-1)+ 
     &		        g(4)*f(i,j)+ 
     &		        g(5)*f(i,j+1)+ 
     &		        g(6)*f(i,j+2)+ 
     &		        g(7)*f(i,j+3)+ 
     &		        g(8)*f(i,j+4))/dx 
		end forall
		endif		
	endif
	endif
	

        end subroutine pder


c --------------------------------------------------------

	subroutine inter(h,f,nx,nz,dir,dim,nop)

	real h(nx,nz)
	real f(nx,nz)
!hpf$	distribute *(*,block), shadow (0,4) :: f
!hpf$   align h with f

	integer dir,dim
	real dx

	if (dim == 1) then
	   if ( dir == -1 )then

		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=0.5*(f(i,j)+f(i-1,j))
		end forall

	else

		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=0.5*(f(i+1,j)+f(i,j))
		end forall

	endif

	else

	if ( dir == -1 )then

		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=0.5*(f(i,j)+f(i,j-1))
		end forall

	else

		forall (i=2:nx-1,j=2:nz-1)
		h(i,j)=0.5*(f(i,j+1)+f(i,j))
		end forall

	endif

	endif
	

        end subroutine inter

	end subroutine psv_evolution	
