Compare commits

...

5 Commits
master ... ppeg

  1. 2
      .gitignore
  2. 7
      src/Makefile.am
  3. 47
      src/cbanddpasss.f
  4. 37
      src/cbutterworth.f
  5. 2
      src/fsimpson.f
  6. 28
      src/legendre.f
  7. 22
      src/qpalloc.f
  8. 6
      src/qpdifmat0.f
  9. 7
      src/qpdifmatl.f
  10. 24
      src/qpdifmats.f
  11. 241
      src/qpfftinv.f
  12. 20
      src/qpgetinp.f
  13. 234
      src/qpgrnspec.f
  14. 10
      src/qpmain.f
  15. 2
      src/qppsvkern.f
  16. 6
      src/qppsvkerng.f
  17. 187
      src/qpsprop.f
  18. 79
      src/qpsprop0.f
  19. 1949
      src/qpsprop1.f
  20. 61
      src/qpspropg.f
  21. 80
      src/qpspropg0.f
  22. 2459
      src/qpspropg1.f
  23. 2
      src/qpstart2t.f
  24. 5
      src/qpstart4.f
  25. 2
      src/qpstart4a.f
  26. 27
      src/qpsublayer.f
  27. 5
      src/qptprop.f
  28. 87
      src/qpwvint.f
  29. 25
      src/ruku.f
  30. 24
      src/swavelet.f
  31. 56
      src/transfs2t.f

2
.gitignore

@ -12,4 +12,4 @@
/config.status
*.o
*.mod
/src/fomosto_qssp2017
/src/fomosto_qsspppeg2017

7
src/Makefile.am

@ -1,11 +1,12 @@
bin_PROGRAMS = fomosto_qssp2017
fomosto_qssp2017_SOURCES = banddpasss.f brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f disazi.f \
bin_PROGRAMS = fomosto_qsspppeg2017
fomosto_qsspppeg2017_SOURCES = banddpasss.f brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f disazi.f \
four1w.f fsimpson.f moments.f qpalloc.f qpdifmat0.f qpdifmatl.f qpdifmats.f \
qpfftinv.f qpgetinp.f qpgrnspec.f qplegendre.f qpmain.f qppsvkern.f qppsvkerng.f \
qpqmodel.f qpshkern.f qpsmat0.f qpsmatc.f qpsmat.f qpsprop0.f qpsprop.f qpspropg0.f \
qpspropg.f qpstart0a.f qpstart0.f qpstart0g.f qpstart2t.f qpstart4a.f qpstart4.f \
qpstart4g.f qpstart6a.f qpstart6.f qpstart6g.f qpsublayer.f qptmat.f qptprop.f \
qpwvint.f ruku.f skipdoc.f spbdphj.f spbdphy.f spbh12.f spbjh.f spbphj.f spbphy.f \
spbpsj.f spbpsy.f swavelet.f taper.f transfs2t.f
spbpsj.f spbpsy.f swavelet.f taper.f transfs2t.f cbanddpasss.f cbutterworth.f \
qpsprop1.f qpspropg1.f legendre.f
clean-local:
-rm -f *.mod

47
src/cbanddpasss.f

@ -0,0 +1,47 @@
subroutine cbandpass(n,f1,f2,df,fi,nf,h)
implicit none
integer*4 n,nf
real*8 f1,f2,df,fi
complex*16 h(nf)
c
integer*4 l,k
real*8 arg
complex*16 cw,cw1,cw2,cdelt
c
real*8 PI
complex*16 ci
data PI/3.14159265358979d0/
data ci/(0.d0,1.d0)/
c
if(n.le.0.or.f1.ge.f2)then
do l=1,nf
h(l)=(1.d0,0.d0)
enddo
else if(f1.le.0.d0)then
cw2=dcmplx(-2.d0*PI*fi,2.d0*PI*f2)
do l=1,nf
cw=dcmplx(-2.d0*PI*fi,2.d0*PI*dble(l-1)*df)
cdelt=cw2-cw1
h(l)=(1.d0,0.d0)
do k=1,n
arg=PI*dble(2*k-1)/dble(2*n)
h(l)=-h(l)*ci*cdelt
& /(cw-cdelt*cdexp(dcmplx(0.d0,arg)))
enddo
enddo
else
cw1=dcmplx(-2.d0*PI*fi,2.d0*PI*f1)
cw2=dcmplx(-2.d0*PI*fi,2.d0*PI*f2)
do l=1,nf
cw=dcmplx(-2.d0*PI*fi,2.d0*PI*dble(l-1)*df)
cdelt=cw*(cw2-cw1)
h(l)=(1.d0,0.d0)
do k=1,n
arg=PI*dble(2*k-1)/dble(2*n)
h(l)=-h(l)*ci*cdelt
& /(cw*cw-cw1*cw2-cdelt*cdexp(dcmplx(0.d0,arg)))
enddo
enddo
endif
return
end

37
src/cbutterworth.f

@ -0,0 +1,37 @@
subroutine cbutterworth(n,fc,df,fi,nf,lpf)
implicit none
integer*4 n,nf
real*8 fc,df,fi
complex*16 lpf(nf)
c
integer*4 l,k
complex*16 s,cs
c
real*8 PI
data PI/3.14159265358979d0/
c
if(n.le.0.or.fc.le.0.d0)then
do l=1,nf
lpf(l)=(1.d0,0.d0)
enddo
else if(n.eq.2*(n/2))then
do l=1,nf
s=dcmplx(-fi/fc,df*dble(l-1)/fc)
lpf(l)=(1.d0,0.d0)
do k=1,n/2
cs=dcmplx(2.d0*dcos(PI*dble(2*k+n-1)/dble(2*n)),0.d0)
lpf(l)=lpf(l)/(s*s-s*cs+(1.d0,0.d0))
enddo
enddo
else
do l=1,nf
s=dcmplx(-fi/fc,df*dble(l-1)/fc)
lpf(l)=(1.d0,0.d0)/(s+(1.d0,0.d0))
do k=1,(n-1)/2
cs=dcmplx(2.d0*dcos(PI*dble(2*k+n-1)/dble(2*n)),0.d0)
lpf(l)=lpf(l)/(s*s-s*cs+(1.d0,0.d0))
enddo
enddo
endif
return
end

