Commit b19c2fac authored by Rongjiang Wang's avatar Rongjiang Wang Committed by Sebastian Heimann
Browse files

eliminate numerical phases in late body waves > 120 deg

parent b5800fb7
......@@ -261,12 +261,11 @@ c
read(unit,*)nlpf,f1corner,f2corner
call skip_comments(unit)
read(unit,*)slwlwcut,slwupcut
if(slwlwcut.gt.slwupcut)then
slwlwcut=slwlwcut/KM2M
slwupcut=slwupcut/KM2M
if(slwupcut.gt.slwmax.or.slwlwcut.ge.slwupcut)then
slwlwcut=0.d0
slwupcut=slwmax
else
slwlwcut=slwlwcut/KM2M
slwupcut=slwupcut/KM2M
endif
call skip_comments(unit)
read(unit,*)nr
......
......@@ -25,9 +25,9 @@ c nfmax: max. number of frequency samples
c ldegmax: max. degree of Legendre polynomials
c
integer lymax,nrmax,nsmax,ngrnmax,nfmax,ldegmax,ndmax
parameter(nrmax=51,nsmax=500,ngrnmax=50)
parameter(nfmax=8192,ldegmax=20000)
parameter(lymax=250+ngrnmax)
parameter(nrmax=91,nsmax=500,ngrnmax=50)
parameter(nfmax=16384,ldegmax=25000)
parameter(lymax=150+ngrnmax)
parameter(ndmax=4)
c
c GREEN FUNCTION PARAMETERS
......
......@@ -4,7 +4,7 @@
c
include 'qpglobal.h'
c
integer i,il,istp,ly,lf,ldeg,ldeg0,ldegf,ldegup,ldegpre
integer i,il,istp,ly,lf,ldeg,ldeg0,ldegf,ldegup
double precision f,ksp2,dll1,omi,expo
double precision fl,depst1,depst2,dys2,dxs
double precision kcut1(4),kcut2(4)
......@@ -16,9 +16,9 @@ c
double complex ypsv(6,4),ypsvg(6,4),ysh(2,2)
double complex ysgr(4),yspr(4),ystt(4),yttt(4)
c
double precision expos,expoa,expof
double precision expos
double complex c2,c3,c4
data expos,expoa,expof/12.d0,12.d0,6.d0/
data expos/24.d0/
data c2,c3,c4/(2.d0,0.d0),(3.d0,0.d0),(4.d0,0.d0)/
c
c Initiation
......@@ -90,35 +90,19 @@ c
open(27,file=wgrnfile(ig),
& form='unformatted',status='unknown')
c
ldeg0=1+ndmax
do ldeg0=1+ndmax,ldegmin
ldeg0=10+ndmax
do ldeg0=10+ndmax,ldegmin
dll1=dble(ldeg0)*dble(ldeg0+1)
expo=0.d0
do ly=min0(lys,lyr),max0(lys,lyr)-1
expo=expo+dsqrt(dll1)*dlog(rrup(ly)/rrlw(ly))
enddo
if(expo.gt.expos+dlog(dsqrt(dll1)*fdisk(ldeg0)))goto 10
if(expo.gt.expos)goto 10
enddo
10 continue
c
omi=PI2*fcut
ldegup=ldeg0
do ldegup=ldeg0,min0(ldegcut+ndmax,ldeg0+idint(REARTH*omi*slwmax))
dll1=dble(ldegup)*dble(ldegup+1)
expo=0.d0
do ly=min0(lys,lyr),max0(lys,lyr)-1
if(vsup(ly).gt.0.d0)then
ksp2=(omi*rrup(ly)/vsup(ly))**2
else
ksp2=(omi*rrup(ly)/vpup(ly))**2
endif
if(dll1.gt.ksp2)then
expo=expo+dsqrt(dll1-ksp2)*dlog(rrup(ly)/rrlw(ly))
endif
enddo
if(expo.gt.expof+dlog(dsqrt(dll1)*fdisk(ldegup)))goto 20
enddo
20 continue
ldegup=min0(ldegcut+ndmax,ldeg0+idint(REARTH*omi*slwmax))
c
if(ldegup.ge.ldegmax-ndmax-1)then
print *,' Warning in qpgrnspec: cutoff degree = ',ldegup,
......@@ -152,29 +136,12 @@ c
enddo
enddo
c
ldegpre=ldeg0
do lf=1,nfcut
f=dble(lf-1)*df
omi=PI2*f
call qpqmodel(f)
c
do ldegf=ldegpre,min0(ldegup,ldeg0+idint(REARTH*omi*slwmax))
dll1=dble(ldegf)*dble(ldegf+1)
expo=0.d0
do ly=min0(lys,lyr),max0(lys,lyr)-1
if(vsup(ly).gt.0.d0)then
ksp2=(omi*rrup(ly)/vsup(ly))**2
else
ksp2=(omi*rrup(ly)/vpup(ly))**2
endif
if(dll1.gt.ksp2)then
expo=expo+dsqrt(dll1-ksp2)*dlog(rrup(ly)/rrlw(ly))
endif
enddo
if(expo.gt.expos+dlog(dsqrt(dll1)*fdisk(ldegf)))goto 50
enddo
50 ldegf=min0(ldegf,ldegup)
ldegpre=ldegf
ldegf=min0(ldegup,ldeg0+idint(REARTH*omi*slwmax))
c
do ldeg=0,ldegf+1
dll1=dble(ldeg)*dble(ldeg+1)
......@@ -211,7 +178,7 @@ c
goto 200
endif
enddo
do ly=lycm,min0(lycc,ly0)-1
do ly=max0(lys,lyr,lycm-1)+1,min0(lycc,ly0)-1
ksp2=(omi*rrup(ly)/vpup(ly))**2
if(dll1.gt.ksp2)then
expo=expo+dsqrt(dll1-ksp2)*dlog(rrup(ly)/rrlw(ly))
......@@ -221,7 +188,7 @@ c
goto 200
endif
enddo
do ly=lycc,ly0-1
do ly=max0(lys,lyr,lycc-1)+1,ly0-1
ksp2=(omi*rrup(ly)/vsup(ly))**2
if(dll1.gt.ksp2)then
expo=expo+dsqrt(dll1-ksp2)*dlog(rrup(ly)/rrlw(ly))
......
......@@ -4,63 +4,43 @@
c
include 'qpglobal.h'
c
integer i,id,is,ir,ig,nd,nt0,nf0,ntcut0,nfcut0,ly,ishift
integer lf,lf1,istp,ldeg,ldeg0,ldegf,ldegf0
integer i,id,is,ir,ig,nd,nt0,nf0,ntcut0,nfcut0,ishift
integer lf,lf1,istp,ldeg,ldegf
integer istat,ldegup,ldeglw,ldegneed
integer ldegtap(4,nsmax,nrmax)
double precision depsarc,rrs,slws,slwr,twin,anorm
double precision f,dt0,df0,rn,re,azi,bazi,bazi0,f20
double precision tap(0:ldegmax),ldf(nsmax,nrmax)
double precision f0(nsmax,nrmax)
double precision depsarc,anorm
double precision f,dt0,df0,rn,re,azi,bazi,bazi0
double precision tap(0:ldegmax)
double complex cp0,cp1,cp2,wavelet
double complex cfac,ca,cb,dur,dut,dup,dus,dug,duv,duw
double complex wvf(nfmax,nsmax)
double complex expl(nsmax),clvd(nsmax),ss12(nsmax)
double complex ss11(nsmax),ds31(nsmax),ds23(nsmax)
logical fullwave
c
do ldeg=0,ldegmax
tap(ldeg)=1.d0
enddo
c
fullwave=ldegmin.gt.10*(1+ndmax)
& .and.(nlpf.le.0.or.nlpf.gt.0.and.f1corner.le.0.d0)
c
do is=1,ns
depsarc=deps(is)/REARTH
rrs=REARTH-deps(is)
lys=lyob
do ly=lyob+1,min0(lycm-1,ly0)
if(rrs.le.rrup(ly))then
lys=ly
endif
enddo
if(slwupcut.gt.0.d0)then
slws=dmin1(slwupcut,slwmax)*REARTH
else
slws=slwmax*REARTH
endif
c
do ir=1,nr
call disazi(1.d0,lats(is),lons(is),
& latr(ir),lonr(ir),rn,re)
c
dis(is,ir)=dsqrt(rn**2+re**2)
twin=dble(ntcutout)*dt+tred(ir)-togs(is)
if(dis(is,ir).gt.0.d0)then
slwr=twin/dis(is,ir)
else
slwr=slws
endif
c
if(slwr.ge.slws.or.slwupcut.le.0.d0)then
ldf(is,ir)=0.d0
c determine order of differential transform
c
if(fullwave.or.dis(is,ir).le.dmax1(5.d0*depsarc,PI/90.d0))then
idr(is,ir)=0
else
ldf(is,ir)=PI2*slwr
idr(is,ir)=min0(ndmax,
& idnint(dlog(dis(is,ir)/depsarc)/dlog(10.d0)))
idr(is,ir)=idnint(dble(ndmax)*(dis(is,ir)-5.d0*depsarc)/PI)
endif
c
c f0 is the frequency, at which the distance is equal to 10 times
c of cut-off wavelength
c
f0(is,ir)=10.d0/(slwr*dsqrt(dis(is,ir)**2+depsarc**2))
c
ssd(is,ir)=dcmplx(dsin(dis(is,ir)),0.d0)
ssf(is,ir)=dcmplx(2.d0*dsin(0.5d0*dis(is,ir))**2,0.d0)
......@@ -186,7 +166,6 @@ c
stop
endif
c
ldegf0=0
do lf=1,nfcut
f=dble(lf-1)*df
read(21)ldegf
......@@ -205,28 +184,27 @@ c
read(26)((yv(ldeg,istp,0),ldeg=0,ldegf),istp=1,4)
read(27)((yw(ldeg,istp,0),ldeg=0,ldegf),istp=2,3)
c
if(lf.eq.1)ldegf0=ldegf
ldegneed=0
ldeglw=ldegf
ldeg0=ldegf
ldeglw=ldegmax
do is=isg1(ig),isg2(ig)
do ir=1,nr
id=idr(is,ir)
if(ldf(is,ir).le.0.d0.or.dis(is,ir).le.0.d0)then
ldegtap(4,is,ir)=max0(0,ldegf-id)
if(fullwave)then
ldegtap(4,is,ir)=max0(0,ldegf-idr(is,ir))
ldegtap(3,is,ir)=ldegtap(4,is,ir)*4/5
ldegtap(2,is,ir)=0
ldegtap(1,is,ir)=0
pause
else
ldegtap(4,is,ir)=max0(0,min0(ldegf-id,
& ldegf0+idnint(dsqrt(f**2+f0(is,ir)**2)*ldf(is,ir))))
ldegtap(4,is,ir)=idnint(REARTH*PI2*f*slwupcut)
ldegtap(3,is,ir)=ldegtap(4,is,ir)*4/5
ldegtap(4,is,ir)=min0(ldegf-idr(is,ir),
& 10+ldegtap(4,is,ir))
ldegtap(2,is,ir)=min0(ldegtap(3,is,ir),
& idnint(REARTH*PI2*f*slwlwcut))
ldegtap(1,is,ir)=ldegtap(2,is,ir)*4/5
endif
ldegtap(3,is,ir)=max0(0,min0(ldegtap(4,is,ir)*4/5,
& max0(0,ldegtap(4,is,ir)-40)))
ldegneed=max0(ldegneed,ldegtap(4,is,ir))
ldeglw=min0(ldeglw,ldegtap(4,is,ir))
ldegtap(2,is,ir)=min0(ldegtap(3,is,ir),
& idnint(REARTH*PI2*f*slwlwcut))
ldegtap(1,is,ir)=min0(ldegtap(3,is,ir),
& idnint(0.8d0*REARTH*PI2*f*slwlwcut))
ldeg0=min0(ldeg0,ldegtap(1,is,ir))
ldeglw=max0(0,min0(ldeglw,ldegtap(1,is,ir)-idr(is,ir)))
enddo
enddo
ldegneed=min0(ldegneed+nd,ldegf)
......
......@@ -9,11 +9,11 @@ c
double precision PI
data PI/3.14159265358979d0/
c
do l=0,ldeg(1)
do l=0,ldeg(1)-1
tap(l)=0.d0
enddo
fac=0.5d0*PI/dble(1+ldeg(2)-ldeg(1))
do l=ldeg(1)+1,ldeg(2)
do l=ldeg(1),ldeg(2)-1
tap(l)=dsin(fac*dble(l-ldeg(1)))**2
enddo
do l=ldeg(2),ldeg(3)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment