Commit 845d2900 authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

dos2unix line endings

parent cf6cec8b
subroutine qsam2ve(n,ck,y,c,nj,rsite)
implicit none
c
c converting wave amplitudes to displacement-stress vectors
c by Langer block-diagonal decomposition technique for com-
c putational efficiency
c
c converting wave amplitudes to displacement-stress vectors
c by Langer block-diagonal decomposition technique for com-
c putational efficiency
c
integer n,nj
double complex ck
......@@ -11,37 +11,37 @@ c
logical rsite
c
include 'qsglobal.h'
c
integer j
c
integer j
double complex bp,bs,swap(4)
c
if(rsite)then
bp=wbrs(n)*kprs(n)
bs=wbrs(n)*ksrs(n)
do j=1,nj
swap(1)= c(1,j)-c(2,j)
swap(2)= c(1,j)+c(2,j)
swap(3)=-c(3,j)+c(4,j)
swap(4)= c(3,j)+c(4,j)
c
y(1,j)=kprs(n)*swap(1)+ck*swap(4)
y(2,j)=wars(n)*swap(2)-bs*swap(3)
y(3,j)=ck*swap(2)-ksrs(n)*swap(3)
y(4,j)=bp*swap(1)+wars(n)*swap(4)
if(rsite)then
bp=wbrs(n)*kprs(n)
bs=wbrs(n)*ksrs(n)
do j=1,nj
swap(1)= c(1,j)-c(2,j)
swap(2)= c(1,j)+c(2,j)
swap(3)=-c(3,j)+c(4,j)
swap(4)= c(3,j)+c(4,j)
c
y(1,j)=kprs(n)*swap(1)+ck*swap(4)
y(2,j)=wars(n)*swap(2)-bs*swap(3)
y(3,j)=ck*swap(2)-ksrs(n)*swap(3)
y(4,j)=bp*swap(1)+wars(n)*swap(4)
enddo
else
bp=wb(n)*kp(n)
bs=wb(n)*ks(n)
do j=1,nj
swap(1)= c(1,j)-c(2,j)
swap(2)= c(1,j)+c(2,j)
swap(3)=-c(3,j)+c(4,j)
swap(4)= c(3,j)+c(4,j)
c
y(1,j)=kp(n)*swap(1)+ck*swap(4)
y(2,j)=wa(n)*swap(2)-bs*swap(3)
y(3,j)=ck*swap(2)-ks(n)*swap(3)
y(4,j)=bp*swap(1)+wa(n)*swap(4)
else
bp=wb(n)*kp(n)
bs=wb(n)*ks(n)
do j=1,nj
swap(1)= c(1,j)-c(2,j)
swap(2)= c(1,j)+c(2,j)
swap(3)=-c(3,j)+c(4,j)
swap(4)= c(3,j)+c(4,j)
c
y(1,j)=kp(n)*swap(1)+ck*swap(4)
y(2,j)=wa(n)*swap(2)-bs*swap(3)
y(3,j)=ck*swap(2)-ks(n)*swap(3)
y(4,j)=bp*swap(1)+wa(n)*swap(4)
enddo
endif
c
......
......@@ -25,9 +25,9 @@ c
else
bsj(ik,3,ir)=bessj(3,x)
endif
bsj(ik,-1,ir)=-bsj(ik,1,ir)
do i=-1,3
bsj(ik,i,ir)=bsj(ik,i,ir)*geospr(ir)
bsj(ik,-1,ir)=-bsj(ik,1,ir)
do i=-1,3
bsj(ik,i,ir)=bsj(ik,i,ir)*geospr(ir)
enddo
enddo
enddo
......
......@@ -13,17 +13,17 @@ c
double precision rr,z,pi,taunorm,rnow,rlast,tnow,tlast
double precision v00,depth,hpmin,vsliquid
double precision s1,s2,s,ds,smin,shead,twinmin
double precision r1,r2,dr,dm,t1,rdis
double precision r1,r2,dr,dm,t1,rdis
double precision mis,mcl,mdc,st,di,ra,deg2rad
double precision suppress,ros,vps,vss,fcut
double precision rot(3,3),sm(3,3),swap(3,3)
double precision rot(3,3),sm(3,3),swap(3,3)
double precision resolut(3),t0(nrmax)
character*80 outfile0(7)
c
c source parameters
c =================
c
pi=4.d0*datan(1.d0)
pi=4.d0*datan(1.d0)
deg2rad=pi/180.d0
call skip_comments(unit)
read(unit,*)zs
......@@ -201,75 +201,75 @@ c
if (ssel(7) .ne. 0) then
backspace(unit)
endif
if(ssel(7).eq.1)then
read(unit,*)ssel(7),(mtensor(i),i=1,6),outfile0(7)
else if(ssel(7).eq.2)then
read(unit,*)ssel(7),mis,mcl,mdc,st,di,ra,outfile0(7)
st=st*deg2rad
di=di*deg2rad
ra=ra*deg2rad
c
c use principal stress coordinates:
c x: along T-axis
c y: along N-axis
c z: along P-axis (symmetry axis of CLVD)
c
do i=1,3
do j=1,3
sm(i,j)=0.d0
enddo
enddo
sm(1,1)=mis-0.5d0*mcl+mdc
sm(2,2)=mis-0.5d0*mcl
sm(3,3)=mis+mcl-mdc
c
c construct the rotation matrix:
c 1. around x by -45 deg;
c 2. around z by rake angle
c 3. around x -dip angle;
c 4. around z by -strike angle.
c
rot(1,1)=(dcos(st)*dcos(ra)
& +dsin(st)*(dcos(di)*dsin(ra)-dsin(di)))/dsqrt(2.d0)
rot(1,2)=dcos(st)*dsin(ra)-dsin(st)*dcos(di)*dcos(ra)
rot(1,3)=(dcos(st)*dcos(ra)
& +dsin(st)*(dcos(di)*dsin(ra)+dsin(di)))/dsqrt(2.d0)
rot(2,1)=(dsin(st)*dcos(ra)
& -dcos(st)*(dcos(di)*dsin(ra)-dsin(di)))/dsqrt(2.d0)
rot(2,2)=dsin(st)*dsin(ra)+dcos(st)*dcos(di)*dcos(ra)
rot(2,3)=(dsin(st)*dcos(ra)
& -dcos(st)*(dcos(di)*dsin(ra)+dsin(di)))/dsqrt(2.d0)
rot(3,1)=(-dsin(di)*dsin(ra)-dcos(di))/dsqrt(2.d0)
rot(3,2)=dsin(di)*dcos(ra)
rot(3,3)=(-dsin(di)*dsin(ra)+dcos(di))/dsqrt(2.d0)
c
do i=1,3
do j=1,3
swap(i,j)=0.d0
do l=1,3
swap(i,j)=swap(i,j)+rot(i,l)*sm(l,j)
enddo
enddo
enddo
do i=1,3
do j=1,3
sm(i,j)=0.d0
do l=1,3
sm(i,j)=sm(i,j)+swap(i,l)*rot(j,l)
enddo
enddo
enddo
mtensor(1)=sm(1,1)
mtensor(2)=sm(2,2)
mtensor(3)=sm(3,3)
mtensor(4)=sm(1,2)
mtensor(5)=sm(2,3)
mtensor(6)=sm(3,1)
else
do i=1,6
mtensor(i)=0.d0
enddo
ssel(7)=0
if(ssel(7).eq.1)then
read(unit,*)ssel(7),(mtensor(i),i=1,6),outfile0(7)
else if(ssel(7).eq.2)then
read(unit,*)ssel(7),mis,mcl,mdc,st,di,ra,outfile0(7)
st=st*deg2rad
di=di*deg2rad
ra=ra*deg2rad
c
c use principal stress coordinates:
c x: along T-axis
c y: along N-axis
c z: along P-axis (symmetry axis of CLVD)
c
do i=1,3
do j=1,3
sm(i,j)=0.d0
enddo
enddo
sm(1,1)=mis-0.5d0*mcl+mdc
sm(2,2)=mis-0.5d0*mcl
sm(3,3)=mis+mcl-mdc
c
c construct the rotation matrix:
c 1. around x by -45 deg;
c 2. around z by rake angle
c 3. around x -dip angle;
c 4. around z by -strike angle.
c
rot(1,1)=(dcos(st)*dcos(ra)
& +dsin(st)*(dcos(di)*dsin(ra)-dsin(di)))/dsqrt(2.d0)
rot(1,2)=dcos(st)*dsin(ra)-dsin(st)*dcos(di)*dcos(ra)
rot(1,3)=(dcos(st)*dcos(ra)
& +dsin(st)*(dcos(di)*dsin(ra)+dsin(di)))/dsqrt(2.d0)
rot(2,1)=(dsin(st)*dcos(ra)
& -dcos(st)*(dcos(di)*dsin(ra)-dsin(di)))/dsqrt(2.d0)
rot(2,2)=dsin(st)*dsin(ra)+dcos(st)*dcos(di)*dcos(ra)
rot(2,3)=(dsin(st)*dcos(ra)
& -dcos(st)*(dcos(di)*dsin(ra)+dsin(di)))/dsqrt(2.d0)
rot(3,1)=(-dsin(di)*dsin(ra)-dcos(di))/dsqrt(2.d0)
rot(3,2)=dsin(di)*dcos(ra)
rot(3,3)=(-dsin(di)*dsin(ra)+dcos(di))/dsqrt(2.d0)
c
do i=1,3
do j=1,3
swap(i,j)=0.d0
do l=1,3
swap(i,j)=swap(i,j)+rot(i,l)*sm(l,j)
enddo
enddo
enddo
do i=1,3
do j=1,3
sm(i,j)=0.d0
do l=1,3
sm(i,j)=sm(i,j)+swap(i,l)*rot(j,l)
enddo
enddo
enddo
mtensor(1)=sm(1,1)
mtensor(2)=sm(2,2)
mtensor(3)=sm(3,3)
mtensor(4)=sm(1,2)
mtensor(5)=sm(2,3)
mtensor(6)=sm(3,1)
else
do i=1,6
mtensor(i)=0.d0
enddo
ssel(7)=0
endif
call skip_comments(unit)
read(unit,*)iazi
......@@ -308,7 +308,7 @@ c
read(unit,*)iflat
call skip_comments(unit)
read(unit,*)(resolut(i),i=1,3)
do i=1,3
do i=1,3
if(resolut(i).le.0.d0)resolut(i)=0.1d0
resolut(i)=1.d-02*resolut(i)
enddo
......@@ -583,8 +583,8 @@ c
c compare direct p, reflected p and head wave phase
c
lcut=0
do l=max0(ls,lzr),lp
c
do l=max0(ls,lzr),lp
c
twinmin=twindow
c
c 1. direct or reflected p wave
......@@ -658,17 +658,17 @@ c
enddo
smin=shead
hpmin=hp(l)
endif
c
endif
c
if(twinmin.lt.twindow)lcut=l
enddo
if(lcut.lt.1)then
stop ' time window too small!'
else if(lcut.lt.lp)then
lp=lcut
hp(lp)=0.d0
n0=nno(lp)
write(*,'(a,i3)')' actually used number of layers: ',n0
enddo
if(lcut.lt.1)then
stop ' time window too small!'
else if(lcut.lt.lp)then
lp=lcut
hp(lp)=0.d0
n0=nno(lp)
write(*,'(a,i3)')' actually used number of layers: ',n0
endif
c
c for partial solution only
......@@ -736,12 +736,12 @@ c
fsel(i,istp)=0
endif
enddo
enddo
c
c no toroidal component if ms = 0
c
do istp=1,6
if(ms(istp).eq.0)fsel(3,istp)=0
enddo
c
c no toroidal component if ms = 0
c
do istp=1,6
if(ms(istp).eq.0)fsel(3,istp)=0
enddo
c
calsh=ssel(2).eq.1.or.ssel(3).eq.1.or.ssel(6).eq.1
......
......@@ -8,7 +8,7 @@ c
integer nzmax,lmax,nrmax,nfmax,ndtransmax
parameter(lmax=1500)
parameter(nzmax=lmax+2)
parameter(nrmax=101,nfmax=2048)
parameter(nrmax=101,nfmax=2048)
parameter(ndtransmax=4)
c
c INDEX PARAMETERS FOR BESSEL FUNCTION TABLES
......@@ -88,29 +88,29 @@ c
double precision qp1rs(lmax),qp2rs(lmax),qs1rs(lmax),qs2rs(lmax)
common /rsmodel0/z1rs,z2rs,ro1rs,ro2rs,vp1rs,vp2rs,vs1rs,vs2rs,
+ qp1rs,qp2rs,qs1rs,qs2rs,l0rs
c
c layered model parameter:
c n0: number of homogeneous layers
c
integer n0
double precision h(lmax),ro(lmax),vp(lmax),vs(lmax)
double precision qp(lmax),qs(lmax)
common /model/ h,ro,vp,vs,qp,qs,n0
c
double complex acc(lmax),kp(lmax),ks(lmax),cla(lmax),cmu(lmax)
c
c layered model parameter:
c n0: number of homogeneous layers
c
integer n0
double precision h(lmax),ro(lmax),vp(lmax),vs(lmax)
double precision qp(lmax),qs(lmax)
common /model/ h,ro,vp,vs,qp,qs,n0
c
double complex acc(lmax),kp(lmax),ks(lmax),cla(lmax),cmu(lmax)
double complex cvp(lmax),cvs(lmax),wa(lmax),wb(lmax)
common /cpara/ acc,kp,ks,cla,cmu,cvp,cvs,wa,wb
c
c layered model parameter:
c n0rs: number of homogeneous layers (receiver site)
c
integer n0rs
double precision hrs(lmax),rors(lmax),vprs(lmax),vsrs(lmax)
double precision qprs(lmax),qsrs(lmax)
common /rsmodel/ hrs,rors,vprs,vsrs,qprs,qsrs,n0rs
c
c layered model parameter:
c n0rs: number of homogeneous layers (receiver site)
c
integer n0rs
double precision hrs(lmax),rors(lmax),vprs(lmax),vsrs(lmax)
double precision qprs(lmax),qsrs(lmax)
common /rsmodel/ hrs,rors,vprs,vsrs,qprs,qsrs,n0rs
c
double complex accrs(lmax),kprs(lmax),ksrs(lmax),
+ clars(lmax),cmurs(lmax)
+ clars(lmax),cmurs(lmax)
double complex cvprs(lmax),cvsrs(lmax),wars(lmax),wbrs(lmax)
common /rscpara/ accrs,kprs,ksrs,clars,cmurs,cvprs,cvsrs,wars,wbrs
c
......
......@@ -19,9 +19,9 @@ c
double complex cinc(4,6)
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
c
double complex c2
logical rsite
c
double complex c2
data c2/(2.d0,0.d0)/
c
ck=dcmplx(k,0.d0)
......@@ -34,26 +34,26 @@ c determination of starting upper sublayer
c
rsite=.false.
c
if(lup.eq.1.and.isurf.eq.0)then
do j=1,2
do i=1,4
yup(i,j)=(0.d0,0.d0)
enddo
if(lup.eq.1.and.isurf.eq.0)then
do j=1,2
do i=1,4
yup(i,j)=(0.d0,0.d0)
enddo
enddo
yup(1,1)=(1.d0,0.d0)
if(vpair.gt.0.d0)yup(2,1)=-accair/kpair
yup(3,2)=(1.d0,0.d0)
else
n=nno(lup)
c
yup(1,1)=kp(n)
yup(2,1)=wa(n)
yup(3,1)=ck
yup(4,1)=wb(n)*kp(n)
c
yup(1,2)=ck
yup(2,2)=wb(n)*ks(n)
yup(3,2)=ks(n)
n=nno(lup)
c
yup(1,1)=kp(n)
yup(2,1)=wa(n)
yup(3,1)=ck
yup(4,1)=wb(n)*kp(n)
c
yup(1,2)=ck
yup(2,2)=wb(n)*ks(n)
yup(3,2)=ks(n)
yup(4,2)=wa(n)
endif
if(lup.eq.lzr)call cmemcpy(yup,y0,8)
......@@ -69,13 +69,13 @@ c
swave=cdexp(-ks(n)*ch0)
c
c orthonormalization of the p-sv modes
c
cfac=(1.d0,0.d0)/(c0(3,2)*c0(1,1)-c0(1,2)*c0(3,1))
orth(1,1)=c0(3,2)*cfac
orth(1,2)=-c0(1,2)*cfac
orth(2,1)=-c0(3,1)*cfac
c
cfac=(1.d0,0.d0)/(c0(3,2)*c0(1,1)-c0(1,2)*c0(3,1))
orth(1,1)=c0(3,2)*cfac
orth(1,2)=-c0(1,2)*cfac
orth(2,1)=-c0(3,1)*cfac
orth(2,2)=c0(1,1)*cfac
call caxcb(c0,orth,4,2,2,c1)
call caxcb(c0,orth,4,2,2,c1)
c
if(l.gt.lzr)then
c
......@@ -142,7 +142,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,yup,c1,2,rsite)
c
if(l.eq.lzr)call cmemcpy(yup,y0,8)
enddo
......@@ -151,15 +151,15 @@ c===============================================================================
c
c matrix propagation from half-space to source
c
n=nno(llw)
ylw(1,1)=-kp(n)
ylw(2,1)=wa(n)
ylw(3,1)=ck
ylw(4,1)=-wb(n)*kp(n)
c
ylw(1,2)=ck
ylw(2,2)=-wb(n)*ks(n)
ylw(3,2)=-ks(n)
n=nno(llw)
ylw(1,1)=-kp(n)
ylw(2,1)=wa(n)
ylw(3,1)=ck
ylw(4,1)=-wb(n)*kp(n)
c
ylw(1,2)=ck
ylw(2,2)=-wb(n)*ks(n)
ylw(3,2)=-ks(n)
ylw(4,2)=wa(n)
c
if(llw.gt.ls.and.llw.eq.lzr)call cmemcpy(ylw,y0,8)
......@@ -175,13 +175,13 @@ c
swave=cdexp(-ks(n)*ch0)
c
c orthonormalization of the p-sv modes
c
cfac=(1.d0,0.d0)/(c0(4,2)*c0(2,1)-c0(2,2)*c0(4,1))
orth(1,1)=c0(4,2)*cfac
orth(1,2)=-c0(2,2)*cfac
orth(2,1)=-c0(4,1)*cfac
c
cfac=(1.d0,0.d0)/(c0(4,2)*c0(2,1)-c0(2,2)*c0(4,1))
orth(1,1)=c0(4,2)*cfac
orth(1,2)=-c0(2,2)*cfac
orth(2,1)=-c0(4,1)*cfac
orth(2,2)=c0(2,1)*cfac
call caxcb(c0,orth,4,2,2,c1)
call caxcb(c0,orth,4,2,2,c1)
c
if(l.lt.lzr)then
c
......@@ -295,15 +295,15 @@ c
endif
c
if(ipath.eq.1)then
n=nno(lpath)
ylw(1,1)=-kp(n)
ylw(2,1)=wa(n)
ylw(3,1)=ck
ylw(4,1)=-wb(n)*kp(n)
c
ylw(1,2)=ck
ylw(2,2)=-wb(n)*ks(n)
ylw(3,2)=-ks(n)
n=nno(lpath)
ylw(1,1)=-kp(n)
ylw(2,1)=wa(n)
ylw(3,1)=ck
ylw(4,1)=-wb(n)*kp(n)
c
ylw(1,2)=ck
ylw(2,2)=-wb(n)*ks(n)
ylw(3,2)=-ks(n)
ylw(4,2)=wa(n)
c
if(lpath.gt.ls.and.lpath.eq.lzr)call cmemcpy(ylw,y0,8)
......@@ -318,12 +318,12 @@ c
swave=cdexp(-ks(n)*ch0)
c
c orthonormalization of the p-sv modes
c
cfac=(1.d0,0.d0)/(c0(4,2)*c0(2,1)-c0(2,2)*c0(4,1))
orth(1,1)=c0(4,2)*cfac
orth(1,2)=-c0(2,2)*cfac
orth(2,1)=-c0(4,1)*cfac
orth(2,2)=c0(2,1)*cfac
c
cfac=(1.d0,0.d0)/(c0(4,2)*c0(2,1)-c0(2,2)*c0(4,1))
orth(1,1)=c0(4,2)*cfac
orth(1,2)=-c0(2,2)*cfac
orth(2,1)=-c0(4,1)*cfac
orth(2,2)=c0(2,1)*cfac
c
call caxcb(c0,orth,4,2,2,c1)
if(l.lt.lzr)then
......@@ -428,7 +428,7 @@ c
if(n0rs.gt.0)then
c
c for receiver-site structure different from source-site structure
c
c
n=nno(lzr)
do istp=1,6
call qsve2am(n,ck,y(1,istp),cinc(1,istp),1,rsite)
......@@ -438,26 +438,26 @@ c determination of starting upper sublayer
c
rsite=.true.
c
if(isurf.eq.0)then
do j=1,2
do i=1,4
yup(i,j)=(0.d0,0.d0)
enddo
if(isurf.eq.0)then
do j=1,2