2
src/fsimpson.f

@ -11,8 +11,6 @@ c
real*8 eps
data nmax/2048/
data eps/1.0d-03/
c
dx=(xb-xa)/dble(n)
c
n=1
dx=xb-xa

28
src/legendre.f

@ -0,0 +1,28 @@
subroutine legendre(raddis,plm,ldeg,ldegmax)
implicit none
integer*4 ldeg,ldegmax
real*8 raddis
real*8 plm(0:ldegmax,0:2)
c
c calculate Plm(l,m,x)/(1-x^2)^(m/2)
c where Plm are the associated Legendre polynomials
c
integer*4 l,m
real*8 x
c
x=dcos(raddis)
plm(0,0)=1.d0
plm(0,1)=0.d0
plm(1,1)=1.d0
plm(0,2)=0.d0
plm(1,2)=0.d0
plm(2,2)=3.d0
do m=0,2
plm(m+1,m)=dble(2*m+1)*x*plm(m,m)
do l=m+2,ldeg
plm(l,m)=(dble(2*l-1)*x*plm(l-1,m)
& -dble(l+m-1)*plm(l-2,m))/dble(l-m)
enddo
enddo
return
end

22
src/qpalloc.f

@ -3,13 +3,13 @@ c
c CONSTANTS
c =========
integer*4 ndmax
parameter(ndmax=2)
parameter(ndmax=4)
real*8 PI,PI2
parameter(PI=3.14159265358979d0,PI2=6.28318530717959d0)
real*8 DEG2RAD,KM2M
parameter(DEG2RAD=1.745329251994328d-02,KM2M=1.0d+03)
real*8 BIGG
parameter(BIGG=6.6732d-11)
real*8 BIGG0
parameter(BIGG0=6.6732d-11)
real*8 RESOLUT
parameter(RESOLUT=0.01d0)
real*8 FLTAPER
@ -17,20 +17,20 @@ c =========
real*8 FSBUP,FSBLW,FSBREF
parameter(FSBUP=1.0d+02,FSBLW=2.5d-04,FSBREF=1.0d+00)
real*8 fbvatm,fbvocean,fbvcore
parameter(fbvatm=1.0d-06,fbvocean=0.0d+00,fbvcore=0.0d+00)
parameter(fbvatm=0.d0,fbvocean=0.0d+00,fbvcore=0.0d+00)
c
logical*2 selpsv,selsh,setzh12a,setzh12t
logical*2 nogravity,freesurf
logical*2 dispersion
c
integer*4 ngrn,nt,ntcut,ntcutout,nf,nfcut,nlpf
integer*4 ngrn,nt,ntcut,ntcutout,nf,nfcut,nbpf
integer*4 lyadd,ipatha,ipathb,ldeggr,ldegmin,ldegcut,ldegmax
integer*4 nr,igfirst,iglast
integer*4 ns,l0,lymax
integer*4 lys,lyr,lylwa,lylwb,lyos,lyob,lycm,lycc,ly0
integer*4 icmp(11)
c
real*8 dt,df,fi,fcut,fgr,rratmos,depatmos
real*8 dt,dtout,df,fi,fcut,fgr,rratmos,depatmos
real*8 rearth,rr0,minpath,maxpath
real*8 slwmax,slwlwcut,slwupcut,f1corner,f2corner
real*8 dpr,freeairgrd
@ -47,7 +47,6 @@ c
logical*2,allocatable:: ksmallp(:),ksmalls(:),ksmallt(:)
c
integer*4,allocatable:: lygrn(:),grnsel(:),ldegpsv(:),ldegsh(:)
integer*4,allocatable:: nruku(:,:)
integer*4,allocatable:: isg1(:),isg2(:),nsg(:)
integer*4,allocatable:: idr(:,:)
integer*4,allocatable:: lylwp(:),lylws(:),lylwt(:),
@ -93,7 +92,8 @@ c
& cro(:),croup(:),crolw(:),cla(:),claup(:),clalw(:),
& cmu(:),cmuup(:),cmulw(:),cvp(:),cvpup(:),cvplw(:),
& cvs(:),cvsup(:),cvslw(:),cgr(:),cgrup(:),cgrlw(:),
& cga(:),cgaup(:),cgalw(:)
& cga(:),cgaup(:),cgalw(:),cgr0(:),cgrup0(:),cgrlw0(:),
& cga0(:),cgaup0(:),cgalw0(:)
complex*16,allocatable:: cps(:,:),cpt(:,:),kp(:),ks(:),kt(:)
complex*16,allocatable:: mat2x2up(:,:,:),mat2x2lw(:,:,:),
& mat2x2inv(:,:,:),mas3x3up(:,:,:),mas3x3lw(:,:,:),
@ -115,5 +115,11 @@ c
& wspecfile(:),especfile(:),fspecfile(:),gspecfile(:),
& pspecfile(:),qspecfile(:)
character*10,allocatable:: rname(:)
c
c ppeg parameters
c
integer*4 icruku
integer*4,allocatable:: nstep(:)
real*8 BIGG
c
end module

