Commit 71568945 authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

upgrade to qseis6a

parent 845d2900
......@@ -9,4 +9,4 @@ Makefile
config.log
config.status
*.o
src/qseis
src/fomosto_qseis2006a
AC_INIT([qseis], [1.2])
AC_INIT([fomosto-qseis], [2006a])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AC_PROG_F77
AC_CONFIG_FILES([
......
This diff is collapsed.
bin_PROGRAMS = qseis
qseis_SOURCES = bessj0.f bessj1.f bessj.f caxcb.f cdgemp.f cmemcpy.f four1.f qsam2ve.f qsbsj.f qsfftinv.f qsgetinp.f qshksh.f qskern.f qslayer.f qsmain.f qsmultis.f qspsv.f qsqmodel.f qssh.f qssource.f qssublay.f qsve2am.f qswavelet.f qswaveno.f qswvint.f taper.f wavelet.f skip_comments.f qsglobal.h
bin_PROGRAMS = fomosto_qseis2006a
fomosto_qseis2006a_SOURCES = qsmatrix0.f bessj0.f bessj1.f bessj.f caxcb.f cdgemp.f cmemcpy.f four1.f qsam2ve.f qsbsj.f qsfftinv.f qsgetinp.f qshksh.f qskern.f qslayer.f qsmain.f qsmultis.f qspsv.f qsqmodel.f qssh.f qssource.f qssublay.f qsve2am.f qswavelet.f qswaveno.f qswvint.f taper.f wavelet.f skip_comments.f qsglobal.h
subroutine qsam2ve(n,ck,y,c,nj,rsite)
subroutine qsam2ve(n,ck,z,dynamic,y,c,nj,rsite)
implicit none
c
c converting wave amplitudes to displacement-stress vectors
......@@ -6,14 +6,16 @@ c by Langer block-diagonal decomposition technique for com-
c putational efficiency
c
integer n,nj
double precision z
double complex ck
double complex y(4,nj),c(4,nj)
logical rsite
logical dynamic,rsite
c
include 'qsglobal.h'
c
integer j
double complex bp,bs,swap(4)
integer i,j,m
double complex bp,bs
double complex swap(4),a(4,4)
c
if(rsite)then
bp=wbrs(n)*kprs(n)
......@@ -29,7 +31,7 @@ c
y(3,j)=ck*swap(2)-ksrs(n)*swap(3)
y(4,j)=bp*swap(1)+wars(n)*swap(4)
enddo
else
else if(dynamic)then
bp=wb(n)*kp(n)
bs=wb(n)*ks(n)
do j=1,nj
......@@ -43,6 +45,16 @@ c
y(3,j)=ck*swap(2)-ks(n)*swap(3)
y(4,j)=bp*swap(1)+wa(n)*swap(4)
enddo
else
call qsmatrix0(a,ck,z,n)
do j=1,nj
do i=1,4
y(i,j)=(0.d0,0.d0)
do m=1,4
y(i,j)=y(i,j)+a(i,m)*c(m,j)
enddo
enddo
enddo
endif
c
return
......
......@@ -48,10 +48,10 @@ c
c
c seismometer filtering
c
if(asm.ne.1.d0)then
if(cdabs(asm).ne.1.d0)then
do lf=1,nf
do ir=1,nr
cy(lf,ir)=dcmplx(asm,0.d0)*cy(lf,ir)
cy(lf,ir)=asm*cy(lf,ir)
enddo
enddo
endif
......
......@@ -560,7 +560,7 @@ c
vps=vps*rr/rr0
vss=vss*rr/rr0
endif
call qssource(ros,vps,vss)
call qssource(ros,dcmplx(vps,0.d0),dcmplx(vss,0.d0))
c
if(v0.gt.0.d0)then
v00=1.d0/v0
......
......@@ -6,16 +6,16 @@ c nrmax: max. number of traces;
c nfmax: max. number of frequency samples
c
integer nzmax,lmax,nrmax,nfmax,ndtransmax
parameter(lmax=1500)
parameter(lmax=2500)
parameter(nzmax=lmax+2)
parameter(nrmax=101,nfmax=2048)
parameter(nrmax=101,nfmax=8192)
parameter(ndtransmax=4)
c
c INDEX PARAMETERS FOR BESSEL FUNCTION TABLES
c ===========================================
c
integer nbsjmax
parameter(nbsjmax=30000)
parameter(nbsjmax=150000)
c
c INDEX PARAMETERS FOR SEISMOMETER CHARACTERISTICS
c ================================================
......@@ -132,7 +132,7 @@ c source parameters
c
integer ls,ms(6),ics(6)
double precision zs
double precision sfct0(6,6),sfct1(6,6)
double complex sfct0(6,6),sfct1(6,6)
common /source/ zs,sfct0,sfct1,ls,ms,ics
c
c path filtering
......@@ -159,7 +159,7 @@ c
c seismometer filtering
c
integer nroot,npole
double precision asm
double complex asm
double complex root(nrootmax),pole(npolemax)
common /seismometer/ root,pole,asm,nroot,npole
c
......
......@@ -112,7 +112,7 @@ c
c determine layer no of each depth
c
li=1
zswap=h(1)
zswap=hrs(1)
nnors(1)=1
do l=2,lprs
if(z(l).ge.zswap.and.li.lt.n0rs)then
......
......@@ -23,7 +23,7 @@ c
print *,'# Q QQ S E I S #'
print *,'# QQQQ SSSS EEEEE III SSSS #'
print *,'# #'
print *,'# (Version 2006) #'
print *,'# (Version 2006a) #'
print *,'# #'
print *,'# #'
print *,'# by #'
......@@ -64,7 +64,7 @@ c
runtime=time()-runtime
write(*,'(a)')' #############################################'
write(*,'(a)')' # #'
write(*,'(a)')' # End of computations with qseis06 #'
write(*,'(a)')' # End of computations with qseis06a #'
write(*,'(a)')' # #'
write(*,'(a,i10,a)')' # Run time: ',runtime,
+ ' sec #'
......
subroutine qsmatrix0(a,ck,z,n)
implicit none
c
integer n
double precision z
double complex ck
double complex a(4,4)
c
include 'qsglobal.h'
c
double complex c0,c1,c2
data c0,c1,c2/(0.d0,0.d0),(1.d0,0.d0),(2.d0,0.d0)/
c
integer i
double complex ck2,cz,cxi,cfac,b(4,4)
double complex cxps,cdkp,cdks,cps,csp,cdeps,cdesp
c
ck2=ck*ck
c
a(1,1)=kp(n)
a(2,1)=c2*cmu(n)*ck2-acc(n)
a(3,1)=ck
a(4,1)=c2*ck*cmu(n)*kp(n)
a(1,2)=-kp(n)
a(2,2)=a(2,1)
a(3,2)=ck
a(4,2)=-a(4,1)
c
b(1,3)=ck
b(2,3)=c2*ck*cmu(n)*ks(n)
b(3,3)=ks(n)
b(4,3)=c2*cmu(n)*ck2-acc(n)
b(1,4)=ck
b(2,4)=-b(2,3)
b(3,4)=-ks(n)
b(4,4)=b(4,3)
c
cz=dcmplx(z,0.d0)
cxi=cla(n)+c2*cmu(n)
c
c cxps=(kp-ks)*z
c cdkp=(k-kp)*k^2*vp^2/omiga^2
c cdks=(k-ks)*k^2*vp^2/omiga^2
c cdeps=(1-exp(cxps))*k^2*vp^2/omiga^2
c cdesp=(1-exp(-cxps))*k^2*vp^2/omiga^2
c
cxps=acc(n)*(cla(n)+cmu(n))*cz/(cmu(n)*cxi*(kp(n)+ks(n)))
cdkp=ck2/(ck+kp(n))
cdks=ck2*cxi/((ck+ks(n))*cmu(n))
if(dabs(z).gt.0.d0)then
i=1
cps=(1.d0,0.d0)
csp=(1.d0,0.d0)
cdeps=(1.d0,0.d0)
cdesp=(1.d0,0.d0)
10 i=i+1
cfac=cxps/dcmplx(dble(i),0.d0)
cps=cps*cfac
csp=-csp*cfac
cdeps=cdeps+cps
cdesp=cdesp+csp
if(cdabs(cps).gt.1.0d-12)goto 10
cfac=(cla(n)+cmu(n))*ck2*cz/(cmu(n)*(kp(n)+ks(n)))
cdeps=-cdeps*cfac
cdesp=cdesp*cfac
else
cdeps=0.d0
cdesp=0.d0
endif
c
c define stable propagator for omiga -> 0:
c
c a(i,3)=[a(i,3)-a(i,1)*exp(cxps)]*k^2*vp^2/omiga^2
c a(i,4)=[a(i,4)+a(i,2)*exp(-cxps)]*k^2*vp^2/omiga^2
c
a(1,3)=cdkp+kp(n)*cdeps
a(2,3)=c2*cmu(n)*ck*(ck*cdeps-cdks)+cxi*ck2*cdexp(cxps)
a(3,3)=-cdks+ck*cdeps
a(4,3)=c2*cmu(n)*ck*(cdkp+kp(n)*cdeps)-cxi*ck2
a(1,4)=cdkp+kp(n)*cdesp
a(2,4)=c2*cmu(n)*ck*(cdks-ck*cdesp)-cxi*ck2*cdexp(-cxps)
a(3,4)=cdks-ck*cdesp
a(4,4)=c2*cmu(n)*ck*(cdkp+kp(n)*cdesp)-cxi*ck2
c
return
end
\ No newline at end of file
......@@ -51,13 +51,13 @@ c
open(unit,file=outfile(icmp,istp),status='old')
read(unit,'(a1)')textline
do lf=1,nf
read(unit,*)t(2*lf-1),(y1(ir),ir=1,nr)
read(unit,*)t(2*lf),(y2(ir),ir=1,nr)
read(unit,*,end=100)t(2*lf-1),(y1(ir),ir=1,nr)
read(unit,*,end=100)t(2*lf),(y2(ir),ir=1,nr)
do ir=1,nr
grns(lf,icmp,ir,istp)=dcmplx(y1(ir),y2(ir))
enddo
enddo
close(unit)
100 close(unit)
endif
enddo
enddo
......
......@@ -20,6 +20,7 @@ c
double complex y1(4,2),yup(4,2),ylw(4,2),orth(2,2)
double complex coef(4,4),cnorm(2),coefrs(2,2),brs(2,6)
logical rsite
logical dynamic(nzmax)
c
double complex c2
data c2/(2.d0,0.d0)/
......@@ -34,6 +35,13 @@ c determination of starting upper sublayer
c
rsite=.false.
c
do l=1,lp
n=nno(l)
dynamic(l)=isurf.ne.0.or.n0rs.gt.0.or.
& cdabs(kp(n)-ks(n)).gt.0.1d0*k.or.
& cdabs(kp(n)-ks(n))*hp(l).gt.0.1d0.or.
& .not.(pup(l).and.pdw(l).and.svup(l).and.svdw(l))
enddo
if(lup.eq.1.and.isurf.eq.0)then
do j=1,2
do i=1,4
......@@ -64,7 +72,7 @@ c
c
c determination of propagation matrix
c
call qsve2am(n,ck,yup,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l-1),yup,c0,2,rsite)
pwave=cdexp(-kp(n)*ch0)
swave=cdexp(-ks(n)*ch0)
c
......@@ -113,20 +121,20 @@ c
c1(2,1)=(0.d0,0.d0)
c1(2,2)=(0.d0,0.d0)
if(l-1.eq.lzr)then
call qsve2am(n,ck,y0,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l-1),y0,c0,2,rsite)
c0(2,1)=(0.d0,0.d0)
c0(2,2)=(0.d0,0.d0)
call qsam2ve(n,ck,y0,c0,2,rsite)
call qsam2ve(n,ck,0.d0,dynamic(l-1),y0,c0,2,rsite)
endif
endif
if(.not.svdw(l-1))then
c1(4,1)=(0.d0,0.d0)
c1(4,2)=(0.d0,0.d0)
if(l-1.eq.lzr)then
call qsve2am(n,ck,y0,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l-1),y0,c0,2,rsite)
c0(4,1)=(0.d0,0.d0)
c0(4,2)=(0.d0,0.d0)
call qsam2ve(n,ck,y0,c0,2,rsite)
call qsam2ve(n,ck,0.d0,dynamic(l-1),y0,c0,2,rsite)
endif
endif
c
......@@ -142,7 +150,7 @@ c
c c1(3,2)=c1(3,2)
c1(4,2)=c1(4,2)*swave*swave
c
call qsam2ve(n,ck,yup,c1,2,rsite)
call qsam2ve(n,ck,hp(l-1),dynamic(l-1),yup,c1,2,rsite)
c
if(l.eq.lzr)call cmemcpy(yup,y0,8)
enddo
......@@ -170,7 +178,7 @@ c
c
c determination of propagation matrix
c
call qsve2am(n,ck,ylw,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l),ylw,c0,2,rsite)
pwave=cdexp(-kp(n)*ch0)
swave=cdexp(-ks(n)*ch0)
c
......@@ -201,20 +209,20 @@ c
c1(1,1)=(0.d0,0.d0)
c1(1,2)=(0.d0,0.d0)
if(l+1.eq.lzr)then
call qsve2am(n,ck,y0,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
c0(1,1)=(0.d0,0.d0)
c0(1,2)=(0.d0,0.d0)
call qsam2ve(n,ck,y0,c0,2,rsite)
call qsam2ve(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
endif
endif
if(.not.svup(l))then
c1(3,1)=(0.d0,0.d0)
c1(3,2)=(0.d0,0.d0)
if(l+1.eq.lzr)then
call qsve2am(n,ck,y0,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
c0(3,1)=(0.d0,0.d0)
c0(3,2)=(0.d0,0.d0)
call qsam2ve(n,ck,y0,c0,2,rsite)
call qsam2ve(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
endif
endif
if(.not.pdw(l))then
......@@ -248,7 +256,7 @@ c
c1(3,2)=c1(3,2)*swave*swave
c c1(4,2)=c1(4,2)
c
call qsam2ve(n,ck,ylw,c1,2,rsite)
call qsam2ve(n,ck,-hp(l),dynamic(l),ylw,c1,2,rsite)
if(l.gt.ls.and.l.eq.lzr)call cmemcpy(ylw,y0,8)
enddo
c
......@@ -258,7 +266,7 @@ c===============================================================================
c
do istp=1,6
do i=1,4
b(i,istp)=dcmplx(sfct0(i,istp)+k*sfct1(i,istp),0.d0)
b(i,istp)=sfct0(i,istp)+ck*sfct1(i,istp)
enddo
enddo
do i=1,4
......@@ -313,7 +321,7 @@ c
c
c determination of propagation matrix
c
call qsve2am(n,ck,ylw,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l),ylw,c0,2,rsite)
pwave=cdexp(-kp(n)*ch0)
swave=cdexp(-ks(n)*ch0)
c
......@@ -344,20 +352,20 @@ c
c1(1,1)=(0.d0,0.d0)
c1(1,2)=(0.d0,0.d0)
if(l+1.eq.lzr)then
call qsve2am(n,ck,y0,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
c0(1,1)=(0.d0,0.d0)
c0(1,2)=(0.d0,0.d0)
call qsam2ve(n,ck,y0,c0,2,rsite)
call qsam2ve(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
endif
endif
if(.not.svup(l))then
c1(3,1)=(0.d0,0.d0)
c1(3,2)=(0.d0,0.d0)
if(l+1.eq.lzr)then
call qsve2am(n,ck,y0,c0,2,rsite)
call qsve2am(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
c0(3,1)=(0.d0,0.d0)
c0(3,2)=(0.d0,0.d0)
call qsam2ve(n,ck,y0,c0,2,rsite)
call qsam2ve(n,ck,0.d0,dynamic(l),y0,c0,2,rsite)
endif
endif
if(.not.pdw(l))then
......@@ -391,7 +399,7 @@ c
c1(3,2)=c1(3,2)*swave*swave
c c1(4,2)=c1(4,2)
c
call qsam2ve(n,ck,ylw,c1,2,rsite)
call qsam2ve(n,ck,-hp(l),dynamic(l),ylw,c1,2,rsite)
if(l.gt.ls.and.l.eq.lzr)call cmemcpy(ylw,y0,8)
enddo
do i=1,4
......@@ -431,7 +439,7 @@ c for receiver-site structure different from source-site structure
c
n=nno(lzr)
do istp=1,6
call qsve2am(n,ck,y(1,istp),cinc(1,istp),1,rsite)
call qsve2am(n,ck,0.d0,.true.,y(1,istp),cinc(1,istp),1,rsite)
enddo
c
c determination of starting upper sublayer
......@@ -469,7 +477,7 @@ c
c
c determination of propagation matrix
c
call qsve2am(n,ck,yup,c0,2,rsite)
call qsve2am(n,ck,0.d0,.true.,yup,c0,2,rsite)
pwave=cdexp(-kprs(n)*ch0)
swave=cdexp(-ksrs(n)*ch0)
c
......@@ -503,11 +511,11 @@ c
c c1(3,2)=c1(3,2)
c1(4,2)=c1(4,2)*swave*swave
c
call qsam2ve(n,ck,yup,c1,2,rsite)
call qsam2ve(n,ck,hprs(l-1),.true.,yup,c1,2,rsite)
if(l.eq.lzrrs)call cmemcpy(yup,y0,8)
enddo
n=nnors(lprs)
call qsve2am(n,ck,yup,c0,2,rsite)
call qsve2am(n,ck,0.d0,.true.,yup,c0,2,rsite)
do istp=1,6
brs(1,istp)=cinc(1,istp)
brs(2,istp)=cinc(3,istp)
......@@ -579,7 +587,7 @@ c
endif
c
do istp=1,6
call qsve2am(1,ck,y(1,istp),b0(1,istp),1,rsite)
call qsve2am(1,ck,0.d0,.true.,y(1,istp),b0(1,istp),1,rsite)
b(1,istp)=-yup(2,1)*b0(1,istp)-yup(2,2)*b0(3,istp)
b(2,istp)=-yup(4,1)*b0(1,istp)-yup(4,2)*b0(3,istp)
......@@ -595,7 +603,7 @@ c
endif
b0(2,istp)=b(1,istp)
b0(4,istp)=b(2,istp)
call qsam2ve(1,ck,y(1,istp),b0(1,istp),1,rsite)
call qsam2ve(1,ck,0.d0,.true.,y(1,istp),b0(1,istp),1,rsite)
enddo
endif
return
......
......@@ -101,7 +101,7 @@ c===============================================================================
c
do istp=1,6
do i=1,2
b(i,istp)=dcmplx(sfct0(i+4,istp)+k*sfct1(i+4,istp),0.d0)
b(i,istp)=sfct0(i+4,istp)+dcmplx(k,0.d0)*sfct1(i+4,istp)
enddo
enddo
do i=1,2
......
subroutine qssource(ros,vps,vss)
subroutine qssource(ros,cvps,cvss)
implicit none
c
double precision ros,vps,vss
double precision ros
double complex cvps,cvss
c
include 'qsglobal.h'
c
......@@ -23,15 +24,15 @@ c explosion source (m11=m22=m33=1)
c
ms(1)=0
ics(1)=1
sfct0(1,1)=-1.d0/(pi2*ros*vps*vps)
sfct1(4,1)=-(vss/vps)**2/pi
sfct0(1,1)=-dcmplx(1.d0/(pi2*ros),0.d0)/(cvps*cvps)
sfct1(4,1)=-(cvss/cvps)**2/dcmplx(pi,0.d0)
c
c istype = 2
c strike-slip (m12=m21=1)
c
ms(2)=2
ics(2)=-1
sfct1(4,2)=1.d0/pi2
sfct1(4,2)=dcmplx(1.d0/pi2,0.d0)
sfct1(6,2)=-sfct1(4,2)
c
c istype = 3
......@@ -39,7 +40,7 @@ c dip-slip (m13=m31=1)
c
ms(3)=1
ics(3)=1
sfct0(3,3)=-1.d0/(pi2*ros*vss*vss)
sfct0(3,3)=-dcmplx(1.d0/(pi2*ros),0.d0)/(cvss*cvss)
sfct0(5,3)=sfct0(3,3)
c
c istp = 4
......@@ -47,22 +48,23 @@ c compensated linear vector dipole (CLVD) (m11=m22=-1/2, M33=1)
c
ms(4)=0
ics(4)=1
sfct0(1,4)=-1.d0/(pi2*ros*vps*vps)
sfct1(4,4)=(3.d0-4.d0*(vss/vps)**2)/(2.d0*pi2)
sfct0(1,4)=-dcmplx(1.d0/(pi2*ros),0.d0)/(cvps*cvps)
sfct1(4,4)=((3.d0,0.d0)-(4.d0,0.d0)*(cvss/cvps)**2)
& /dcmplx(2.d0*pi2,0.d0)
c
c istp = 5
c vertical-single-force (fz=1)
c
ms(5)=0
ics(5)=1
sfct0(2,5)=1.d0/pi2
sfct0(2,5)=dcmplx(1.d0/pi2,0.d0)
c
c istp = 6
c horizontal-single-force (fx=1)
c
ms(6)=1
ics(6)=1
sfct0(4,6)=1.d0/pi2
sfct0(4,6)=dcmplx(1.d0/pi2,0.d0)
sfct0(6,6)=sfct0(4,6)
c
return
......
subroutine qsve2am(n,ck,y,c,nj,rsite)
subroutine qsve2am(n,ck,z,dynamic,y,c,nj,rsite)
implicit none
c
c converting displacement-stress vectors to wave amplitudes
......@@ -6,14 +6,16 @@ c by Langer block-diagonal decomposition technique for com-
c putational efficiency
c
integer n,nj
double precision z
double complex ck
double complex y(4,nj),c(4,nj)
logical rsite
logical dynamic,rsite
c
include 'qsglobal.h'
c
integer j
double complex ca,cadp,cads,swap(4)
integer i,j,m,key
double complex ca,cadp,cads
double complex swap(4),a(4,4),b(4,4)
c
if(rsite)then
ca=(0.5d0,0.d0)/accrs(n)
......@@ -30,7 +32,7 @@ c
c(3,j)=-swap(3)+swap(4)
c(4,j)= swap(3)+swap(4)
enddo
else
else if(dynamic)then
ca=(0.5d0,0.d0)/acc(n)
cadp=ca/kp(n)
cads=ca/ks(n)
......@@ -45,6 +47,28 @@ c
c(3,j)=-swap(3)+swap(4)
c(4,j)= swap(3)+swap(4)
enddo
else
call qsmatrix0(b,ck,z,n)