
	subroutine psv_input

	include 'common.h'

	character*15 junk
	character*50 file

	real, dimension(101) :: rch,much,lamch,rhoch

c file input
	
	Rmax=6371000.

	read(*,*)junk
	read(*,1)junk,seisfile
	read(*,1)junk,recfile
	read(*,1)junk,snapfile
	read(*,*)junk
	read(*,2)junk,nt
	read(*,3)junk,dt
	read(*,*)junk
	read(*,2)junk,izfree
	read(*,2)junk,iord
	read(*,2)junk,nop
	read(*,3)junk,aabsorb
	read(*,2)junk,nabsorb
	read(*,3)junk,Dmax
	read(*,3)junk,Thetamax
	read(*,*)junk
	read(*,2)junk,model_type
	read(*,3)junk,vp0
	read(*,3)junk,vs0
	read(*,3)junk,rho0
	read(*,*)junk
	read(*,2)junk,isx
	read(*,3)junk,sdepth
	read(*,3)junk,aa
	read(*,2)junk,it0
	read(*,3)junk,alpha
	read(*,2)junk,source_type
	read(*,*)junk
	read(*,2)junk,idisp
	read(*,2)junk,isnap
	read(*,2)junk,icolor
	read(*,2)junk,icomp
	read(*,*)junk
	read(*,2)junk,ixa
	read(*,2)junk,ixe
	read(*,2)junk,irec
	read(*,2)junk,izrec
	read(*,2)junk,isamp
	read(*,*)junk
	read(*,2)junk,ifork
	read(*,2)junk,icheck
	read(*,2)junk,icheckpoint
	read(*,3)junk,rmax
 1	format(a,a)
 2	format(a,i8)
 3	format(a,f10.6)

c read receiver coords

	open(unit=10,file=recfile)
	read(10,*)nr
	do i=1,nr
	read(10,*)rec_loc(i)
	enddo
	close(10)

c the space grid

	dr=Dmax/(nz-2*izfree)
	dphi=(Thetamax/180*pi)/nx


c initialize Radius

	forall(iz=1:nz)vr2(iz)=Rmax-(iz-izfree)*dr
	forall(iz=1:nz)vr1(iz)=vr2(iz)-.5*dr

	forall(ix=1:nx,iz=1:nz)r1(ix,iz)=vr1(iz)
	forall(ix=1:nx,iz=1:nz)r2(ix,iz)=vr2(iz)


c read model

	if(model_type==1)then

	vrho=rho0
	vmu=vrho*vs0**2
	vlam=vrho*vp0**2-2*vmu

        forall(j=1:nz)rho(:,j)=vrho(j)
        forall(j=1:nz)mu(:,j)=vmu(j)
        forall(j=1:nz)lam(:,j)=vlam(j)


	elseif(model_type==2)then


	call prem

	vrho=vrho*1000.		
	vmu=vmu*1000.
	vlam=vlam*1000.

	vmu=vrho*vmu**2
	vlam=vrho*vlam**2-2*vmu

        forall(j=1:nz)rho(:,j)=vrho(j)
	forall(j=1:nz)mu(:,j)=vmu(j)
	forall(j=1:nz)lam(:,j)=vlam(j)

	elseif(model_type==3)then

	call homo

	vrho=vrho*1000.		
	vmu=vmu*1000.
	vlam=vlam*1000.

	vmu=vrho*vmu**2
	vlam=vrho*vlam**2-2*vmu

	forall(j=1:nz)rho(:,j)=vrho(j)	
	forall(j=1:nz)mu(:,j)=vmu(j)	
	forall(j=1:nz)lam(:,j)=vlam(j)	

	elseif(model_type==4)then

	nzch=101
	open(unit=10,file='DATAV/homo100p_model')
	do i=101,1,-1
	read(10,*)rch(i),rhoch(i),lamch(i),much(i)
	enddo
	close(10)

	write(*,*)' Interpolating ... '
	do i=1,nz