6
src/qpdifmat0.f

@ -27,11 +27,9 @@ c
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else if(ly.ge.lyos)then
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rolw(ly)-dro*rrlw(ly)
@ -43,14 +41,14 @@ c
rorr=rolw(ly)*dexp(dlog(roup(ly)/rolw(ly))
& *(rr-rrlw(ly))/(rrup(ly)-rrlw(ly)))
crorr=dcmplx(rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
mass=2.d0*PI2/beta**3
& *(rorr*(beta*rr*(beta*rr-2.d0)+2.d0)
& -rolw(ly)*(beta*rrlw(ly)*(beta*rrlw(ly)-2.d0)+2.d0))
endif
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
& +dcmplx(BIGG0*mass/rr**2,0.d0)
c
c stabilization at the static limit
c

7
src/qpdifmatl.f

@ -30,11 +30,9 @@ c
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else if(ly.ge.lyos)then
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rorr-dro*rrlw(ly)
@ -46,14 +44,15 @@ c
rorr=rolw(ly)*dexp(dlog(roup(ly)/rolw(ly))
& *(rr-rrlw(ly))/(rrup(ly)-rrlw(ly)))
crorr=dcmplx(rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
mass=2.d0*PI2/beta**3
& *(rorr*(beta*rr*(beta*rr-2.d0)+2.d0)
& -rolw(ly)*(beta*rrlw(ly)*(beta*rrlw(ly)-2.d0)+2.d0))
endif
c
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
& +dcmplx(BIGG0*mass/rr**2,0.d0)
c
c stabilization at the static limit
c

24
src/qpdifmats.f

@ -7,8 +7,10 @@
c
c 6x6 coefficient matrix for spheroidal mode l > 0 in solid media
c
real*8 dro,rorr,ro1,mass
complex*16 cldeg,clp1,cll1,cup,clw
complex*16 crr,drr,crorr,clarr,cmurr,cxirr,cgarr,cgrrr
complex*16 cgarr0
c
complex*16 c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
@ -19,12 +21,28 @@ c
crr=dcmplx(rr,0.d0)
cup=(crr-rrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
crorr=cup*croup(ly)+clw*crolw(ly)
clarr=cup*claup(ly)+clw*clalw(ly)
cmurr=cup*cmuup(ly)+clw*cmulw(ly)
cgarr=cup*cgaup(ly)+clw*cgalw(ly)
cgrrr=cup*cgrup(ly)+clw*cgrlw(ly)
cxirr=clarr+c2*cmurr
c
if(rr.le.rrlw(ly))then
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
else
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rorr-dro*rrlw(ly)
mass=PI*(rr-rrlw(ly))*((4.d0/3.d0)*ro1
& *(rr**2+rr*rrlw(ly)+rrlw(ly)**2)
& +dro*(rr**3+rr**2*rrlw(ly)
& +rr*rrlw(ly)**2+rrlw(ly)**3))
endif
cgarr0=dcmplx(2.d0*PI2*BIGG0*rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG0*mass/rr**2,0.d0)
c
mat(1,1)=(c1-c2*clarr/cxirr)/crr
mat(1,2)=c1/cxirr/crr

241
src/qpfftinv.f

@ -4,10 +4,10 @@
integer*4 ierr
c
integer*4 lf,mf,ir
real*8 f
real*8 f,t
c
real*8,allocatable:: y0(:),y(:),dswap(:)
complex*16,allocatable:: cy(:,:),lpf(:),cswap(:)
complex*16,allocatable:: cy(:,:),bpf(:),cswap(:)
c
allocate(y0(nr),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: y0 not allocated!'
@ -15,8 +15,8 @@ c
if(ierr.ne.0)stop 'Error in qpfftinv: y not allocated!'
allocate(cy(nf,nr),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: cy not allocated!'
allocate(lpf(nf),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: lpf not allocated!'
allocate(bpf(nf),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: bpf not allocated!'
allocate(dswap(4*nf),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: dswap not allocated!'
allocate(cswap(2*nf),stat=ierr)
@ -33,41 +33,42 @@ c
c
c low-pass filter
c
if(nlpf.gt.0)then
if(nbpf.gt.0)then
if(f1corner.le.0.d0)then
call butterworth(nlpf,f2corner,df,nf,lpf)
call cbutterworth(nbpf,f2corner,df,fi,nf,bpf)
else
call bandpass(nlpf,f1corner,f2corner,df,nf,lpf)
call cbandpass(nbpf,f1corner,f2corner,df,fi,nf,bpf)
endif
c
do lf=1,nf
do ir=1,nr
ue(lf,ir)=ue(lf,ir)*lpf(lf)
un(lf,ir)=un(lf,ir)*lpf(lf)
uz(lf,ir)=uz(lf,ir)*lpf(lf)
c
ge(lf,ir)=ge(lf,ir)*lpf(lf)
gn(lf,ir)=gn(lf,ir)*lpf(lf)
gz(lf,ir)=gz(lf,ir)*lpf(lf)
c
roe(lf,ir)=roe(lf,ir)*lpf(lf)
ron(lf,ir)=ron(lf,ir)*lpf(lf)
roz(lf,ir)=roz(lf,ir)*lpf(lf)
c
uee(lf,ir)=uee(lf,ir)*lpf(lf)
uen(lf,ir)=uen(lf,ir)*lpf(lf)
uez(lf,ir)=uez(lf,ir)*lpf(lf)
unn(lf,ir)=unn(lf,ir)*lpf(lf)
unz(lf,ir)=unz(lf,ir)*lpf(lf)
uzz(lf,ir)=uzz(lf,ir)*lpf(lf)
c
see(lf,ir)=see(lf,ir)*lpf(lf)
sen(lf,ir)=sen(lf,ir)*lpf(lf)
sez(lf,ir)=sez(lf,ir)*lpf(lf)
snn(lf,ir)=snn(lf,ir)*lpf(lf)
snz(lf,ir)=snz(lf,ir)*lpf(lf)
szz(lf,ir)=szz(lf,ir)*lpf(lf)
c
gm(lf,ir)=gm(lf,ir)*lpf(lf)
ue(lf,ir)=ue(lf,ir)*bpf(lf)
un(lf,ir)=un(lf,ir)*bpf(lf)
uz(lf,ir)=uz(lf,ir)*bpf(lf)
c
ge(lf,ir)=ge(lf,ir)*bpf(lf)
gn(lf,ir)=gn(lf,ir)*bpf(lf)
gz(lf,ir)=gz(lf,ir)*bpf(lf)
c
roe(lf,ir)=roe(lf,ir)*bpf(lf)
ron(lf,ir)=ron(lf,ir)*bpf(lf)
roz(lf,ir)=roz(lf,ir)*bpf(lf)
c
uee(lf,ir)=uee(lf,ir)*bpf(lf)
uen(lf,ir)=uen(lf,ir)*bpf(lf)
uez(lf,ir)=uez(lf,ir)*bpf(lf)
unn(lf,ir)=unn(lf,ir)*bpf(lf)
unz(lf,ir)=unz(lf,ir)*bpf(lf)
uzz(lf,ir)=uzz(lf,ir)*bpf(lf)
c
see(lf,ir)=see(lf,ir)*bpf(lf)
sen(lf,ir)=sen(lf,ir)*bpf(lf)
sez(lf,ir)=sez(lf,ir)*bpf(lf)
snn(lf,ir)=snn(lf,ir)*bpf(lf)
snz(lf,ir)=snz(lf,ir)*bpf(lf)
szz(lf,ir)=szz(lf,ir)*bpf(lf)
c
gm(lf,ir)=gm(lf,ir)*bpf(lf)
enddo
enddo
endif
@ -76,148 +77,148 @@ c
c
c output seismograms: displacement vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ue,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,dispout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,un,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,dispout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,dispout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,ue,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,dispout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,un,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,dispout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,dispout(3))
endif
if(icmp(2).eq.1)then
c
c output seismograms: velocity vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ue,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,veloout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,un,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,veloout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,veloout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,ue,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,veloout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,un,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,veloout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uz,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,veloout(3))
endif
if(icmp(3).eq.1)then
c
c output seismograms: acceleration vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ue,cy,cswap,dswap,
& 2,ntcutout,rname,y0,y,acceout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,un,cy,cswap,dswap,
& 2,ntcutout,rname,y0,y,acceout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uz,cy,cswap,dswap,
& 2,ntcutout,rname,y0,y,acceout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,ue,cy,cswap,
& dswap,2,ntcutout,rname,y0,y,acceout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,un,cy,cswap,
& dswap,2,ntcutout,rname,y0,y,acceout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uz,cy,cswap,
& dswap,2,ntcutout,rname,y0,y,acceout(3))
endif
c
if(icmp(4).eq.1)then
c
c output seismograms: strain tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uee,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uen,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uez,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unn,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uzz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(6))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uee,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,strainout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uen,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,strainout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uez,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,strainout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,unn,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,strainout(4))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,unz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,strainout(5))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uzz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,strainout(6))
endif
c
if(icmp(5).eq.1)then
c
c output seismograms: strain rate tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uee,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uen,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uez,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unn,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uzz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(6))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uee,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,strainrateout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uen,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,strainrateout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uez,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,strainrateout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,unn,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,strainrateout(4))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,unz,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,strainrateout(5))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,uzz,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,strainrateout(6))
endif
if(icmp(6).eq.1)then
c
c output seismograms: stress tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,see,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sen,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sez,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snn,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,szz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(6))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,see,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,stressout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,sen,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,stressout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,sez,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,stressout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,snn,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,stressout(4))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,snz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,stressout(5))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,szz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,stressout(6))
endif
c
if(icmp(7).eq.1)then
c
c output seismograms: stress rate tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,see,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sen,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sez,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snn,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,szz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(6))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,see,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,stressrateout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,sen,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,stressrateout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,sez,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,stressrateout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,snn,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,stressrateout(4))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,snz,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,stressrateout(5))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,szz,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,stressrateout(6))
endif
c
if(icmp(8).eq.1)then
c
c output seismograms: rotation vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roe,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,rotaout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ron,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,rotaout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,rotaout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,roe,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,rotaout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,ron,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,rotaout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,roz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,rotaout(3))
endif
c
if(icmp(9).eq.1)then
c
c output seismograms: rotation rate vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roe,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,rotarateout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ron,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,rotarateout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,rotarateout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,roe,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,rotarateout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,ron,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,rotarateout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,roz,cy,cswap,
& dswap,1,ntcutout,rname,y0,y,rotarateout(3))
endif
if(icmp(10).eq.1)then
c
c output seismograms: local gravitational acceleration vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ge,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,gravout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,gn,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,gravout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,gz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,gravout(3))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,ge,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,gravout(1))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,gn,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,gravout(2))
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,gz,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,gravout(3))
endif
c
if(icmp(11).eq.1)then
c
c output seismograms: gravimeter response
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,gm,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,grmout)
call transfs2t(nf,nr,df,togsmin,dt,dtout,fi,tred,gm,cy,cswap,
& dswap,0,ntcutout,rname,y0,y,grmout)
endif
return

