Compare commits

...

6 Commits

  1. 2
      configure.ac
  2. 4
      doc/examples/qssp2017-test.inp
  3. 8
      src/Makefile.am
  4. 2
      src/fsimpson.f
  5. 6
      src/legendre.f
  6. 16
      src/qpalloc.f
  7. 27
      src/qpdifmats.f
  8. 256
      src/qpfftinv.f
  9. 37
      src/qpgetinp.f
  10. 63
      src/qpgrnspec.f
  11. 5
      src/qpmain.f
  12. 2
      src/qppsvkern.f
  13. 46
      src/qpqmodel.f
  14. 187
      src/qpsprop.f
  15. 79
      src/qpsprop0.f
  16. 25
      src/qpspropg.f
  17. 77
      src/qpspropg0.f
  18. 2
      src/qpstart2t.f
  19. 5
      src/qpstart4.f
  20. 2
      src/qpstart4a.f
  21. 5
      src/qptprop.f
  22. 87
      src/qpwvint.f
  23. 9
      src/ruku.f
  24. 24
      src/swavelet.f
  25. 56
      src/transfs2t.f

2
configure.ac

@ -1,4 +1,4 @@
AC_INIT([fomosto-qssp], [2017])
AC_INIT([fomosto-qssp], [2020])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AC_PROG_F77
AC_CONFIG_FILES([

4
doc/examples/qssp2017-test.inp

@ -130,9 +130,9 @@
# [deg] [deg] [sec]
# (Note: Time_reduction = start time of the time window)
#--------------------------------------------------------------------------------------------------------
# disp | velo | acce | strain | strain_rate | stress | stress_rate | rotation | rotation_rate | gravity
# disp | velo | acce | rotation | rot_rate |strain | strain_rate | stress | stress_rate | gravity |gravimeter
#--------------------------------------------------------------------------------------------------------
1 1 1 0 0 0 0 0 0 1
1 1 1 0 0 0 0 0 0 0 0
'test'
8000.0
0 0.05 0.25

8
src/Makefile.am

@ -1,9 +1,11 @@
bin_PROGRAMS = fomosto_qssp2017
fomosto_qssp2017_SOURCES = banddpasss.f brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f disazi.f \
bin_PROGRAMS = fomosto_qssp2020
fomosto_qssp2020_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 \
qpfftinv.f qpgetinp.f qpgrnspec.f legendre.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
clean-local:
-rm -f *.mod

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

6
src/qplegendre.f → src/legendre.f

@ -1,8 +1,8 @@
subroutine qplegendre(ldeg,raddis)
use qpalloc
subroutine legendre(raddis,plm,ldeg,ldegmax)
implicit none
integer*4 ldeg
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

16
src/qpalloc.f

@ -23,14 +23,14 @@ c
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
@ -38,11 +38,11 @@ c
c
complex*16 comi,comi2
c
character*80 dispout(3),veloout(3),acceout(3),
& rotaout(3),rotarateout(3),
& strainout(6),strainrateout(6),
& stressout(6),stressrateout(6),
& gravout(3),grmout
character*800 dispout(3),veloout(3),acceout(3),
& rotaout(3),rotarateout(3),
& strainout(6),strainrateout(6),
& stressout(6),stressrateout(6),
& gravout(3),grmout
c
logical*2,allocatable:: ksmallp(:),ksmalls(:),ksmallt(:)
c
@ -111,7 +111,7 @@ c
& expl(:),clvd(:),ss12(:),
& ss11(:),ds31(:),ds23(:)
c
character*80,allocatable:: specfile(:),uspecfile(:),vspecfile(:),
character*800,allocatable:: specfile(:),uspecfile(:),vspecfile(:),
& wspecfile(:),especfile(:),fspecfile(:),gspecfile(:),
& pspecfile(:),qspecfile(:)
character*10,allocatable:: rname(:)

27
src/qpdifmats.f

@ -7,6 +7,7 @@
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
c
@ -19,12 +20,32 @@ 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
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
if(rr.le.rrlw(ly))then
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
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
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
c
mat(1,1)=(c1-c2*clarr/cxirr)/crr
mat(1,2)=c1/cxirr/crr

256
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,57 @@ 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 butterworth(nbpf,f2corner,df,nf,cswap)
else
call bandpass(nlpf,f1corner,f2corner,df,nf,lpf)
call bandpass(nbpf,f1corner,f2corner,df,nf,cswap)
endif
mf=1
do lf=2*nf,nf+2,-1
mf=mf+1
cswap(lf)=dconjg(cswap(mf))
enddo
cswap(nf+1)=(0.d0,0.d0)
call four1w(cswap,dswap,2*nf,+1)
do lf=1,2*nf
t=dble(lf-1)*dt
cswap(lf)=cswap(lf)*dcmplx(dexp(PI2*fi*t)*df,0.d0)
enddo
call four1w(cswap,dswap,2*nf,-1)
do lf=1,nf
bpf(lf)=cswap(lf)*dcmplx(dt,0.d0)
enddo
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 +92,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

37
src/qpgetinp.f

@ -7,9 +7,9 @@ c work space
c
integer*4 i,j,l,ir,ig,isg,is,is1,flen,iswap,nhypo,sdfsel,ierr
integer*4 ipath,isurf
real*8 twindow,twinout,suppress,munit,sfe,sfn,sfz
real*8 twindow,twinout,suppress,munit,sfe,sfn,sfz,omi
real*8 strike,dip,rake,depdif,dswap(14)
character*80 grndir,outfile,fswap
character*800 grndir,outfile,fswap
c
c uniform receiver depth
c ======================
@ -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
@ -102,7 +103,20 @@ c
endif
if(ldegcut.lt.ldegmin)ldegcut=ldegmin
ldegmin=min0(max0(1+ndmax,ldegmin),ldegcut)
ldegmax=ldegcut+ndmax+1
ldegmax=ldegcut+1+ndmax
c
omi=PI*fcut
if(idint(rearth*omi*slwmax).gt.ldegmax)then
write(*,'(a)')' Warning from qpgrnspec:'
write(*,'(a)')' Maximum harmonic degree not large enough!'
write(*,'(a,i6)')' Suggestion: increase the max. degree to >= ',
& 1+idint(rearth*omi*slwmax)
write(*,'(a,f10.5,a)')' ... or decrease frequency cutoff to <=',
& dble(ldegcut)/(rearth*PI2*slwmax),' Hz'
write(*,'(a,f10.5,a)')' ... or decrease slowness cutoff to <=',
& dble(ldegcut)*KM2M/(rearth*omi),' s/km'
stop
endif
c
allocate(plm(0:ldegmax,0:2),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: plm not allocated!'
@ -260,7 +274,7 @@ c
enddo
enddo
c
do flen=80,1,-1
do flen=800,1,-1
if(grndir(flen:flen).ne.' ')goto 200
enddo
200 continue
@ -453,9 +467,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
@ -481,7 +495,7 @@ c
read(unit,*)latr(ir),lonr(ir),rname(ir),tred(ir)
enddo
c
do flen=80,1,-1
do flen=800,1,-1
if(outfile(flen:flen).ne.' ')goto 300
enddo
300 continue
@ -619,6 +633,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

63
src/qpgrnspec.f

@ -138,16 +138,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 +185,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,9 +196,7 @@ c
slwcut=dmin1(slwcut,fac/vpup(lylwa))
endif
endif
if(idint(rearth*omi*slwcut).gt.ldegcut)then
print *,' Warning from qpgrnspec: ldegcut may be too small!'
endif
c
ldegup=min0(ldegcut,max0(ldeg0,idint(rearth*omi*slwcut)))
c
if(fgr.gt.0.d0.or.ldeggr.gt.0)then
@ -218,6 +228,19 @@ c
enddo
enddo
endif
c
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
enddo
c
do lf=1,nfcut
f=dble(lf-1)*df
@ -249,7 +272,7 @@ c
ksmallt(ly)=.true.
enddo
c
ldegf=min0(ldegup,max0(ldeg0,idint(rearth*omi*slwcut)))
ldegf=min0(ldegup,ldeg0+idint(dble(ldegup-ldeg0)*f/fcut))
setzh12a=.false.
setzh12t=.false.
c
@ -440,12 +463,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 +478,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 +487,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 +500,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
@ -621,7 +644,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)

5
src/qpmain.f

@ -25,7 +25,8 @@ c
print *,'# based on #'
print *,'# a spherically symmetric earth model #'
print *,'# #'
print *,'# (Version 2017) #'
print *,'# (Version 2020) #'
print *,'# Last update (correction of errs): 2020-04-14 #'
print *,'# #'
print *,'# by #'
print *,'# Rongjiang Wang #'
@ -68,7 +69,7 @@ c
runtime=time()-runtime
write(*,'(a)')' #############################################'
write(*,'(a)')' # #'
write(*,'(a)')' # End of computations with atmqssp #'
write(*,'(a)')' # End of computations with qssp2020 #'
write(*,'(a)')' # #'
write(*,'(a,i10,a)')' # Run time: ',runtime,
+ ' sec #'

2
src/qppsvkern.f

@ -11,7 +11,7 @@ c
complex*16 ypsv(6,4)
c
integer i,ly,istp,lyup,lylw,lwup
double precision xmin
real*8 xmin
complex*16 cldeg,cruplw,cxp,cxs
c
complex*16 c1,c2

46
src/qpqmodel.f

@ -10,47 +10,59 @@ c
c
integer*4 ly
real*8 fac,pup,plw,sup,slw
complex*16 cqpup,cqsup,cqplw,cqslw
c
if(.not.dispersion)then
if(.not.dispersion.or.f.ge.FSBREF)then
fac=0.d0
else if(f.le.FSBLW)then
fac=dlog(FSBLW/FSBREF)/PI
else if(f.ge.FSBUP)then
fac=dlog(FSBUP/FSBREF)/PI
else
fac=dlog(f/FSBREF)/PI
endif
c
if(1.d0+fac/qsmin.le.0.d0)then
stop ' Error in qpqmodel: too small Qs value!'
stop 'Error in qpqmodel: too small Qs value!'
endif
do ly=1,ly0
if(f.gt.0.d0)then
cqpup=dcmplx(1.d0,0.5d0/qpup(ly))
cqplw=dcmplx(1.d0,0.5d0/qplw(ly))
else
cqpup=(1.d0,0.d0)
cqplw=(1.d0,0.d0)
endif
pup=1.d0+fac/qpup(ly)
plw=1.d0+fac/qplw(ly)
c
cvpup(ly)=dcmplx(vpup(ly)*pup/cdabs(cqpup),0.d0)*cqpup
cvplw(ly)=dcmplx(vplw(ly)*plw/cdabs(cqplw),0.d0)*cqplw
cvp(ly)=(0.5d0,0.d0)*(cvpup(ly)+cvplw(ly))
c
if(ly.ge.lyob.and.ly.lt.lycm.or.ly.ge.lycc)then
sup=1.d0+fac/qsup(ly)
slw=1.d0+fac/qslw(ly)
cvsup(ly)=dcmplx(vsup(ly),0.d0)
& *dcmplx(sup,0.5d0/qsup(ly))
cvslw(ly)=dcmplx(vslw(ly),0.d0)
& *dcmplx(slw,0.5d0/qslw(ly))
if(f.gt.0.d0)then
cqsup=dcmplx(1.d0,0.5d0/qsup(ly))
cqslw=dcmplx(1.d0,0.5d0/qslw(ly))
else
cqsup=(1.d0,0.d0)
cqslw=(1.d0,0.d0)
endif
cvsup(ly)=dcmplx(vsup(ly)*sup/cdabs(cqsup),0.d0)*cqsup
cvslw(ly)=dcmplx(vslw(ly)*slw/cdabs(cqslw),0.d0)*cqslw
cvs(ly)=(0.5d0,0.d0)*(cvsup(ly)+cvslw(ly))
else
cvsup(ly)=(0.d0,0.d0)
cvslw(ly)=(0.d0,0.d0)
cvs(ly)=(0.d0,0.d0)
endif
pup=1.d0+fac/qpup(ly)
plw=1.d0+fac/qplw(ly)
cvpup(ly)=dcmplx(vpup(ly),0.d0)
& *dcmplx(pup,0.5d0/qpup(ly))
cvplw(ly)=dcmplx(vplw(ly),0.d0)
& *dcmplx(plw,0.5d0/qplw(ly))
cvp(ly)=(0.5d0,0.d0)*(cvpup(ly)+cvplw(ly))
c
cmuup(ly)=croup(ly)*cvsup(ly)**2
claup(ly)=croup(ly)*cvpup(ly)**2-(2.d0,0.d0)*cmuup(ly)
cmulw(ly)=crolw(ly)*cvslw(ly)**2
clalw(ly)=crolw(ly)*cvplw(ly)**2-(2.d0,0.d0)*cmulw(ly)
cmu(ly)=cro(ly)*cvs(ly)**2
cla(ly)=cro(ly)*cvp(ly)**2-(2.d0,0.d0)*cmu(ly)
cmu(ly)=(0.5d0,0.d0)*(cmuup(ly)+cmulw(ly))
cla(ly)=(0.5d0,0.d0)*(claup(ly)+clalw(ly))
enddo
return
end

187
src/qpsprop.f

@ -5,23 +5,22 @@ c
c calculation of spheroidal response
c ypsv(6,4): solution vector (complex)
c
integer*4 lyup,lylw
integer*4 ldeg,lyup,lylw
complex*16 ypsv(6,4)
c
c work space
c
integer*4 i,j,j0,istp,ly,ldeg,key
integer*4 i,j,j0,istp,ly,lyswap,key
real*8 y4max
complex*16 cdet,cyswap,uc
complex*16 c0(6,3),c1(6,3),cc0(4,2),cc1(4,2)
complex*16 y0(6,3),y1(6,3),yc(6,3)
complex*16 y0(6,3),y1(6,3)
complex*16 yup(6,3),ylw(6,3),yupc(4,2),ylwc(4,2)
complex*16 wave(6),orth(3,3),orthc(2,2)
complex*16 coef6(6,6),b6(6,4),coef4(4,4),b4(4,2)
logical*2 comsys
c
complex*16 c3
data c3/(3.d0,0.d0)/
complex*16 c2,c3,c6
data c2,c3,c6/(2.d0,0.d0),(3.d0,0.d0),(6.d0,0.d0)/
c
c===============================================================================
c
@ -29,8 +28,6 @@ c propagation from surface to atmosphere/ocean bottom
c
if(lylwa.gt.lylw.or.cdabs(comi).le.0.d0.and.
& (lyr.lt.lyob.or.lyr.ge.lycm.and.lyr.lt.lycc))return
c
comsys=ldeg.eq.1.and.freesurf.and.lyup.eq.1
c
if(lyob.gt.lyup)then
do j=1,2
@ -41,6 +38,10 @@ c
c
if(freesurf.and.lyup.eq.1)then
yupc(1,1)=(1.d0,0.d0)
if(ldeg.eq.1)then
yupc(1,2)=-(1.d0,0.d0)/cgrup(1)
yupc(2,2)=-(1.d0,0.d0)/cgrup(1)
endif
yupc(3,2)=(1.d0,0.d0)
else
call qpstart4a(ldeg,lyup,2,yupc)
@ -58,19 +59,6 @@ c
y0(i,3)=(0.d0,0.d0)
enddo
endif
c
if(comsys)then
do j=1,2
yc(1,j)=yupc(1,j)
yc(2,j)=yupc(2,j)
yc(3,j)=-yupc(2,j)/(cro(lyup)*comi2*crrup(lyup)**2)
yc(5,j)=yupc(3,j)
yc(6,j)=yupc(4,j)
enddo
do i=1,6
yc(i,3)=(0.d0,0.d0)
enddo
endif
endif
c
do ly=lyup,min0(lyob,lys)-1
@ -104,18 +92,6 @@ c
enddo
enddo
endif
if(comsys)then
c
c orthonormalization of the receiver vectors
c
call caxcb(yc,orthc,6,2,2,y1)
call cmemcpy(y1,yc,12)
do j=1,2
do i=1,6
yc(i,j)=yc(i,j)*wave(2*j)
enddo
enddo
endif
c
cc1(1,1)=cc1(1,1)*wave(1)*wave(2)
cc1(2,1)=(1.d0,0.d0)
@ -172,12 +148,15 @@ c
enddo
yup(1,1)=(1.d0,0.d0)
yup(3,2)=(1.d0,0.d0)
if(ldeg.eq.1)then
yup(1,3)=-(1.d0,0.d0)/cgrup(1)
yup(3,3)=-(1.d0,0.d0)/cgrup(1)
endif
yup(5,3)=(1.d0,0.d0)
else
call qpstart6a(ldeg,lyup,2,yup)
endif
if(lyr.eq.lyup)call cmemcpy(yup,y0,18)
if(comsys)call cmemcpy(yup,yc,18)
endif
endif
c
@ -226,18 +205,6 @@ c
enddo
enddo
endif
if(comsys)then
c
c orthonormalization of the receiver vectors
c
call caxcb(yc,orth,6,3,3,y1)
call cmemcpy(y1,yc,18)
do j=1,3
do i=1,6
yc(i,j)=yc(i,j)*wave(2*j)
enddo
enddo
endif
c
c1(1,1)=c1(1,1)*wave(1)*wave(2)
c1(2,1)=(1.d0,0.d0)
@ -303,21 +270,6 @@ c
y0(i,3)=(0.d0,0.d0)
enddo
endif
if(comsys)then
do i=1,6
cyswap=yc(i,j0)
yc(i,j0)=yc(i,3)
yc(i,3)=cyswap
enddo
do j=1,2
do i=1,6
yc(i,j)=yc(i,j)-yup(4,j)*yc(i,3)/yup(4,3)
enddo
enddo
do i=1,6
yc(i,3)=(0.d0,0.d0)
enddo
endif
else if(lyup.lt.lycc)then
call qpstart4a(ldeg,lyup,2,yupc)
if(lyr.eq.lyup)then
@ -333,19 +285,6 @@ c
y0(i,3)=(0.d0,0.d0)
enddo
endif
if(comsys)then
do j=1,2
yc(1,j)=yupc(1,j)
yc(2,j)=yupc(2,j)
yc(3,j)=-yupc(2,j)/(cro(lyup)*comi2*crrup(lyup)**2)
yc(4,j)=(0.d0,0.d0)
yc(5,j)=yupc(3,j)
yc(6,j)=yupc(4,j)
enddo
do i=1,6
yc(i,3)=(0.d0,0.d0)
enddo
endif
endif
endif
c
@ -380,18 +319,6 @@ c
enddo
enddo
endif
if(comsys)then
c
c orthonormalization of the receiver vectors
c
call caxcb(yc,orthc,6,2,2,y1)
call cmemcpy(y1,yc,12)
do j=1,2
do i=1,6
yc(i,j)=yc(i,j)*wave(2*j)
enddo
enddo
endif
c
cc1(1,1)=cc1(1,1)*wave(1)*wave(2)
cc1(2,1)=(1.d0,0.d0)
@ -445,7 +372,6 @@ c
else if(lyup.lt.ly0)then
call qpstart6a(ldeg,lyup,2,yup)
if(lyr.eq.lyup)call cmemcpy(yup,y0,18)
if(comsys)call cmemcpy(yup,yc,18)
endif
endif
c
@ -490,18 +416,6 @@ c
enddo
enddo
endif
if(comsys)then
c
c orthonormalization of the receiver vectors
c
call caxcb(yc,orth,6,3,3,y1)
call cmemcpy(y1,yc,18)
do j=1,3
do i=1,6
yc(i,j)=yc(i,j)*wave(2*j)
enddo
enddo
endif
c
c1(1,1)=c1(1,1)*wave(1)*wave(2)
c1(2,1)=(1.d0,0.d0)
@ -551,17 +465,6 @@ c
c yup(5,j)=yup(5,j)
yup(6,j)=yup(6,j)/crrup(lys)
enddo
c
if(comsys)then
do j=1,3
yc(1,j)=yc(1,j)/crrup(lyup)
yc(2,j)=yc(2,j)/crrup(lyup)**2
yc(3,j)=yc(3,j)/crrup(lyup)
yc(4,j)=yc(4,j)/crrup(lyup)**2
c yc(5,j)=yc(5,j)
yc(6,j)=yc(6,j)/crrup(lyup)
enddo
endif
c
c===============================================================================
c
@ -1106,19 +1009,6 @@ c
ypsv(i,istp)=(0.d0,0.d0)
enddo
enddo
if(comsys)then
do istp=1,2
uc=(0.d0,0.d0)
do j=1,2
uc=uc+b4(j,istp)*yc(5,j)
enddo
ypsv(5,istp)=ypsv(5,istp)-uc
ypsv(6,istp)=ypsv(6,istp)-uc*c3/crrup(lyup)
uc=uc/cgrup(lyup)
ypsv(1,istp)=ypsv(1,istp)-uc
ypsv(3,istp)=ypsv(3,istp)-uc
enddo
endif
else
do istp=1,4
do i=1,6
@ -1157,19 +1047,6 @@ c
enddo
enddo
endif
if(comsys)then
do istp=1,4
uc=(0.d0,0.d0)
do j=1,3
uc=uc+b6(j,istp)*yc(5,j)
enddo
ypsv(5,istp)=ypsv(5,istp)-uc
ypsv(6,istp)=ypsv(6,istp)-uc*c3/crrup(lyup)
uc=uc/cgrup(lyup)
ypsv(1,istp)=ypsv(1,istp)-uc
ypsv(3,istp)=ypsv(3,istp)-uc
enddo
endif
endif
c
if(lylwa.le.0)return
@ -1182,7 +1059,8 @@ c
c
c lowest layer is within inner core
c
call qpstart6a(ldeg,lylwa,1,ylw)
lyswap=lylwa
call qpstart6a(ldeg,lyswap,1,ylw)
if(lylwa.eq.lyr.and.lylwa.gt.lys)call cmemcpy(ylw,y0,18)
endif
c
@ -1299,7 +1177,8 @@ c
c
c lowest layer is within the liquid core
c
call qpstart4a(ldeg,lylwa,1,ylwc)
lyswap=lylwa
call qpstart4a(ldeg,lyswap,1,ylwc)
c
if(lylwa.eq.lyr.and.lylwa.gt.lys)then
do j=1,2
@ -1406,7 +1285,8 @@ c
c
c lowest layer is within the mantle
c
call qpstart6a(ldeg,lylwa,1,ylw)
lyswap=lylwa
call qpstart6a(ldeg,lyswap,1,ylw)
if(lylwa.eq.lyr.and.lylwa.gt.lys)then
call cmemcpy(ylw,y0,18)
endif
@ -1526,7 +1406,8 @@ c
c
c lowest layer is within the atmosphere
c
call qpstart4a(ldeg,lylwa,1,ylwc)
lyswap=lylwa
call qpstart4a(ldeg,lyswap,1,ylwc)
c
if(lylwa.eq.lyr.and.lylwa.gt.lys)then