c	write(*,*)' Interpolating at r = ',vr1(i)/1000,'km'
	if(vr1(i).gt.6371000.)then
	vrho(i)=rhoch(101)
	vlam(i)=lamch(101)
	vmu(i)=much(101)
	else
	call rint(rch,nzch,vr1(i),rhoch,vrho(i))
	call rint(rch,nzch,vr1(i),much,vmu(i))
	call rint(rch,nzch,vr1(i),lamch,vlam(i))
	endif
	enddo

	forall(j=1:nz)rho(:,j)=vrho(j)	
	forall(j=1:nz)mu(:,j)=vmu(j)	
	forall(j=1:nz)lam(:,j)=vlam(j)	

	open(unit=10,file=seisfile(1:lnblnk(seisfile))//'_model')
	do i=1,nz
	write(10,*)vr1(i),vrho(i),vlam(i),vmu(i)
	enddo
	close(10)


	elseif(model_type==5)then

c acoustic model

	call premnoc

	vmu=0.
	vrho=vrho*1000.		
	vmu=vmu*1000.
	vlam=vlam*1000.

	vmu=vrho*vmu**2
	vlam=vrho*vlam**2-2*vmu

        forall(j=1:nz)rho(:,j)=vrho(j)
	forall(j=1:nz)mu(:,j)=vmu(j)
	forall(j=1:nz)lam(:,j)=vlam(j)

	mu=0.

	endif

	call fd_absorb_model

	end

c -------------------------------------------
        subroutine add_num_to_char(file,i)
c -------------------------------------------


        implicit none

        character*50 file
        character*4 ichar

        integer i,in,lnblnk

        if(i.lt.10)then
        write(ichar,'(i1)') i
        in=1
        elseif(i.lt.100)then
        write(ichar,'(i2)') i
        in=2
        elseif(i.lt.1000)then
        write(ichar,'(i3)') i
        in=3
        elseif(i.lt.10000)then
        write(ichar,'(i4)') i
        in=4
        else
        write(*,*)' Warning : number too large !!! '
        endif
 
        file = file(1:lnblnk(file))//ichar(1:in)
 
        return
        end

c -----------------------------------
c
c progam prem to generate radially symmetric earth models
c
c izfree included  4-8-94

	subroutine prem

	include 'common.h'

	ddr=0.001*dmax/(nz-2*izfree)
	re=6371.0

	write(*,*)' Generating Earth model .... '
	write(*,*)' No. of radii  :  ',nz
	write(*,*)' izfree     : ',izfree
	write(*,*)' Max depth :      ',dmax
	write(*,*)' '

	do i=1,nz

	depth=(i-izfree)*ddr
	dist=re-depth
	x=dist/re




c	if(dist.lt.6371.and.dist.ge.6368)then
c	vrho(i)=1.020
c	vlam(i) =1.450
c	vmu(i) =0.
	if(dist.ge.6356)then
	vrho(i)=2.6
	vlam(i)=5.8
	vmu(i)=3.2
	elseif(dist.lt.6356.and.dist.ge.6346.6)then 
	vrho(i)=2.9
	vlam(i)=6.8
	vmu(i)=3.9
	elseif(dist.lt.6346.6.and.dist.ge.6151)then 
	vrho(i)=2.691+.6924*x
	vlam(i)=4.1875+3.9382*x
	vmu(i)=2.1519+2.3481*x
	elseif(dist.lt.6151.and.dist.ge.5971)then 
	vrho(i)=7.1089-3.8045*x
	vlam(i)=20.3926-12.2569*x
	vmu(i)=8.9496-4.4597*x
	elseif(dist.ge.5771)then 
	vrho(i)=11.2494-8.0298*x
	vlam(i)=39.7027-32.6166*x
	vmu(i)=22.3512-18.5856*x
	elseif(dist.lt.5771.and.dist.ge.5701)then 
	vrho(i)=5.3197-1.4836*x
	vlam(i)=19.0957-9.8672*x
	vmu(i)=9.9839-4.9324*x
	elseif(dist.lt.5701.and.dist.ge.5600)then 
	vrho(i)=7.9565-6.4761*x+5.5283*x**2-3.0807*x**3
	vlam(i)=29.2766-23.6027*x+5.5242*x**2-2.5514*x**3
	vmu(i)=22.3459-17.2473*x-2.0834*x**2+0.9783*x**3
	elseif(dist.lt.5600.and.dist.ge.3630)then 
	vrho(i)=7.9565-6.4761*x+5.5283*x**2-3.0807*x**3
	vlam(i)=24.9520-40.4673*x+51.4832*x**2-26.6419*x**3
	vmu(i)=11.1671-13.7818*x+17.4575*x**2-9.2777*x**3
	elseif(dist.lt.3630.and.dist.ge.3480)then 
	vrho(i)=7.9565-6.4761*x+5.5283*x**2-3.0807*x**3 
	vlam(i)=15.3891-5.3181*x+5.5242*x**2-2.5514*x**3
	vmu(i)=6.9254+1.4672*x-2.0834*x**2+.9783*x**3
	elseif(dist.lt.3480.and.dist.ge.1221.5)then 
	vrho(i)=12.5815-1.2638*x-3.6426*x**2-5.5281*x**3
	vlam(i)=11.0487-4.0362*x+4.8023*x**2-13.5732*x**3
	vmu(i)=0.
	elseif(dist.lt.1221.5)then
	vrho(i)=13.0885-8.8381*x**2
	vlam(i)=11.2622-6.3640*x**2
	vmu(i)=3.6678-4.4475*x**2
	endif

	enddo


	do i=nz-izfree,nz
	vmu(i)=vmu(nz-izfree-1)
	vrho(i)=vrho(nz-izfree-1)
	vlam(i)=vlam(nz-izfree-1)
	enddo

	end
c --------------------------------------------------
c -----------------------------------
c
c progam prem to generate radially symmetric earth models
c
c izfree included  4-8-94

	subroutine premnoc

	include 'common.h'

	ddr=0.001*dmax/(nz-2*izfree)
	re=6371.0

	write(*,*)' Generating Earth model .... '
	write(*,*)' No. of radii  :  ',nz
	write(*,*)' izfree     : ',izfree
	write(*,*)' Max depth :      ',dmax
	write(*,*)' '

	do i=1,nz

	depth=(i-izfree)*ddr
	dist=re-depth
	x=dist/re




c	if(dist.lt.6371.and.dist.ge.6368)then
c	vrho(i)=1.020
c	vlam(i) =1.450
c	vmu(i) =0.
c	if(dist.ge.6356)then
c	vrho(i)=2.6
c	vlam(i)=5.8
c	vmu(i)=3.2
c	elseif(dist.lt.6356.and.dist.ge.6346.6)then 
c	vrho(i)=2.9
c	vlam(i)=6.8
c	vmu(i)=3.9
c	elseif(dist.lt.6346.6.and.dist.ge.6151)then 
	if(dist.ge.6151)then 
	vrho(i)=2.691+.6924*x
	vlam(i)=4.1875+3.9382*x
	vmu(i)=2.1519+2.3481*x
	elseif(dist.lt.6151.and.dist.ge.5971)then 
	vrho(i)=7.1089-3.8045*x
	vlam(i)=20.3926-12.2569*x
	vmu(i)=8.9496-4.4597*x
	elseif(dist.ge.5771)then 
	vrho(i)=11.2494-8.0298*x
	vlam(i)=39.7027-32.6166*x
	vmu(i)=22.3512-18.5856*x
	elseif(dist.lt.5771.and.dist.ge.5701)then 
	vrho(i)=5.3197-1.4836*x
	vlam(i)=19.0957-9.8672*x
	vmu(i)=9.9839-4.9324*x
	elseif(dist.lt.5701.and.dist.ge.5600)then 
	vrho(i)=7.9565-6.4761*x+5.5283*x**2-3.0807*x**3
	vlam(i)=29.2766-23.6027*x+5.5242*x**2-2.5514*x**3
	vmu(i)=22.3459-17.2473*x-2.0834*x**2+0.9783*x**3
	elseif(dist.lt.5600.and.dist.ge.3630)then 
	vrho(i)=7.9565-6.4761*x+5.5283*x**2-3.0807*x**3
	vlam(i)=24.9520-40.4673*x+51.4832*x**2-26.6419*x**3
	vmu(i)=11.1671-13.7818*x+17.4575*x**2-9.2777*x**3
	elseif(dist.lt.3630.and.dist.ge.3480)then 
	vrho(i)=7.9565-6.4761*x+5.5283*x**2-3.0807*x**3 
	vlam(i)=15.3891-5.3181*x+5.5242*x**2-2.5514*x**3
	vmu(i)=6.9254+1.4672*x-2.0834*x**2+.9783*x**3
	elseif(dist.lt.3480.and.dist.ge.1221.5)then 
	vrho(i)=12.5815-1.2638*x-3.6426*x**2-5.5281*x**3
	vlam(i)=11.0487-4.0362*x+4.8023*x**2-13.5732*x**3
	vmu(i)=0.
	elseif(dist.lt.1221.5)then
	vrho(i)=13.0885-8.8381*x**2
	vlam(i)=11.2622-6.3640*x**2
	vmu(i)=3.6678-4.4475*x**2
	endif

	enddo


	do i=nz-izfree,nz
	vmu(i)=vmu(nz-izfree-1)
	vrho(i)=vrho(nz-izfree-1)
	vlam(i)=vlam(nz-izfree-1)
	enddo

	end
c --------------------------------------------------
 
c -----------------------------------
c
c progam prem to generate radially symmetric earth models
c
c izfree included  4-8-94

	subroutine homo

	include 'common.h'

	ddr=0.001*dmax/(nz-2*izfree)
	re=6371.0

	write(*,*)' Generating Earth model .... '
	write(*,*)' No. of radii  :  ',nz
	write(*,*)' izfree     : ',izfree
	write(*,*)' Max depth :      ',dmax
	write(*,*)' '

	do i=1,nz

	depth=(i-izfree)*ddr
	dist=re-depth
	x=dist/re




	if(dist.ge.3480.)then 
	vrho(i)=5.
	vlam(i)=13.163
	vmu(i)=6.5
	else
	vrho(i)=9.5
	vlam(i)=8.06
	vmu(i)=0.
	endif

	enddo


	do i=nz-izfree,nz
	vmu(i)=vmu(nz-izfree-1)
	vrho(i)=vrho(nz-izfree-1)
	vlam(i)=vlam(nz-izfree-1)
	enddo

	end
c --------------------------------------------------
 
c -----------------------------------
c
c progam prem to generate radially symmetric earth models
c
c izfree included  4-8-94

	subroutine prem_mhw

	include 'common.h'

	real md(71)
	real mvp(71)
	real mvs(71)
	real mrho(71)


	ddr=0.001*dmax/(nz-2*izfree)
	re=6371.0

	write(*,*)' Read model from file '
	write(*,*)' '
	open(unit=49,file='MODEL/model_mhw',status='unknown')
	read(49,*)
	do i=1,71
	read(49,*)md(i),mvp(i),mvs(i),mrho(i)
c	write(*,*)md(i),mvp(i),mvs(i),mrho(i)
	enddo
	close(unit=49)

	do i=1,nz

	depth=(i-izfree)*ddr


		if(depth.le.md(1))then
		vlam(i)=mvp(1)
		vmu(i)=mvs(1)
		vrho(i)=mrho(1)
		endif

		if(depth.ge.md(1))then
		
		  do j=1,70
		    if(depth.gt.md(j).and.depth.lt.md(j+1))then
		    ddepth=depth-md(j)
		    rmdepth=md(j+1)-md(j)
		    vlam(i)=mvp(j)+(mvp(j+1)-mvp(j))/rmdepth*ddepth 
		    vmu(i)=mvs(j)+(mvs(j+1)-mvs(j))/rmdepth*ddepth 
		    vrho(i)=mrho(j)+(mrho(j+1)-mrho(j))/rmdepth*ddepth 
		    goto 99
		    endif
		  enddo
                
 99		continue
		endif	  	

	enddo

	do i=nz-izfree,nz
	vmu(i)=vmu(nz-izfree-1)
	vrho(i)=vrho(nz-izfree-1)
	vlam(i)=vlam(nz-izfree-1)
	enddo

	end
c --------------------------------------------------
c -----------------------------------
c
c progam prem to generate radially symmetric earth models
c
c izfree included  4-8-94

	subroutine prem_mhw_pwdk(dddp)

	include 'common.h'

	real md(71)
	real mvp(71)
	real mvs(71)
	real mrho(71)

c	write(*,*)' Upper boundary of DDP ',dddp


	ddr=0.001*dmax/(nz-2*izfree)
	re=6371.0

c	write(*,*)' Read model from file '
c	write(*,*)' '
	open(unit=49,file='MODEL/model_mhw',status='unknown')
	read(49,*)
	do i=1,71
	read(49,*)md(i),mvp(i),mvs(i),mrho(i)
c	write(*,*)md(i),mvp(i),mvs(i),mrho(i)
	enddo
	close(unit=49)

	do i=1,nz

	depth=(i-izfree)*ddr


		if(depth.le.md(1))then
		vlam(i)=mvp(1)
		vmu(i)=mvs(1)
		vrho(i)=mrho(1)
		endif

		if(depth.ge.md(1))then
		
		  do j=1,70
		    if(depth.gt.md(j).and.depth.lt.md(j+1))then
		    ddepth=depth-md(j)
		    rmdepth=md(j+1)-md(j)
		    vlam(i)=mvp(j)+(mvp(j+1)-mvp(j))/rmdepth*ddepth 
		    vmu(i)=mvs(j)+(mvs(j+1)-mvs(j))/rmdepth*ddepth 
		    vrho(i)=mrho(j)+(mrho(j+1)-mrho(j))/rmdepth*ddepth 
		    goto 99
		    endif
		  enddo
                
 99		continue
		endif	  	

	enddo

	do i=1,nz
	depth=(i-izfree)*ddr

c setting the DDP top according to dddp

	if(depth>=2300.and.depth<=dddp)then
	vlam(i)=13.163+(depth-2300)*(13.257-13.143)/(2600.-2300.)
	elseif(depth>dddp.and.depth<=2891)then
	vlam(i)=13.714-(2891.-depth)*(13.714-13.657)/(2891.-2600.)
	endif

	enddo

	do i=nz-izfree,nz
	vmu(i)=vmu(nz-izfree-1)
	vrho(i)=vrho(nz-izfree-1)
	vlam(i)=vlam(nz-izfree-1)
	enddo

	end

C ------------------------------------------------

	subroutine fd_absorb_model

	include 'common.h'
	real taper(nx)

c taper parameters

	nabs=nabsorb

	aabs=3.*nabs
	
c initialise spatial window

	window=1.
	do i=1,nabs
	taper(nabs+1-i)=exp(-1./aabs**2*(i-1)**2)	
	enddo
	
	do i=1,nabs
	window(:,nz+1-i)=taper(i)
	enddo
	write(*,*)minval(window),' < window < ',maxval(window)

	return
	end


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

	subroutine  rint(x,nx,xr,f,out)

c this routine interpolates linearly a scalar field known at 
c x(nx) and y(ny) locations to an arbitrary location
c within xmax and ymax.
c The receiver locations are xr,  the scalar input
c field is plane and the function value at xr  is
c out.

	real x(nx)
	real xr,out
	real f(nx)

	do i=1,nx-1
	if(xr .ge. x(i) .and. xr .lt. x(i+1) ) then
	ilow=i
	iup =i+1
	dx=x(i+1)-x(i)
	hx=xr-x(i)
	endif
	enddo

	f1=f(ilow)
	f2=f(iup)

	out = f1 + (f2-f1)/dx*hx

	return
	end