20
src/qpgetinp.f

@ -24,11 +24,12 @@ c time (frequency) sampling
c =========================
c
call skipdoc(unit)
read(unit,*)twindow,dt
ntcut=1+idnint(twindow/dt)
read(unit,*)twindow,dtout
ntcut=1+idnint(twindow/dtout)
nt=2
100 nt=2*nt
if(nt.lt.ntcut)goto 100
dt=twindow/dble(nt-1)
nf=nt/2
df=1.d0/(dble(nt)*dt)
c
@ -103,6 +104,12 @@ c
if(ldegcut.lt.ldegmin)ldegcut=ldegmin
ldegmin=min0(max0(1+ndmax,ldegmin),ldegcut)
ldegmax=ldegcut+ndmax+1
c
c modification for qsspppeg
c
fgr=2.d0*(dble(nf)*df+fcut)
ldeggr=2*ldegmax
selsh=.false.
c
allocate(plm(0:ldegmax,0:2),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: plm not allocated!'
@ -453,9 +460,9 @@ c
read(unit,*)outfile
call skipdoc(unit)
read(unit,*)twinout
ntcutout=min0(nt,1+idnint(twinout/dt))
ntcutout=1+idnint(twinout/dtout)
call skipdoc(unit)
read(unit,*)nlpf,f1corner,f2corner
read(unit,*)nbpf,f1corner,f2corner
call skipdoc(unit)
read(unit,*)slwlwcut,slwupcut
slwlwcut=slwlwcut/KM2M
@ -619,6 +626,11 @@ c
endif
enddo
c
if(abs(dp0(l)-rratmos)/rratmos.lt.1.0d-6)then
write(*,*) ' Snapping lowest layer depth to earth model radius.'
dp0(l)=rratmos
endif
if(dp0(l).gt.rratmos)then
stop ' Error in qpgetinp: earth radius larger than pre-defined!'
else if(dp0(l).lt.rratmos)then

234
src/qpgrnspec.f

@ -3,21 +3,23 @@
implicit none
integer*4 ig
c
integer*4 i,nn,il,istp,ly,lf,ldeg,ldeg0,ldegf,ldegup,ierr
real*8 f,ksp,xlw,xup,dll,omi,expo,rpath,slwcut
integer*4 i,nn,il,istp,ly,lf,ldeg,ldeg0,ldegf,ldegup
integer*4 ierr,nly,ncruku
real*8 f,h,ksp,xlw,xup,dll,omi,expo,rpath,slwcut
real*8 fac,fl,depst1,depst2,dys2,rr0a,a,b,rrs,x
real*8 kcut1(4),kcut2(4)
complex*16 ca,cb,cag,cs1,cs2,cs3,cs4,ct1,ct2,cll1
complex*16 cmur,ys3d,yt1d,cg1,cg5
complex*16 ypsv(6,4),ypsvg(6,4),ysh(2,2)
logical*2 forruku
complex*16 cmur,ys3d,yt1d,cg1,cg5,cgfac
complex*16 ypsv(6,4),ypsv0(6,4),ypsvg(6,4),ysh(2,2),y1aux(4)
c
real*8 fsimpson
c
real*8 expos
complex*16 c2,c3,c4
data expos/48.d0/
data c2,c3,c4/(2.d0,0.d0),(3.d0,0.d0),(4.d0,0.d0)/
real*8 expos,gfac
complex*16 c1,c2,c3,c4
data expos,gfac/24.d0,0.1d0/
data c1,c2,c3,c4/(1.d0,0.d0),(2.d0,0.d0),(3.d0,0.d0),(4.d0,0.d0)/
c
cgfac=dcmplx(gfac,0.d0)
c
if(ig.eq.igfirst)then
allocate(disk(0:ldegmax),stat=ierr)
@ -138,16 +140,28 @@ c
open(27,file=pspecfile(ig),form='unformatted',status='unknown')
open(28,file=qspecfile(ig),form='unformatted',status='unknown')
c
ldeg0=1+ndmax
do ldeg0=1+ndmax,ldegmin
dll=dsqrt(dble(ldeg0)*dabs(dble(ldeg0-1)))
expo=0.d0
do ly=min0(lys,lyr),max0(lys,lyr)-1
expo=expo+dll*dlog(rrup(ly)/rrlw(ly))
ldeg0=10+ndmax
if(vsup(lys).gt.0.d0)then
slwcut=1.d0/vsup(lys)
else
slwcut=1.d0/vpup(lys)
endif
if(vsup(lyr).gt.0.d0)then
slwcut=dmin1(slwcut,1.d0/vsup(lyr))
else
slwcut=dmin1(slwcut,1.d0/vpup(lyr))
endif
if(slwmax.ge.slwcut)then
do ldeg0=10+ndmax,ldegmin
dll=dsqrt(dble(ldeg0)*dabs(dble(ldeg0-1)))
expo=0.d0
do ly=min0(lys,lyr),max0(lys,lyr)-1
expo=expo+dll*dlog(rrup(ly)/rrlw(ly))
enddo
if(expo.gt.expos)goto 10
enddo
if(expo.gt.expos)goto 10
enddo
10 continue
10 continue
endif
c
lylwa=0
if(ipatha.eq.1)then
@ -173,7 +187,7 @@ c
30 continue
endif
c
omi=PI2*fcut
omi=PI2*dble(nfcut-1)*df
slwcut=slwmax
c
if(lylwa.gt.lys+1)then
@ -184,18 +198,12 @@ c
slwcut=dmin1(slwcut,fac/vpup(lylwa))
endif
endif
c
if(idint(rearth*omi*slwcut).gt.ldegcut)then
print *,' Warning from qpgrnspec: ldegcut may be too small!'
endif
ldegup=min0(ldegcut,max0(ldeg0,idint(rearth*omi*slwcut)))
c
if(fgr.gt.0.d0.or.ldeggr.gt.0)then
do ly=1,ly0
do ldeg=0,ldegup
nruku(ldeg,ly)=10
enddo
enddo
endif
ldegup=min0(ldegcut,max0(ldeg0,idint(rearth*omi*slwcut)))
c
write(*,*)' '
write(*,'(a,i3,a,f7.2,a)')' ... calculate Green functions for ',
@ -211,13 +219,20 @@ c
write(27)nt,ntcut,dt,nf,nfcut,df,ldegup
write(28)nt,ntcut,dt,nf,nfcut,df,ldegup
c
if(.not.nogravity)then
do ly=1,ly0
do ldeg=0,ldegup
nruku(ldeg,ly)=0
enddo
do istp=1,6
do ldeg=0,ldegmax
ul0(ldeg,istp)=(0.d0,0.d0)
vl0(ldeg,istp)=(0.d0,0.d0)
wl0(ldeg,istp)=(0.d0,0.d0)
el0(ldeg,istp)=(0.d0,0.d0)
fl0(ldeg,istp)=(0.d0,0.d0)
gl0(ldeg,istp)=(0.d0,0.d0)
pl0(ldeg,istp)=(0.d0,0.d0)
ql0(ldeg,istp)=(0.d0,0.d0)
enddo
endif
enddo
c
icruku=0
c
do lf=1,nfcut
f=dble(lf-1)*df
@ -256,10 +271,10 @@ c
do ldeg=0,ldegf
if(nogravity.or.lylwa.gt.0.or.lylwb.le.ly0.or.
& .not.freesurf)then
forruku=.false.
c forruku=.false.
else
fac=dsqrt((f/fgr)**2+(dble(ldeg)/dble(ldeggr))**2)
forruku=fac.lt.1.d0+FLTAPER
c forruku=fac.lt.1.d0+FLTAPER
endif
dll=dsqrt(dble(ldeg)*dabs(dble(ldeg-1)))
c
@ -440,12 +455,12 @@ c
if(selpsv.or.lys.lt.lyob.or.lys.ge.lycm.and.lys.lt.lycc)then
ly=max0(lylwa,min0(lylwb,lylwp(0),lylws(0)))
depst1=(rratmos-depatmos-rrup(ly))/KM2M
ly=max0(lylwa,min0(lylwb,lylwp(ldegf+1),lylws(ldegf+1)))
ly=max0(lylwa,min0(lylwb,lylwp(ldegf),lylws(ldegf)))
depst2=(rratmos-depatmos-rrup(ly))/KM2M
else
ly=max0(lylwa,min0(lylwb,lylwt(1)))
depst1=(rratmos-depatmos-rrup(ly))/KM2M
ly=max0(lylwa,min0(lylwb,lylwt(ldegf+1)))
ly=max0(lylwa,min0(lylwb,lylwt(ldegf)))
depst2=(rratmos-depatmos-rrup(ly))/KM2M
endif
c
@ -455,7 +470,7 @@ c
if(lys.ge.lyob.and.lys.le.min0(lycm-1,ly0))then
do ly=lyob,min0(lycm-1,ly0)
ldegsh(ly)=1
do ldeg=1,ldegf+1
do ldeg=1,ldegf
if(ly.ge.lyupt(ldeg).and.ly.le.lylwt(ldeg))then
ldegsh(ly)=ldeg
endif
@ -464,7 +479,7 @@ c
else if(lys.ge.lycc)then
do ly=lycc,ly0
ldegsh(ly)=1
do ldeg=1,ldegf+1
do ldeg=1,ldegf
if(ly.ge.lyupt(ldeg).and.ly.le.lylwt(ldeg))then
ldegsh(ly)=ldeg
endif
@ -477,7 +492,7 @@ c of psv solution
c
do ly=1,ly0
ldegpsv(ly)=0
do ldeg=0,ldegf+1
do ldeg=0,ldegf
if(ly.ge.min0(lyupp(ldeg),lyups(ldeg)).and.
& ly.le.max0(lylwp(ldeg),lylws(ldeg)))then
ldegpsv(ly)=ldeg
@ -490,29 +505,113 @@ c
endif
c
do ldeg=0,ldegf
c
ncruku=0
do ly=1,ly0-1
h=rrup(ly)-rrlw(ly)
nly=1+idint(h*(f/vpup(ly)+0.5d0/rrlw(ly)))
if(vpup(ly).gt.0.d0)then
nly=max0(nly,1+idint(h*dmax1(f/vpup(ly),
& dble(ldeg)/rrlw(ly)))/5)
else
nly=max0(nly,1+idint(h*dmax1(f/vsup(ly),
& dble(ldeg)/rrlw(ly)))/5)
endif
ncruku=ncruku+nly
enddo
ncruku=ncruku+1000
if(icruku.gt.0)deallocate(nstep)
allocate(nstep(ncruku),stat=ierr)
if(ierr.ne.0)stop 'Error in qpgrnspec: nstep not allocated!'
c
if(.not.selpsv.or.ldeg.gt.ldegcut.or.
& max0(lylwp(ldeg),lylws(ldeg)).lt.max0(lylwa,lyr,lys).or.
& min0(lyupp(ldeg),lyups(ldeg)).gt.min0(lyr,lys))then
do istp=1,4
y1aux(istp)=(0.d0,0.d0)
do i=1,6
ypsv(i,istp)=(0.d0,0.d0)
enddo
enddo
else if(nogravity)then
call qppsvkern(f,ldeg,ypsv)
do istp=1,4
y1aux(istp)=ypsv(1,istp)
do i=1,4
ypsv(i,istp)=(0.d0,0.d0)
enddo
enddo
else
fac=dsqrt((f/fgr)**2+(dble(ldeg)/dble(ldeggr))**2)
if(fac.le.1.d0)then
do i=1,ncruku
nstep(i)=0
enddo
icruku=0
call qppsvkerng(f,ldeg,ypsv)
icruku=0
BIGG=BIGG0*gfac
do ly=1,ly0
cgaup(ly)=cgaup0(ly)*cgfac
cgalw(ly)=cgalw0(ly)*cgfac
cga(ly)=cga0(ly)*cgfac
enddo
call qppsvkerng(f,ldeg,ypsvg)
BIGG=BIGG0
do ly=1,ly0
cgaup(ly)=cgaup0(ly)
cgalw(ly)=cgalw0(ly)
cga(ly)=cga0(ly)
enddo
do istp=1,4
y1aux(istp)=ypsv(1,istp)
do i=1,4
ypsv(i,istp)=(ypsv(i,istp)-ypsvg(i,istp))/(c1-cgfac)
enddo
enddo
else if(fac.ge.1.d0+FLTAPER)then
call qppsvkern(f,ldeg,ypsv)
do istp=1,4
y1aux(istp)=ypsv(1,istp)
do i=1,4
ypsv(i,istp)=(0.d0,0.d0)
enddo
enddo
else
do i=1,ncruku
nstep(i)=0
enddo
icruku=0
call qppsvkerng(f,ldeg,ypsv)
icruku=0
BIGG=BIGG0*gfac
do ly=1,ly0
cgaup(ly)=cgaup0(ly)*cgfac
cgalw(ly)=cgalw0(ly)*cgfac
cga(ly)=cga0(ly)*cgfac
enddo
call qppsvkerng(f,ldeg,ypsvg)
call qppsvkern(f,ldeg,ypsv)
BIGG=BIGG0
do ly=1,ly0
cgaup(ly)=cgaup0(ly)
cgalw(ly)=cgalw0(ly)
cga(ly)=cga0(ly)
enddo
do istp=1,4
y1aux(istp)=ypsv(1,istp)
do i=1,4
ypsvg(i,istp)=(ypsv(i,istp)-ypsvg(i,istp))/(c1-cgfac)
enddo
enddo
c
ca=dcmplx(dsin(0.5d0*PI*(fac-1.d0)/FLTAPER)**2,0.d0)
cb=(1.d0,0.d0)-ca
call qppsvkern(f,ldeg,ypsv)
do istp=1,4
y1aux(istp)=ca*ypsv(1,istp)+cb*y1aux(istp)
do i=1,4
ypsv(i,istp)=(0.d0,0.d0)
enddo
do i=1,6
ypsv(i,istp)=ca*ypsv(i,istp)+cb*ypsvg(i,istp)
enddo
@ -520,18 +619,23 @@ c
endif
endif
c
if(.not.selsh.or.ldeg.gt.ldegcut.or.ldeg.eq.0.or.
& lys.lt.lyob.or.lys.ge.lycm.and.lys.lt.lycc.or.
& lylwt(ldeg).lt.max0(lylwa,lyr,lys).or.
& lyupt(ldeg).gt.min0(lyr,lys))then
do istp=1,2
do i=1,2
ysh(i,istp)=(0.d0,0.d0)
enddo
c if(.not.selsh.or.ldeg.gt.ldegcut.or.ldeg.eq.0.or.
c & lys.lt.lyob.or.lys.ge.lycm.and.lys.lt.lycc.or.
c & lylwt(ldeg).lt.max0(lylwa,lyr,lys).or.
c & lyupt(ldeg).gt.min0(lyr,lys))then
c do istp=1,2
c do i=1,2
c ysh(i,istp)=(0.d0,0.d0)
c enddo
c enddo
c else
c call qpshkern(f,ldeg,ysh)
c endif
do istp=1,2
do i=1,2
ysh(i,istp)=(0.d0,0.d0)
enddo
else
call qpshkern(f,ldeg,ysh)
endif
enddo
c
if(lyr.gt.1)then
cg1=cgalw(lyr-1)
@ -560,7 +664,7 @@ c
fl0(ldeg,1)=cs2*ypsv(4,2)
gl0(ldeg,1)=(0.d0,0.d0)
pl0(ldeg,1)=cs2*ypsv(5,2)
ql0(ldeg,1)=cs2*(cg1*ypsv(1,2)-cg5*ypsv(5,2)+ypsv(6,2))
ql0(ldeg,1)=cs2*(cg1*y1aux(2)-cg5*ypsv(5,2)+ypsv(6,2))
endif
c
c 2. Explosion (M11=M22=M33=1)
@ -568,8 +672,7 @@ c
cs1=dcmplx( disk(ldeg)/roup(lys),0.d0)/cvpup(lys)**2
cs2=dcmplx(-disk(ldeg)*4.d0/rrup(lys),0.d0)
& *(cvsup(lys)/cvpup(lys))**2
cs4=dcmplx( disk(ldeg)*2.d0/rrup(lys),0.d0)
& *(cvsup(lys)/cvpup(lys))**2
cs4=-(0.5d0,0.d0)*cs2
ul0(ldeg,2)=cs1*ypsv(1,1)+cs2*ypsv(1,2)+cs4*ypsv(1,4)
vl0(ldeg,2)=cs1*ypsv(3,1)+cs2*ypsv(3,2)+cs4*ypsv(3,4)
wl0(ldeg,2)=(0.d0,0.d0)
@ -577,9 +680,9 @@ c
fl0(ldeg,2)=cs1*ypsv(4,1)+cs2*ypsv(4,2)+cs4*ypsv(4,4)
gl0(ldeg,2)=(0.d0,0.d0)
pl0(ldeg,2)=cs1*ypsv(5,1)+cs2*ypsv(5,2)+cs4*ypsv(5,4)
ql0(ldeg,2)=cs1*(cg1*ypsv(1,1)-cg5*ypsv(5,1)+ypsv(6,1))
& +cs2*(cg1*ypsv(1,2)-cg5*ypsv(5,2)+ypsv(6,2))
& +cs4*(cg1*ypsv(1,4)-cg5*ypsv(5,4)+ypsv(6,4))
ql0(ldeg,2)=cs1*(cg1*y1aux(1)-cg5*ypsv(5,1)+ypsv(6,1))
& +cs2*(cg1*y1aux(2)-cg5*ypsv(5,2)+ypsv(6,2))
& +cs4*(cg1*y1aux(4)-cg5*ypsv(5,4)+ypsv(6,4))
c
c 3. CLVD (M33=1,M11=M22=-0.5)
c
@ -604,9 +707,9 @@ c
fl0(ldeg,3)=cs1*ypsv(4,1)+cs2*ypsv(4,2)+cs4*ypsv(4,4)
gl0(ldeg,3)=(0.d0,0.d0)
pl0(ldeg,3)=cs1*ypsv(5,1)+cs2*ypsv(5,2)+cs4*ypsv(5,4)
ql0(ldeg,3)=cs1*(cg1*ypsv(1,1)-cg5*ypsv(5,1)+ypsv(6,1))
& +cs2*(cg1*ypsv(1,2)-cg5*ypsv(5,2)+ypsv(6,2))
& +cs4*(cg1*ypsv(1,4)-cg5*ypsv(5,4)+ypsv(6,4))
ql0(ldeg,3)=cs1*(cg1*y1aux(1)-cg5*ypsv(5,1)+ypsv(6,1))
& +cs2*(cg1*y1aux(2)-cg5*ypsv(5,2)+ypsv(6,2))
& +cs4*(cg1*y1aux(4)-cg5*ypsv(5,4)+ypsv(6,4))
endif
c
c 4. Horizontal single force (F1=1)
@ -621,7 +724,7 @@ c
pl0(ldeg,4)=(0.d0,0.d0)
ql0(ldeg,4)=(0.d0,0.d0)
else
cs4=dcmplx(-0.5d0*disk(ldeg)/(dble(ldeg)*dble(ldeg+1)),0.d0)
cs4=dcmplx(-disk(ldeg)/(dble(ldeg)*dble(ldeg+1)),0.d0)
ct2=cs4
ul0(ldeg,4)=cs4*ypsv(1,4)
vl0(ldeg,4)=cs4*ypsv(3,4)
@ -630,7 +733,7 @@ c
fl0(ldeg,4)=cs4*ypsv(4,4)
gl0(ldeg,4)=ct2*ysh(2,2)
pl0(ldeg,4)=cs4*ypsv(5,4)
ql0(ldeg,4)=cs4*(cg1*ypsv(1,4)-cg5*ypsv(5,4)+ypsv(6,4))
ql0(ldeg,4)=cs4*(cg1*y1aux(4)-cg5*ypsv(5,4)+ypsv(6,4))
endif
c
c 5. Dip-slip (M13=M31=1)
@ -655,7 +758,7 @@ c
fl0(ldeg,5)=cs3*ypsv(4,3)
gl0(ldeg,5)=ct1*ysh(2,1)
pl0(ldeg,5)=cs3*ypsv(5,3)
ql0(ldeg,5)=cs3*(cg1*ypsv(1,3)-cg5*ypsv(5,3)+ypsv(6,3))
ql0(ldeg,5)=cs3*(cg1*y1aux(3)-cg5*ypsv(5,3)+ypsv(6,3))
endif
c
c 6. Strike-slip (M12=M21=1)
@ -680,7 +783,7 @@ c
fl0(ldeg,6)=cs4*ypsv(4,4)
gl0(ldeg,6)=ct2*ysh(2,2)
pl0(ldeg,6)=cs4*ypsv(5,4)
ql0(ldeg,6)=cs4*(cg1*ypsv(1,4)-cg5*ypsv(5,4)+ypsv(6,4))
ql0(ldeg,6)=cs4*(cg1*y1aux(4)-cg5*ypsv(5,4)+ypsv(6,4))
endif
enddo
c
@ -725,5 +828,6 @@ c
& zh12p,zh12sv,zh12sh,
& zjupg,zjlwg,zhupg,zhlwg,wjg,whg)
endif
if(icruku.gt.0)deallocate(nstep)
return
end

10
src/qpmain.f

@ -21,11 +21,15 @@ c
print *,'# QQQQQ SSSSS SSSSS P #'
print *,'# #'
print *,'# Complete synthetic seismograms #'
print *,'# (displacement/strain/stress/rotation) #'
print *,'# (displacement/strain/stress/rotation/gravity) #'
print *,'# based on #'
print *,'# a spherically symmetric earth model #'
print *,'# #'
print *,'# (Version 2017) #'
print *,'# Last update: 2019-5-27 #'
print *,'# #'
print *,'# This version (QSSPPPEG) only provides #'
print *,'# (pre-P-wave) elasto-gravitational effect #'
print *,'# #'
print *,'# by #'