#!/bin/sh cat > Makefile << EOF MCHDEP = mchdep70.o MCHCMD = numarg.o DSPEC8 = dspec832.o FCMP = g77 -O3 CCMP = gcc .c.o: \$(CCMP) -c \$< .f.o: \$(FCMP) \$(FFLAG) -c \$< all: hspec9 rhfoc10 dsplt10 clean hspec9: hspec9.o \$(MCHDEP) \$(FCMP) \$(FFLAG) -o hspec9 hspec9.o \$(MCHDEP) rm hspec9.o rhfoc10: rhfoc10.o \$(MCHCMD) \$(MCHDEP) \$(FCMP) \$(FFLAG) -o rhfoc10 rhfoc10.o \$(MCHDEP) \$(MCHCMD) rm rhfoc10.o dsplt10: dsplt10.o \$(MCHCMD) \$(MCHDEP) \$(FCMP) \$(FFLAG) -o dsplt10 dsplt10.o \$(MCHDEP) \$(MCHCMD) rm dsplt10.o clean: rm *.o EOF cat > numarg.f << EOF function numarg() integer numarg numarg = iargc() return end EOF cat > numarg_.c << EOF long int numarg_() { extern int xargc; return ((long)xargc -1); } EOF cat > hspec9.f << EOF program hspec9 c---------------------------------------------------------------------c c c c COMPUTER PROGRAMS IN SEISMOLOGY c c VOLUME VI c c c c PROGRAM: HSPEC9 c c This program is a merge of HSPEC8 and RHWVINTA which c does away with the intermediate file c c c COPYRIGHT 1986 c c R. B. Herrmann c c Department of Earth and Atmospheric Sciences c c Saint Louis University c c 221 North Grand Boulevard c c St. Louis, Missouri 63013 c c U. S. A. c c c c---------------------------------------------------------------------c parameter(LER=0, LIN=5, LOT=6) parameter(NDST=100) parameter(NL=100) common/source/depth,lmax,dph ,vamin,vamax,vbmin,vbmax,hvert common/damp/alpha common/model/d(NL),a(NL),b(NL),rho(NL),mmax,qa(NL),qb(NL) common/jout/jsrc(10) , jbdry character*50 name2 character*50 names character*50 names1, names2 character*3 istat3 logical ixst3 real*4 ffreq(8) real*4 r(NDST), t0(NDST) real*4 ar(10),ai(10) character*4 icchar(11) common/c/cmax,c1,c2,cmin integer*4 lsrc(10) data lsrc/1,2,3,4,9,5,6,10,7,8/ data icchar/' Z1 ',' R1 ',' Z2 ',' R2 ',' T2 ', 1 ' Z3 ',' R3 ',' T3 ',' Z4 ',' R4 ',' XX '/ c----- c call machine dependent initialization c----- call mchdep() c----- c read name of hspec8 data file c----- read(LIN,'(a)')names call gethsp(names,name2,fl,fu,dt,n,n1,n2,df,nyq2,xleng,xfac) c----- c----- c read in rhwvinta data file c----- read(LIN,'(a)')names call getrhw(names,r,t0,ndist,names1,names2) c----- c process c----- df = 1./(n*dt) nyq = nyq2/2 write(LOT,2) fl,fu,df,n1,n2,depth,n 2 format(4hfl =,f10.5,5x,4hfu =,f10.5,5x,4hdf =,f10.5,/ 1 4x,4hn1 =,i4,5x,4hn2 =,i4,5x,7hdepth =,f10.2,4h n =,i5) write(LOT,5)alpha,dt 5 format(7halpha =,f10.5,5x,4hdt =,f10.3) write(LOT,4) 4 format('frequencies for which response computed ') c----- c open output file for hspec8 c----- open(unit=2,file='hspec9.tmp',status='unknown',form= 1 'unformatted',access='sequential') rewind 2 c----- c open temporary distance multiplexed file which c will be examined in case it exists already for c work already done c----- inquire(file=names1,exist=ixst3) if(ixst3)then istat3 = 'old' else istat3 = 'new' endif open(unit=3,file=names1,status=istat3,form= 1 'unformatted',access='sequential') rewind 3 ifreq = 0 if(ixst3)then call reset3(ifreq,jsrc,lsrc,ndist,n1,n2) endif do 101 i=1,8 101 ffreq(i)=-1.0 n11 = n1 if(ifreq.gt.n1)n11 = ifreq + 1 do 100 i = n11,n2 freq=(i-1)*df if(freq.lt.df) freq = 0.01*df rewind 2 call excit(freq,xleng,xfac) do 200 jd=1,ndist rewind 2 call setup(depth,r(jd)) call wvint(r(jd),ar,ai,depth) fac = 6.2831853*freq*t0(jd) xr = cos(fac) xi = sin(fac) do 301 jj=1,10 datar=xr*ar(jj)-xi*ai(jj) datai=xr*ai(jj)+xi*ar(jj) if(jsrc(lsrc(jj)).eq.1)write(3)datar,datai 301 continue 200 continue index=mod(i,8) if(index.eq.0)index=8 ffreq(index)=freq if (index.eq.8) then write(LOT,3)ffreq do 102 ii=1,8 102 ffreq(ii)=-1. endif 3 format(8f10.5) 100 continue write(LOT,3)ffreq 9999 continue c----- c output the final spectrum as a function of distance c----- open(unit=4,file=names2,status='unknown',form='unformatted', 1 access='sequential') rewind 4 write(4) alpha,depth,fl,fu,dt,n,n1,n2,df,nyq2 write(4)jsrc c-----write out lsrc indexing since rhwvinta changes order write(4)lsrc write(4)d,a,b,rho,mmax,qa,qb c----- c no output the spectrum for each distance c----- do 5000 jd=1,ndist write(4)r(jd),t0(jd) do 5100 jj=1,10 5100 write(4)icchar(jj) rewind 3 do 5200 i=n1,n2 do 5300 jjd=1,ndist do 5400 jj=1,10 if(jsrc(lsrc(jj)).eq.1)then read(3)datar,datai if(jjd.eq.jd)then write(4)datar,datai endif endif 5400 continue 5300 continue 5200 continue 5000 continue rr = -1.0 tt0 = 0.0 write(4)rr,tt0 close (4) close(2,status='delete') close(3,status='delete') stop end subroutine gethsp(names,name2,fl,fu,dt,n,n1,n2,df,nyq2,xleng,xfac) character names*(*) c----- c read in data file of hspec8 commands c----- parameter(LER=0, LIN=5, LOT=6) parameter(NL=100) common/source/depth,lmax,dph ,vamin,vamax,vbmin,vbmax,hvert common/damp/alpha common/model/d(NL),a(NL),b(NL),rho(NL),mmax,qa(NL),qb(NL) common/jout/jsrc(10) , jbdry character*50 name2 character*3 istat2 open(2,file=names,status='old',form='formatted', 1 access='sequential') rewind 2 read(2,18)name2 read(2,18)istat2 read(2,19)ierr,ifreq 18 format(a) 19 format(2i5) read(2,20) alpha,depth,fl,fu,dt,n,n1,n2,df,nyq2 ,mmax 20 format(5e15.7/3i10,e15.7,i10,i10) read(2,21)jsrc , jbdry 21 format(11i5) read(2,22)(d(i),a(i),b(i),rho(i),qa(i),qb(i),i=1,mmax) 22 format(6e11.4) read(2,23)lmax,dph,vamin,vamax,vbmin,vbmax,hvert 23 format(i10,6e11.4) read(2,24)xleng, xfac 24 format(2e15.7) close (2) return end subroutine getrhw(names,r,t0,ndist,names1,names2) parameter(LER=0, LIN=5, LOT=6) character names*(*) character names1*(*), names2*(*) real*4 r(1), t0(1) integer*4 ndist common/c/cmax,c1,c2,cmin c----- c get command parameters from rhwvinta.dat c----- open(unit=2,file=names,status='unknown',form='formatted', 1 access='sequential') rewind 2 read(2,'(a)')names1 read(2,'(a)')names2 read(2,1)cmax,c1,c2,cmin 1 format(4f10.5) ndist = 0 500 continue read(2,1)rr,tshift,vred if(rr.le.0.0)goto 9998 ndist = ndist + 1 r(ndist)=rr if(vred.eq.0.0)vred = 1.0e+20 t0(ndist) = tshift + r(ndist)/vred if(r(ndist).le.0.0)goto 9998 goto 500 9998 continue close (2) return end subroutine reset3(ifreq,jsrc,lsrc,ndist,n1,n2) c----- c routine to reset temporary output file on unit 3 c if it already exists due to an aborted run c c the file is read until an error is found, which c indicates the total number of correct frequencies on the c output file. The file is rewound and the correct frequencies c are read in to reposition the output file. The c parameter ifreq indicates the number of complete frequency sets c----- integer*4 lsrc(1), jsrc(10) c----- c find the correct number by duplicating the writes of the main program c----- ifreq = 0 rewind 3 do 1000 i = n1, n2 do 1100 jd=1,ndist do 1200 jj=1,10 if(jsrc(lsrc(jj)).eq.1)read(3,err=9998,end=9999)xx,yy 1200 continue 1100 continue ifreq = ifreq + 1 1000 continue 9998 continue 9999 continue c----- c there are now ifreq known data sets out there c position write pointer on the output file c----- rewind 3 do 5000 i = 1,ifreq do 5100 jd=1,ndist do 5200 jj=1,10 if(jsrc(lsrc(jj)).eq.1)read(3,err=9998,end=9999)xx,yy 5200 continue 5100 continue 5000 continue ifreq = ifreq -1 + n1 return end subroutine aten(omega,qa,qb,xka,xkb,alpha,a,b,atna,atnb) real*4 qa,qb,alpha,a,b double complex omega,at,atna,atnb,xka,xkb real*8 pi, om1, oml, fac c reference frequency is 1.0 hz om1=6.2831853d+00 c low frequency cutoff is 0.01 hz oml=0.062831853d+00 pi=3.1415927d+00 at=dcmplx(0.0d+00,0.0d+00) if(zabs(omega).gt.oml) at=zlog(omega/om1)/pi if(zabs(omega).gt.oml) go to 40 fac=dsqrt(oml*oml + dble(alpha*alpha))/oml at=zlog(dcmplx(dble(oml),-dble(alpha))/(om1*fac))/pi 40 continue atna=dcmplx(1.0d+00,0.0d+00) atnb=dcmplx(1.0d+00,0.0d+00) if(qa.gt.0.0) atna=(1.+dble(qa)*at+dcmplx(0.0d+00,dble(qa/2.))) if(qb.gt.0.0) atnb=(1.+dble(qb)*at+dcmplx(0.0d+00,dble(qb/2.))) xka=omega/(dble(a)*atna) xkb=omega/(dble(b)*atnb) return end subroutine cmult(e,ca,exa) double complex ca(5,5) real*8 exa real *8 xnorm double complex e(5) double complex c, ee(5) do 1350 i=1,5 c = dcmplx(0.0d+00,0.0d+00) do 1349 j=1,5 c=c+ e(j) * ca(j,i) 1349 continue ee(i)=c 1350 continue call normc(ee,exa,xnorm) do 1351 i=1,5 e(i) = ee(i)*xnorm 1351 continue return end subroutine dmult(da,aa) double complex aa(4,4) double complex sumd,ea(4,4),da(4,4) do 1360 i=1,4 do 1361 j=1,4 sumd = dcmplx(0.0d+00,0.0d+00) do 1362 jj=1,4 sumd=sumd+da(i,jj) * aa(jj,j) 1362 continue ea(i,j)=sumd 1361 continue 1360 continue do 1363 j=1,4 do 1364 i=1,4 da(i,j)=ea(i,j) 1364 continue 1363 continue return end subroutine dnka(ca,wvno2,gam,rho) double complex gam,ca,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz,gam2, 1 gamm1,gamm2,a0c,xz2,wy2,temp common/ ovrflw / a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz real *8 a0 double complex wvno2 double complex cqww2, cqxw2, g1wy2, gxz2, g2wy2, g2xz2 double complex gg1, a0cgg1 dimension ca(5,5) double complex zrho, zrho2 zrho = dcmplx(dble(rho),0.0d+00) zrho2= dcmplx(dble(rho*rho),0.0d+00) gam2 = gam*gam gamm1 = gam-1. gamm2 = gamm1*gamm1 cqww2 = cqw * wvno2 cqxw2 = cqx / wvno2 gg1 = gam*gamm1 a0c = dcmplx(2.0d+00,0.0d+00)*(dcmplx(a0,0.0d+00)-cpcq) xz2 = xz/wvno2 gxz2 = gam*xz2 g2xz2 = gam2 * xz2 a0cgg1 = a0c*(gam+gamm1) wy2 = wy*wvno2 g2wy2 = gamm2 * wy2 g1wy2 = gamm1 * wy2 temp = wvno2*(a0c + wy2) + xz ca(1,5) = -temp/zrho2 temp = dcmplx(0.5d+00,0.0d+00)*a0cgg1 + gxz2 + g1wy2 ca(1,3) = -temp/zrho temp = a0c*gg1 + g2xz2 + g2wy2 ca(3,3) = a0 + temp + temp ca(1,1) = cpcq-temp temp =dcmplx(0.5d+00,0.0d+00)*a0cgg1*gg1 + gam2*gxz2 + gamm2*g1wy2 ca(3,1) = dcmplx(2.0d+00,0.0d+00)*temp*zrho temp = gamm2*(a0c*gam2 + g2wy2) + gam2*g2xz2 ca(5,1) = -zrho2*temp/wvno2 ca(1,4) = (-cqww2+cpz)/zrho ca(2,1) = (-gamm2*cqw + gam2*cpz/wvno2)*zrho ca(2,3) = -(gamm1*cqww2 - gam*cpz)/wvno2 ca(1,2) = (-cqx + wvno2*cpy)/zrho ca(3,2) = wvno2*(gam*cqxw2 - gamm1*cpy)*dcmplx(2.0d+00,0.0d+00) ca(4,1) = (-gam2*cqxw2 + gamm2*cpy)*zrho ca(2,2) = cpcq ca(2,4) = -wz ca(4,2) = -xy ca(2,5)=ca(1,4) ca(5,5)=ca(1,1) ca(5,4)=ca(2,1) ca(5,3)=-ca(3,1)/(dcmplx(2.0d+00,00d+00)*wvno2) ca(5,2)=ca(4,1) ca(4,3)= -ca(3,2)/(dcmplx(2.0d+00,00d+00)*wvno2) ca(4,5)=ca(1,2) ca(4,4)=ca(2,2) ca(3,4)=-dcmplx(2.0d+00,00d+00)*wvno2*ca(2,3) ca(3,5)=-dcmplx(2.0d+00,00d+00)*wvno2*ca(1,3) return end subroutine excit(freq,xleng,xfac) c----- c sample response for all wavenumbers at a given frequency c using Bouchon equal wavenumber sampling = dk c with offset of 0.218dk c----- parameter(LER=0,LIN=5,LOT=6) parameter(NL=100) complex gg(10) common/source/depth,lmax,dph ,vamin,vamax,vbmin,vbmax,hvert common/model/d(NL),a(NL),b(NL),rho(NL),mmax,qa(NL),qb(NL) common/jout/jsrc(10) , jbdry common/damp/alpha double complex wvn,om c----- c setup wavenumber sampling, choosing two large values to c control the low frequency asymptotic integration technique c c we will permit a depth = 0, and change the sampling in this c case as a fraction of thickness of the first layer c----- do 10 i=mmax-1,1,-1 if(d(i).gt.0.0)deep=d(i) 10 continue if(depth.lt. 0.1*deep)then dpth = 0.1*deep else dpth = depth endif omega=6.2831853*freq wvbm = omega/vbmin dk = 6.2831853/xleng wvmm = (5.0/dpth) + xfac*wvbm wvzmx = wvmm * depth nk = wvmm / dk mk=nk+2 mk1=nk+1 write(2)omega,mk c-----output wavenumber in reverse order c also output first two wavenumbers to control c low frequency asymptotic integration c-----to save space in sequential i/o buffer output stream call bufini(1,ierr) do 3998 ii=mk,1,-1 if(ii.eq.mk)then if(wvzmx.gt.5.0)then wv = 6.0/dpth else wv = wvmm + 10.0*dk endif elseif(ii.eq.mk1)then if(wvzmx.gt.5.0)then wv = 2.5/dpth else wv = wvmm + 5.0*dk endif else wv = (ii-1)*dk + 0.218*dk endif wvn=dcmplx(dble(wv),0.0d+00) om=dcmplx(dble(omega),-dble(alpha)) call rshof(gg,om,wvn,jbdry) call bufwr(wv) do 3998 j=1,10 if(jsrc(j).eq.1)then call bufwr(real(gg(j))) call bufwr(aimag(gg(j))) endif 3998 continue call buflsh return end subroutine hska(aa,w,x,y,z,cosp,cosq,wvno2,gam,gamm1,rho) double complex wvno2 double complex aa(4,4),w,x,y,z,cosp,cosq,gam,gamm1 double complex cpq, gcpq, zw2, gzw2, g1w, g1y, gx real*8 drho drho = dble(rho) cpq = cosp-cosq gcpq = gam*cpq zw2 = z/wvno2 gzw2 = gam*zw2 g1w = gamm1*w g1y = gamm1*y gx = gam*x aa(1,1) = gcpq + cosq aa(1,3) = -cpq/dcmplx(drho,0.0d+00) aa(1,2)=-g1w+gzw2 aa(1,4)=(wvno2*w-z)/dcmplx(drho,0.0d+00) aa(2,1)= gx - wvno2*g1y aa(2,2) = -gcpq + cosp aa(2,3)=(-x+wvno2*y)/dcmplx(drho,0.0d+00) aa(2,4)= - wvno2*aa(1,3) aa(3,1)= dcmplx(drho,0.0d+00)*gamm1*gcpq aa(3,2)=dcmplx(drho,0.0d+00)*((-gamm1*g1w)+(gam*gzw2)) aa(3,4)=-wvno2*aa(1,2) aa(3,3)=aa(2,2) aa(4,1)=dcmplx(drho,0.0d+00)*(((gam*gx)/wvno2) - (gamm1*g1y)) aa(4,2)=-aa(3,1)/wvno2 aa(4,3)=-aa(2,1)/wvno2 aa(4,4)=aa(1,1) return end subroutine lmult(d11,d12,cosql,yl,zl,rho,b,atnb) double complex d11,d12,cosql,yl,zl,atnb,h,e1,e2 h= dcmplx(dble(rho*b*b),0.0d+00)*atnb*atnb yl=yl/h zl=zl*h e1=d11 e2=d12 d11=e1*cosql + e2*zl d12=e1*yl + e2*cosql return end subroutine normc(e,ex,xnorm) real*8 ex real *8 test,testt,x,y,fac,xnorm double complex e(*) test = 0.0D+00 testt = 0.0D+00 do 2 i = 1,5 if(dabs(dreal(e(i))).gt.testt) testt =dabs(dreal(e(i))) if(dabs(dimag(e(i))).gt.testt) testt =dabs(dimag(e(i))) 2 continue if(testt.lt.1.0e-30)testt=1.0 do 1 i =1,5 x=dreal(e(i))/testt y=dimag(e(i))/testt fac = x*x + y*y if(test.lt.fac) test = fac 1 continue test = testt*dsqrt(test) if(test.lt.1.0d-30) test=1.0 xnorm = 1./test ex =-dlog(xnorm) return end subroutine rshof(gg,om,wvno,jbdry) parameter(NL=100) common/model/d(NL),a(NL),b(NL),rho(NL),mmax,qa(NL),qb(NL) common/source/depth,lmax,dph ,vamin,vamax,vbmin,vbmax,hvert common/damp/alpha double complex ggg(10) double complex wvno,wvno2,wvno3 double complex cd(5),da(4,4),fr,y(4,2) complex gg(10) double complex om,fourpi,ka2,kb2 double complex d11,d12,fl double complex s12,s21,s32,s14,s23,s34,s32e,s34e double complex atna,atnb double complex wv4pi real *8 fact,facx,exe,exl,exel,exll,elj fourpi=12.5663706d+00*om*om if(zabs(wvno).lt.1.0d-8) go to 100 wv4pi = 2.0d+00 * wvno / fourpi wvno2=wvno*wvno wvno3 = wvno2*wvno call aten(om,qa(lmax),qb(lmax),ka2,kb2,alpha,a(lmax),b(lmax), 1 atna,atnb) ka2=ka2*ka2 kb2=kb2*kb2 c for frequencies less than 0.02 hz, use haskell formulation c otherwise use the dunkin formulation ifreq=1 call scoef(cd,da,fr,om,wvno,exe,exl, 1 fl,d11,d12,exel,exll,ifreq,jbdry) do 50 j=1,2 y(1,j)= cd(1)*da(2,j) + cd(2)*da(3,j) - wvno2*(cd(3)*da(4,j)) y(2,j)=-cd(1)*da(1,j) + cd(3)*da(3,j) + cd(4)*da(4,j) y(3,j)=-cd(2)*da(1,j) - cd(3)*da(2,j) + cd(5)*da(4,j) y(4,j)= wvno2*(cd(3)*da(1,j)) - cd(4)*da(2,j) - cd(5)*da(3,j) 50 continue s12=dcmplx(0.0d+00,0.0d+00) s14=-wv4pi s21=2.0d+00*kb2/(dble(rho(lmax))*fourpi) s23=s12 s32=wv4pi*2.*ka2/dble(rho(lmax)) s34=wv4pi*dble( (2.*b(lmax)/a(lmax))**2 - 3.) s32e=ka2*wv4pi/dble(rho(lmax)) s34e=2.0d+00*wv4pi*(ka2/kb2) do 60 j=1,2 ggg(j)=s32*y(2,j) + s34*y(4,j) ggg(j+2)=s21*y(1,j) + s23*y(3,j) ggg(j+4)=s12*y(2,j) + s14*y(4,j) ggg(j+6)=s32e*y(2,j)+s34e*y(4,j) 60 continue fact = 0.0D+00 elj=exe-exl if(elj.lt.55.) fact=dexp(-elj) do 70 j=1,8 70 gg(j) = (-ggg(j) * fact/fr) facx = 1./(12.5663706*b(lmax)*b(lmax)) gg(9) = 2.0d+00*d11/dble(rho(lmax)) gg(10) = -2.0d+00*wvno*d12*atnb*atnb*dble(b(lmax)*b(lmax)) elj=exll-exel fact = 0.0 if(elj.gt.-40.) fact = dexp(elj) facx = fact*facx gg(9) =(gg(9)*facx) /(fl*(atnb*atnb)) gg(10)=(gg(10)*facx)/(fl*(atnb*atnb)) return 100 continue do 101 i = 1,10 101 gg(i) = cmplx(0.0,0.0) return end subroutine scoef(cd,da,fr,omega,wvno,exe,exl, 1 fl,d11,d12,exel,exll,ifreq,jbdry) parameter(NL=100) common/model/d(NL),a(NL),b(NL),rho(NL),mmax,qa(NL),qb(NL) common/source/depth,lmax,dph ,vamin,vamax,vbmin,vbmax,hvert common/damp/alpha double complex omega,xka,xkb,ra,rb,gam,gamm1,p,q,w,x,y double complex cosq,z,cosp double complex da,ca,aa dimension da(4,4),ca(5,5),aa(4,4) double complex cd(5),e(5),fr double complex atna,atnb double complex d11,d12,e1,e2,fl double complex yl,zl,cosql real *8 exe,exl,exel,exll,ex,exa,exb double complex wvno,wvno2 c----- c initialize matrices c----- exe=0.0 exl=0.0 do 2 i = 1,4 do 3 j = 1,4 da(i,j)=dcmplx(0.0d+00,0.0d+00) 3 continue da(i,i) = dcmplx(1.0d+00,0.0d+00) 2 continue exel = 0.0 exll = 0.0 c----- c set up halfspace conditions c----- wvno2=wvno*wvno call aten(omega,qa(mmax),qb(mmax),xka,xkb,alpha,a(mmax), 1 b(mmax),atna,atnb) ra=zsqrt(wvno2-xka*xka) rb=zsqrt(wvno2-xkb*xkb) gam = dble(b(mmax))*wvno/omega gam = gam * atnb gam = 2.0d+00 * (gam * gam) gamm1 = gam - dcmplx(1.0d+00,0.0d+00) c----- c set up halfspace boundary conditions c c jbdry = -1 RIGID c = 0 ELASTIC c = +1 FREE SURFACE c c----- if(jbdry.eq.0)then e(1)=dble(rho(mmax)*rho(mmax))* 1 (-gam*gam*ra*rb+wvno2*gamm1*gamm1) e(2)=-dble(rho(mmax))*(wvno2*ra) e(3)=dble(rho(mmax))*(-gam*ra*rb+wvno2*gamm1) e(4)=dble(rho(mmax))*(wvno2*rb) e(5)=wvno2*(wvno2-ra*rb) e1 = dble(rho(mmax))*rb e2 = dcmplx(1.0d+00,0.0d+00) 1 /(dble(b(mmax)*b(mmax))*atnb*atnb) elseif(jbdry.lt.0)then e(1) = dcmplx(1.0d+00,0.0d+00) e(2) = dcmplx(0.0d+00,0.0d+00) e(3) = dcmplx(0.0d+00,0.0d+00) e(4) = dcmplx(0.0d+00,0.0d+00) e(5) = dcmplx(0.0d+00,0.0d+00) e1 = dcmplx(1.0d+00,0.0d+00) e2 = dcmplx(0.0d+00,0.0d+00) elseif(jbdry.gt.0)then e(1) = dcmplx(0.0d+00,0.0d+00) e(2) = dcmplx(0.0d+00,0.0d+00) e(3) = dcmplx(0.0d+00,0.0d+00) e(4) = dcmplx(0.0d+00,0.0d+00) e(5) = dcmplx(1.0d+00,0.0d+00) e1 = dcmplx(0.0d+00,0.0d+00) e2 = dcmplx(1.0d+00,0.0d+00) endif do 11 i=1,5 cd(i)=e(i) 11 continue d11=e1 d12=e2 c----- c matrix multiplication from bottom layer upward c----- mmm1 = mmax - 1 do 1340 k = 1,mmm1 m = mmax - k call aten(omega,qa(m),qb(m),xka,xkb,alpha,a(m),b(m), 1 atna,atnb) gam=dble(b(m))*(wvno/omega) gam = gam * atnb gam = dcmplx(2.0d+00,0.0d+00)*gam*gam gamm1 = gam - dcmplx(1.0d+00,0.0d+00) ra=zsqrt(wvno2-xka*xka) rb=zsqrt(wvno2-xkb*xkb) p=ra*dble(d(m)) q=rb*dble(d(m)) dpth=d(m) call var(p,q,ra,rb,w,x,y,z, 1 cosp,cosq,ex,exa,exb,yl,zl,cosql) call dnka(ca,wvno2,gam,rho(m)) exe=exe+exa call cmult(e,ca,exa) exe=exe + exa exel = exel + exb call lmult(e1,e2,cosql,yl,zl,rho(m),b(m),atnb) l1=lmax+1 if(l1.eq.m) then do 1352 i=1,5 cd(i)=e(i) 1352 continue exl=exe exll = exel d11=e1 d12=e2 endif c----- c source layer contributions c----- if(m.eq.lmax)then dpth=dph p=ra*dble(dpth) q=rb*dble(dpth) call var(p,q,ra,rb,w,x,y,z, 1 cosp,cosq,ex,exa,exb,yl,zl,cosql) call dnka(ca,wvno2,gam,rho(lmax)) exl=exl+exa call cmult(cd,ca,exa) exl=exl+exa call lmult(d11,d12,cosql,yl,zl,rho(lmax),b(lmax),atnb) exll=exll+exb dpth=d(lmax)-dph p=ra*dble(dpth) q=rb*dble(dpth) call var(p,q,ra,rb,w,x,y,z, 1 cosp,cosq,ex,exa,exb,yl,zl,cosql) call hska(aa,w,x,y,z,cosp,cosq,wvno2,gam, 1 gamm1,rho(lmax)) call dmult(da,aa) exl=exl+ex endif if(m.lt.lmax) then call hska(aa,w,x,y,z,cosp,cosq,wvno2,gam,gamm1,rho(m)) call dmult(da,aa) exl=exl+ex endif 1340 continue fl=e1 fr=e(1) return end subroutine var(p,q,ra,rb,w,x,y,z,cosp,cosq,ex, 1exa,exl,yl,zl,cosql) c not modified for negative p,q c this assumes that real p and real q have same signs common/ovrflw/a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz double complex cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz double complex p,q,ra,rb,w,x,y,z,cosp,cosq double complex yl,zl,cosql double complex eqp,eqm,epp,epm,sinp,sinq real *8 a0,pr,pi,qr,qi,fac,qmp,ex,exa,exl ex=0.0d+00 exl = 0.0d+00 a0=0.0d+00 pr=dreal(p) pi=dimag(p) qr=dreal(q) qi=dimag(q) epp=dcmplx(dcos(pi),dsin(pi))/2. epm=dconjg(epp) eqp=dcmplx(dcos(qi),dsin(qi))/2. eqm=dconjg(eqp) ex=pr exl=qr fac=0.0 if(pr.lt.15.) fac=dexp(-2.*pr) cosp=epp + fac*epm sinp=epp - fac*epm w=sinp/ra x=ra*sinp fac=0.0d+00 if(qr.lt.15.) fac=dexp(-2.*qr) cosql=eqp + fac*eqm sinq=eqp - fac*eqm yl=sinq/rb zl=rb*sinq exa=pr + qr cpcq=cosp*cosql cpy=cosp*yl cpz=cosp*zl cqw=cosql*w cqx=cosql*x xy=x*yl xz=x*zl wy=w*yl wz=w*zl fac=0.0d+00 qmp=qr-pr if(qmp.gt.-40.) fac=dexp(qmp) cosq=cosql*fac y=fac*yl z=fac*zl fac=0.0d+00 if(exa.lt.60.) a0=dexp(-exa) return end subroutine bufini(irdwr,ierr) c------initialize buffer pointer c------irdwr = 0 read initialize c------irdwr = 1 write initialize integer BUFMAX parameter(BUFMAX=2000) common/buf/iptr,max,buffer(BUFMAX) c save /buf/ iptr = 1 if(irdwr.eq.0)call getbuf(ierr) return end subroutine buflsh c------flush output buffer integer BUFMAX parameter(BUFMAX=2000) common/buf/iptr,max,buffer(BUFMAX) c save /buf/ ipt = iptr -1 if(ipt.gt.0)write(2)ipt,(buffer(i),i=1,ipt) iptr = 1 return end subroutine bufwr(x) c------fill buffer with floating point variable x, c------flush buffer as necessary integer BUFMAX parameter(BUFMAX=2000) common/buf/iptr,max,buffer(BUFMAX) c save /buf/ buffer(iptr) = x iptr = iptr + 1 if(iptr.gt.BUFMAX)call buflsh return end subroutine getbuf(ierr) c------read in file contents into buffer, taking care not to c------read beyond the contents of the file integer BUFMAX parameter(BUFMAX=2000) common/buf/iptr,max,buffer(BUFMAX) c save /buf/ c------ierr = 0 successful read c------ = 1 read error c------ = 2 end of file c------ read(2,err=1000,end=2000)max,(buffer(i),i=1,max) iptr = 1 ierr = 0 return 1000 ierr = 1 return 2000 ierr = 2 return end subroutine bufrd(x,ierr) c-----retrieve a value from buffer array, red in new array c-----as necessary c-----iptr is here the next array element to be read c-----it is always >= 1. We do not worry the upper limit c-----since the calling program must worry about this c-----because read always follows a complete write integer BUFMAX parameter(BUFMAX=2000) common/buf/iptr,max,buffer(BUFMAX) c save /buf/ c only yank in new data if actually required if(iptr.gt.max)call getbuf(ierr) x = buffer(iptr) iptr = iptr + 1 return end subroutine wvint(r,ar,ai,depth) c-----to work with potentially large disk files, we cannot read in c-----all wavenumbers at once. We only work with neighboring c-----points at any time. The first two are for the DC correction, c-----followed by wavenumbers in decreasing order common/c/cmax,c1,c2,cmin common/jout/jsrc(10) , jbdry common/fct/fact ,fact0 common/asym/j0k0,j0k1,j0k2,j0k3,j1k0,j1k1,j1k2,j1k3, 1 j2k0,j2k1,j2k2,j2k3 real j0k0,j0k1,j0k2,j0k3,j1k0,j1k1,j1k2,j1k3,j2k0,j2k1, 1 j2k2,j2k3 dimension ar(10),ai(10) complex aa(10),bb(10),cc(10) complex gg1,sumc complex g(10,2) complex smm(10) complex sumd dimension wvn(2) real j0,j1,j01,j11 logical iasymp read(2) omega,nk call bufini(0,ierr) nkm2 = nk -2 c choose proper strategy for integration, e.g. c whether or not to use asymptotic technique c-----get the last wavenumber and g values call getgk(g,2,jsrc,wvn(2)) c-----get the second last values call getgk(g,1,jsrc,wvn(1)) c----- c find asymptotic coefficients using these values c----- nlow = 1 nupp = 2 do 102 j=1,10 call solu(g(j,nlow),g(j,nupp),wvn(nlow),wvn(nupp), 1 depth,j,aa(j),bb(j),cc(j) ) 102 continue c----- c we no longer need the first value to shift wvnosv = wvn(2) call shift(g,wvn) c----- c----- c get first usable wavenumber which is the largest c----- call getgk(g,1,jsrc,wvn(1)) c----- c decide whether to invoke asymptotic correction c----- if(wvn(1).gt.wvnosv)then iasymp = .false. else iasymp = .true. endif c----- c initialize integral c----- call intini(smm,aa,bb,cc,iasymp,r) c----- c we now can procede with integration c----- c----- c in the variables below the t0,j0,sum refer to upper limit c of integration and t01,j01,j11 and sum1 refer to the lower limit c----- wvcm=omega*cmax wvc1=omega*c1 wvc2=omega*c2 wvcn=omega*cmin wvno = wvn(1) t0 =wvno*r call hank(t0,j0,j1) if(iasymp)then call gasym(g,aa,bb,cc,wvn,1,depth,jsrc) endif c----- c perform wavenumber filtering c----- call wvfilt(cmax,wvcm,wvc1,wvc2,wvcn,wvno,fact) do 200 ik = nkm2,1,-1 i = 2 if(ik.lt.nkm2)then call shift(g,wvn) call getgk(g,1,jsrc,wvn(1)) endif c------to avoid trouble ignore sets of identical waveumbers if(wvn(1).eq.wvn(2))goto 200 c----- c now g(2) and wvn(2) refer to upper limit c now g(1) and wvn(1) refer to lower limit c----- c-----for the first wavenumber small use asymptotic values if(wvno.lt.1.0e-8.and.ik.eq.2)then g(4,1)=g(4,2) g(9,1)=g(9,2) endif i1 = i-1 wvno1 = wvn(i1) t01 = wvno1 * r dkk = (wvno - wvno1) * r call hank(t01,j01,j11) c----- c if the asymptotic technique is invoked, the functional c values must be changed c----- if(iasymp)then call gasym(g,aa,bb,cc,wvn,1,depth,jsrc) endif c----- c perform windowing in wavenumber domain to pass c certain ranges of phase velocity c----- call wvfilt(cmax,wvcm,wvc1,wvc2,wvcn,wvno,fact0) if(jsrc(1).eq.1)then call fmake(1,0,wvn,0,g,gg1) call wint(sumd,gg1,j01,j11,t01,0,dkk,ik) smm(1) = smm(1) + sumd endif if(jsrc(2).eq.1)then call fmake(2,0,wvn,1,g,gg1) call wint(sumd,gg1,j01,j11,t01,1,dkk,ik) smm(2) = smm(2) + sumd endif if(jsrc(3).eq.1)then call fmake(3,0,wvn,0,g,gg1) call wint(sumd,gg1,j01,j11,t01,1,dkk,ik) smm(3) = smm(3) + sumd endif sumc = cmplx(0.0,0.0) c----- c only include near field term if both SH and P-SV c computed c----- if(jsrc(4).eq.1.and.jsrc(9).eq.1)then call fmake(4,9,wvn,0,g,gg1) call wint(sumc,gg1,j01,j11,t01,1,dkk,ik) endif if(jsrc(4).eq.1)then call fmake(4,0,wvn,1,g,gg1) call wint(sumd,gg1,j01,j11,t01,0,dkk,ik) smm(4) = smm(4) + sumd - sumc/r endif if(jsrc(9).eq.1)then call fmake(9,0,wvn,1,g,gg1) call wint(sumd,gg1,j01,j11,t01,0,dkk,ik) smm(5) = smm(5) + sumd - sumc/r endif if(jsrc(5).eq.1)then call fmake(5,0,wvn,0,g,gg1) call wint(sumd,gg1,j01,j11,t01,2,dkk,ik) smm(6) = smm(6) + sumd endif sumc = cmplx(0.0,0.0) c----- c only include near field term if both SH and P-SV c computed c----- if(jsrc(6).eq.1.and.jsrc(10).eq.1)then call fmake(6,10,wvn,0,g,gg1) call wint(sumc,gg1,j01,j11,t01,2,dkk,ik) endif if(jsrc(6).eq.1)then call fmake(6,0,wvn,1,g,gg1) call wint(sumd,gg1,j01,j11,t01,1,dkk,ik) smm(7) = smm(7) + sumd - 2.*sumc/r endif if(jsrc(10).eq.1)then call fmake(10,0,wvn,1,g,gg1) call wint(sumd,gg1,j01,j11,t01,1,dkk,ik) smm(8) = smm(8) + sumd - 2.*sumc/r endif if(jsrc(7).eq.1)then call fmake(7,0,wvn,0,g,gg1) call wint(sumd,gg1,j01,j11,t01,0,dkk,ik) smm(9) = smm(9) + sumd endif if(jsrc(8).eq.1)then call fmake(8,0,wvn,1,g,gg1) call wint(sumd,gg1,j01,j11,t01,1,dkk,ik) smm(10) = smm(10) + sumd endif 201 continue t0 = t01 wvno = wvno1 200 continue c sign change due to k j(-1) smm(2) = -smm(2) smm(10) = -smm(10) do 300 i = 1,10 smm(i)=smm(i)/r ar(i) = real(smm(i)) ai(i) = aimag(smm(i)) 300 continue return end subroutine wint(smm,g1,j01,j11,t01,n,dkk,ik) common/fct/fact1 ,fact0 complex smm,g1 real j01,j11 real j21 nn = n + 1 c----- c trapezoidal rule c----- c if(ik.eq.1)then c fact = 0.5*fact0 c else c fact = fact0 c endif c----- c if commented out then rectangular rule c----- fact = fact0 if(nn.eq.1)then c----- c integral (c + d z) * j0(z) dz c----- smm = fact * (g1 * (j01 * dkk)) elseif(nn.eq.2)then c----- c integral (c + d z) j1(z) dz c----- smm = fact * (g1 * (j11 * dkk)) elseif(nn.eq.3)then if(t01.eq.0.0)then j21 = 0.0 else j21 = 2.*j11/t01 - j01 endif c----- c integral (c + d z) j2(z) dz c----- smm = fact * (g1 * (j21 * dkk)) endif return end subroutine hank(z,h0,h1) real z,h0,h1 real j1z if(z.eq.0.0)then h0 = 1.0 h1 = 0.0 elseif(z.gt.0.0 .and. z.le.3.0)then x = (z/3.)*(z/3.) h0 = 1.-x*(2.2499997-x*(1.2656208-x*(.3163866-x*( 1 .0444479-x*(.0039444-x*(.0002100)))))) j1z = 0.5-x*(.56249985-x*(.21093573-x*(.03954289-x*( 1 .00443319-x*(.00031761-x*(.00001109)))))) h1 = z * j1z else x = 3./z fac = 1./sqrt(z) f0 = .79788456+x*(-.00000077 + x*(-.00552740 + x*( 1 -.00009512+x*(.00137237+x*(-.00072805+x*(.00014476)))) 2 )) t0 = z - .78539816+x*(-.04166397+x*(-.00003954+x*( 1 .00262573+x*(-.00054125+x*(-.00029333+x*(.00013558)))) 2 )) f1 = .79788456+x*(.00000156+x*(.01659667+x*(.00017105+ 1 x*(-.00249511+x*(.00113653+x*(-.00020033)))))) t1 = z-2.35619449+x*(.12499612+x*(.00005650+x*( 1 -.00637879+x*(.00074348+x*(.00079824+x*(-.00029166))) 2 ))) h0 = fac * f0 * cos(t0) h1 = fac * f1 * cos(t1) endif return end subroutine solu(y1,y2,x1,x2,h,j,a,b,c) c we ndo not solve for a,b,c together, only two at most c thus we only need two values of wavenumber, x1 and x2 complex y1,y2,a,b,c c=cmplx(0.0,0.0) go to (300,200,300,200,300,200,300,200,200,200),j 100 continue c---------aexp(-kh) b=cmplx(0.0,0.0) a=y1*exp(x1*h) return 200 continue c---------[ a + b k ]exp(-kh) u1=x1*h u2=x2*h det=x2-x1 a= x2*(y1*exp(u1))-x1*(y2*exp(u2)) a=a/det b= y2*exp(u2) - y1*exp(u1) b=b/det return 300 continue c---------[ a + b k ] k exp(-kh) u1=x1*h u2=x2*h det=x2-x1 a = cmplx(0.0,0.0) b = x2*(y1*exp(u1))/x1 - x1*(y2*exp(u2))/x2 b = b/det c= y2*exp(u2)/x2 - y1*exp(u1)/x1 c = c/det return 400 continue c-------- a k*k*exp(-kh) a = cmplx(0.0,0.0) b = cmplx(0.0,0.0) c = y1 * exp(x1*h)/(x1)**2 return 500 continue c-------- a k exp(-kh) a = cmplx(0.0,0.0) b = y1 * exp(x1*h)/ x1 return end subroutine setup(zz,rr) c---------------------------------------------------------- c c jnkm = r integral exp(-kh) ksup m j sub n (kr) dk c c the r in front takes into account the 1/r in the c do 300 of subroutine wvint c---------------------------------------------------------- implicit double precision (a-h,o-z) real*4 rr,zz common/asym/j0k0,j0k1,j0k2,j0k3,j1k0,j1k1,j1k2,j1k3, 1 j2k0,j2k1,j2k2,j2k3 real*4 j0k0,j0k1,j0k2,j0k3,j1k0,j1k1,j1k2,j1k3,j2k0,j2k1, 1 j2k2,j2k3 r = dble(rr) z = dble(zz) dist=dsqrt(r*r + z*z) dist3=dist**3 dist5=dist**5 dist7=dist**7 rz=r*z z2=z*z r2=r*r r3=r*r2 z3=z*z2 rz2=r*z2 rz3=r*z3 zor = z/dist zor2= zor*zor zor3= zor*zor2 j0k0 = sngl( r/dist ) j0k1 = sngl( rz/dist3 ) j0k2 = sngl( (2.*rz2 - r3)/dist5 ) j0k3 = sngl( (6.*rz3 - 9.*z*r3)/dist7 ) j1k0 = sngl( 1. -z/dist ) j1k1 = sngl( r2/dist3 ) j1k2 = sngl( 3.*z*r2/dist5 ) j1k3 = sngl( 3.*r2*(4.*z2 - r2)/dist7 ) j2k0 = sngl( ( 1. -zor)*(1.-zor)*(dist/r) ) j2k1 = sngl( (1-zor)*(1-zor)*(2.+zor)/r ) j2k2 = sngl( 3.*r3/dist5 ) j2k3 = sngl( 15.*z*r3/dist7 ) return end subroutine fmake(j,k,wvn,l,g,gg1) dimension wvn(1) complex g(10,2) complex gg1 gg1= g(j,1) if(k.gt.0)then gg1=gg1 + g(k,1) endif if(l.gt.0)then gg1=gg1 * wvn(1) endif return end subroutine getgk(g,j,jsrc,wvno) c----- c read input to obtain elements of g(10,j) array c----- complex g(10,2) dimension jsrc(10) call bufrd(wvno,ierr) do 101 i=1,10 if(jsrc(i).eq.1)then call bufrd(xr,ierr) call bufrd(xi,ierr) g(i,j)=cmplx(xr,xi) endif 101 continue return end subroutine shift(g,wvn) c----- c do interchange g(j,1) -> g(j,2) c wvn(1)-> wvn(2) c----- complex g(10,2) dimension wvn(2) wvn(2)=wvn(1) do 100 i=1,10 100 g(i,2)=g(i,1) return end subroutine intini(smm,aa,bb,cc,iasymp,r) common/asym/j0k0,j0k1,j0k2,j0k3,j1k0,j1k1,j1k2,j1k3, 1 j2k0,j2k1,j2k2,j2k3 common/jout/jsrc(10) , jbdry real j0k0,j0k1,j0k2,j0k3,j1k0,j1k1,j1k2,j1k3,j2k0,j2k1, 1 j2k2,j2k3 complex aa(10),bb(10),cc(10),smm(10),sumd logical iasymp if(iasymp)then c-----set up sum arrays, but put in asymptotic value now c-----of setting to zero and then resetting smm(1)= aa(1)*j0k0 + bb(1)*j0k1 + cc(1)*j0k2 smm(2)= aa(2)*j1k1 + bb(2)*j1k2 + cc(2)*j1k3 smm(3)= aa(3)*j1k0 + bb(3)*j1k1 + cc(3)*j1k2 if(jsrc(9).eq.1 .and. jsrc(4).eq.1)then sumd = (aa(4)+aa(9))*j1k0 + (bb(4)+bb(9))*j1k1 + 1 (cc(4)+cc(9))*j1k2 else sumd = cmplx(0.0,0.0) endif sumd = -sumd/r smm(4)=sumd + aa(4)*j0k1 + bb(4)*j0k2 + cc(4)*j0k3 smm(5)=sumd + aa(9)*j0k1 + bb(9)*j0k2 + cc(9)*j0k3 smm(6)= aa(5)*j2k0 + bb(5)*j2k1 + cc(5)*j2k2 if(jsrc(6).eq.1 .and. jsrc(10).eq.1)then sumd= (aa(6)+aa(10))*j2k0 + (bb(6)+bb(10))*j2k1 + 1 (cc(6)+cc(10))*j2k2 else sumd = cmplx(0.0,0.0) endif sumd = -2.*sumd/r smm(7)=sumd + aa(6)*j1k1 + bb(6)*j1k2 + cc(6)*j1k3 smm(8)=sumd + aa(10)*j1k1+ bb(10)*j1k2+cc(10)*j1k3 smm(9)=aa(7)*j0k0 + bb(7)*j0k1 + cc(7)*j0k2 smm(10)= aa(8)*j1k1 + bb(8)*j1k2 + cc(8)*j1k3 else do 100 i=1,10 100 smm(i)=cmplx(0.0,0.0) endif return end subroutine gasym(g,aa,bb,cc,wvn,i,depth,jsrc) c----- c remove asymptotic trend from integrands c----- complex g(10,2),aa(1),bb(1),cc(1) dimension jsrc(10) real wvn(1) wvno = wvn(i) ex = exp(-wvno*depth) do 130 j=1,10 if(jsrc(j).eq.1)then g(j,i)=g(j,i) - ex*(aa(j)+wvno*(bb(j)+ 1 wvno*(cc(j)))) endif 130 continue return end subroutine wvfilt(cmax,wvcm,wvc1,wvc2,wvcn,wvno,fact) pi = 3.1415927 if(cmax.lt.0.0)then fact = 1.0 elseif(wvno.ge.wvc1.and.wvno.le.wvc2) then fact=1.0 elseif(wvno.ge.wvcm.and.wvno.lt.wvc1)then fact=(1.-cos(pi*(wvno-wvcm)/ (wvc1-wvcm)))/2. elseif(wvno.gt.wvc2.and.wvno.le.wvcn)then fact=(1.-cos(pi*(wvno-wvcn)/ (wvc2-wvcn)))/2. else fact = 0.0 endif return end EOF cat > mchdep70.f << EOF subroutine mchdep() c---------------------------------------------------------------------c c c c COMPUTER PROGRAMS IN SEISMOLOGY c c VOLUME VI c c c c PROGRAM: MCHDEP c c c c COPYRIGHT 1985 c c R. B. Herrmann c c Department of Earth and Atmospheric Sciences c c Saint Louis University c c 221 North Grand Boulevard c c St. Louis, Missouri 63103 c c U. S. A. c c c c---------------------------------------------------------------------c return end EOF cat > rhfoc10.f << EOF program rhfc10 c---------------------------------------------------------------------c c c c COMPUTER PROGRAMS IN SEISMOLOGY c c VOLUME VI c c c c PROGRAM: RHFOC10 c c c c COPYRIGHT 1985 c c R. B. Herrmann c c Department of Earth and Atmospheric Sciences c c Saint Louis University c c 221 North Grand Boulevard c c St. Louis, Missouri 63103 c c U. S. A. c c c c---------------------------------------------------------------------c parameter (LER=0, LIN=5, LOT=6) c----- c c rhfoc10 -f FILE [ -v ] [ -t -o -p -i ] -a ALP -l L [ -D -V -A ] c c This convolves output of rhwinvta with the c far-field displacement c pulses to yield velocity time c histories c The symmetric parabolic pulse ( -p ) c of Wang and Herrmann (1980) is used. c A simple triangular pulse ( -t ) is also available. c as is the Ohnaka source time function ( -o ) c The total duration of the pulse is 4L dt for the parabolic pulse c and 2L dt for the triangular pulse. c c L is the number of dt's per tau, e.g. c tau = L dt . If -n is not specified, c then tau = dt by default c c ALP is the constant in the ohnaka time function whose c far field displacement is alpha*alpha*t*exp(-alpha*t) c and whose Fourier spectrum is c [ ALP / ( j2pif + ALP ] ** 2 c This has a corner frequency of fc=ALP/2pi c c Other pulses might be able to be used. But this c one is stable. Note that for tau = dt, c -v indicates verby output on device LER c c -i = impulse source c -D = output displacement time history c -V = output velocity time history (default) c -A = output acceleration time history c c The pulses are such that the far-field dispacements would c have unit area c c The pulses used correspond to the dislocation c far-field displacement to get velocity history c c c----- parameter(NL=100) common/jout/jsrc(10),lsrc(10) common/source/z(2048) common/model/d(NL),a(NL),b(NL),rho(NL),mmax,qa(NL),qb(NL) complex data(1024),datas(513) dimension ar(10),ai(10) character*50 names character*4 icchar (11) integer idva c----- c call machine dependent initialization c----- call mchdep() call gcmdln(ntau,ipt,iverby,alp,ifile,names,idva) c----- c do some error checking, e.g., we cannot generate velocity c for triangular pulse c----- if(ifile.lt.0)write(LER,*)'no data file' if(ifile.lt.0)write(LER,*)'USAGE: ', 1 'rhfoc10 -f FILE [ -V ] [ -t -o -p -i ] -a ALP -l L [ -d -v -a]' if(ifile.lt.0)stop if(ipt.eq.2.and.alp.ge.0.0)goto 121 if(ipt.lt.2 .or. ipt.eq.3)goto 121 write(LER,*)' alp not given for ohnaka pulse' stop 121 continue if(ipt.ge.0)goto 1000 write(LER,*)' pulse shape not given' stop 1000 continue npuls=1 open(unit=3,file=names,access='sequential',form='unformatted' 1 ,status='old') open(unit=4,status='scratch',form='unformatted') rewind 3 rewind 4 19 format(5f10.5) do 99 i=1,n 99 z(i)=0.0 read(3) alpha,depth,fl,fu,dt,n,n1,n2,df,nyq2 read(3)jsrc read(3)lsrc read(3)d,a,b,rho,mmax,qa,qb nyq=nyq2/2 c default for filter if(ntau.le.0)then tau = dt else tau = ntau * dt endif if(iverby.eq.1)write(LER,*)' tau=',tau,' dt=',dt if(iverby.eq.1)write(LER,*)' ipt=',ipt,' l=',ntau if(ipt.eq.0)then call pultd(tau,dt,n,l) elseif(ipt.eq.1)then call pulpd(tau,dt,n,l) elseif(ipt.eq.2)then call pulod(alp,dt,n,l) elseif(ipt.eq.3)then call puldd(dt,n,l) endif fmax = fu inst = 0 c plot source pulse and spectra c since omega is complex time series has to be damped fac = 1.0 dfac = exp(-alpha*dt) do 200 i = 1,n data(i) = cmplx(fac*z(i),0.0) fac = fac * dfac 200 continue call four(data,n,-1,dt,df) do 310 i = 1,nyq freq=(i-1)*df if(freq.gt.fmax) go to 309 if(i.ge.n1.and.i.le.n2) go to 310 309 continue data(i)=cmplx(0.0,0.0) 310 continue do 206 i = 1,nyq 206 datas(i) = data(i) c----- c beginning of distance loop c----- 9997 continue read(3,err=9999,end=9999)r,t0 if(r.lt.0.0) go to 9999 c----- c set up output stream c----- yr = 0.0 write(LOT,110)r,yr,depth,n,t0,dt,tau 110 format(3e16.9,i10/3e16.9) do 249 jj=1,10 249 read(3,err=9999,end=9999)icchar(jj) rewind 4 do 300 i=n1,n2 do 305 j=1,10 ar(j)=0.0 ai(j)=0.0 if(jsrc(lsrc(j)).eq.1)read(3,err=9999,end=9999)ar(j),ai(j) 305 continue write(4)ar,ai 300 continue c----- c Now that all spectra for a given distance are in c make time history and output in desired order c----- do 1300 jk=1,10 do 1301 i=1,n data(i)=cmplx(0.0,0.0) 1301 continue rewind 4 do 1302 i=n1,n2 read(4)ar,ai data(i)=cmplx(ar(jk),ai(jk)) * datas(i) c----- c if integrate c----- if(idva.eq.0)then omega = 6.2831853*(i-1)*df data(i) = data(i) / cmplx(alpha,omega) elseif(idva.eq.2)then omega = 6.2831853*(i-1)*df data(i) = data(i) * cmplx(alpha,omega) endif c----- c make a real time series c----- if(i.ne.1)data(n+2-i)=conjg(data(i)) 1302 continue data(nyq)=cmplx(real(data(nyq)),0.0) call four(data,n,+1,dt,df) c----- c output the time series c----- c correct for damping fac = exp(alpha*t0) dfac = exp(alpha*dt) do 425 i = 1,n data(i)= data(i) * fac fac = fac * dfac 425 continue write(LOT,111)(real(data(i)),i=1,n) 111 format(5e16.9) 1300 continue 9000 continue go to 9997 9999 continue close(3) close(4) stop end subroutine four(data,nn,isign,dt,df) dimension data(1) n = 2 * nn if(dt.eq.0.0) dt = 1./(nn*df) if(df.eq.0.0) df = 1./(nn*dt) if(dt.ne.(nn*df)) df = 1./(nn*dt) j = 1 do 5 i=1,n,2 if(i-j)1,2,2 1 tempr = data(j) tempi = data(j+1) data(j) = data(i) data(j+1)=data(i+1) data(i) = tempr data(i+1) = tempi 2 m = n/2 3 if(j-m) 5,5,4 4 j = j-m m = m/2 if(m-2)5,3,3 5 j=j+m mmax = 2 6 if(mmax-n) 7,10,10 7 istep= 2 *mmax theta = 6.283185307/float(isign*mmax) sinth=sin(theta/2.) wstpr=-2.*sinth*sinth wstpi=sin(theta) wr=1.0 wi=0.0 do 9 m=1,mmax,2 do 8 i=m,n,istep j=i+mmax tempr=wr*data(j)-wi*data(j+1) tempi=wr*data(j+1)+wi*data(j) data(j)=data(i)-tempr data(j+1)=data(i+1)-tempi data(i)=data(i)+tempr 8 data(i+1) = data(i+1)+tempi tempr = wr wr = wr*wstpr-wi*wstpi + wr 9 wi = wi*wstpr+tempr*wstpi + wi mmax = istep go to 6 10 continue if(isign.lt.0) go to 1002 c frequency to time domain do 1001 iiii = 1,n 1001 data(iiii) = data(iiii) * df return 1002 continue c time to frequency domain do 1003 iiii = 1,n 1003 data(iiii) = data(iiii) * dt return end subroutine pulpd(tau,dt,nt,l) c far field displacement common /source/ f(2048) ltest = 0 tl = tau t1 = 0.01*dt t2 = t1 + tau t3 = t2 + tau t4 = t3 + tau t5 = t4 + tau do 100 i = 1,nt y=(i-1)*dt z = y - t1 f(i) = 0.0 if(y.gt.t1)go to 101 101 if(y.gt.t2) go to 102 f(i) = 0.5*(z/tl)**2 go to 100 102 if(y.gt.t3) go to 103 f(i)= -0.5*(z/tl)**2 + 2.*(z/tl) - 1 go to 100 103 if(y.gt.t4) go to 104 f(i)= -0.5*(z/tl)**2 + 2.*(z/tl) - 1. go to 100 104 if(y.gt.t5) go to 105 f(i)= 0.5*(z/tl)**2 - 4.*(z/tl) + 8. go to 100 105 continue ltest = ltest + 1 if(ltest.eq.1) l = i 100 continue c pulse normalized so first integral has area of unity do 200 i = 1,nt 200 f(i) = f(i)/(2.*tl) return end subroutine pulod(alp,dtt,nt,l) c far field displacement c for the ohnaka slip history c Harkrider (1976) Geophys J. 47, p 97. common/source/v(2048) ltest = 0 do 100 i=1,nt t=(i-1)*dtt v(i)=0.0 arg= alp*t al2=alp*alp v(i)=0.0 if(arg.gt.25.)goto 101 v(i)= al2*t*exp(-arg) goto 100 101 continue ltest = ltest +1 if(ltest.eq.1)l = i 100 continue return end subroutine pultd(tau,dtt,nt,l) c waveform is far-field displacement c triangular pulse c this is really a general trapezodal pulse, but only c the symmetric triangular pulse is invoked common/source/d(2048) ltest = 0 dt1=tau dt2=0.0 dt3=tau t1=dt1 t2=t1+dt2 t3=t2+dt3 fac = 0.5*dt1 + dt2 + 0.5*dt3 fac = 1./fac v1 = fac/dt1 v2 = 0.0 v3 = - fac/dt3 do 100 i=1,nt t = (i-1)*dtt d(i)=0.0 if(t.le.t1)then dt = t - 0.0 d(i) = dt*fac/dt1 elseif(t.gt.t1.and.t.le.t2)then dt = t - t1 d(i) = fac elseif(t.gt.t2.and.t.le.t3)then dt = t - t2 d(i)= fac + v3*dt elseif(t.gt.t3)then d(i)=0.0 ltest = ltest + 1 if(ltest.eq.1)l = i endif 100 continue return end subroutine puldd(dt,n,l) common/source/d(2048) c----- c Dirac Delta Pulse c----- do 100 i=1,n d(i) = 0.0 100 continue d(1) = 1.0/dt return end subroutine gcmdln(ntau,ipt,iverby,alp,ifile,names,idva) c----- c parse command line arguments c requires subroutine getarg() and function numarg() c----- character*(*) names integer*4 numarg character*50 name integer idva ntau = -1 ipt = -1 iverby = 0 alp = -1.0 ifile = -1 idva = 1 nmarg=numarg() 11 i = i + 1 if(i.gt.nmarg)goto 13 call getarg(i,name) if(name(1:2).eq.'-l')then i=i+1 call getarg(i,name) read(name,'(i10)')ntau elseif(name(1:2).eq.'-f')then ifile=1 i=i+1 call getarg(i,names) else if(name(1:2).eq.'-a')then i=i+1 call getarg(i,name) read(name,'(f10.0)')alp else if(name(1:2).eq.'-v')then iverby=1 else if(name(1:2).eq.'-t')then ipt = 0 else if(name(1:2).eq.'-p')then ipt = 1 else if(name(1:2).eq.'-o')then ipt = 2 else if(name(1:2).eq.'-i')then ipt = 3 else if(name(1:2).eq.'-D')then idva = 0 else if(name(1:2).eq.'-V')then idva = 1 else if(name(1:2).eq.'-A')then idva = 2 endif goto 11 13 continue return end EOF cat > dsplt10.f << EOF c program dsplot10 c---------------------------------------------------------------------c c c c COMPUTER PROGRAMS IN SEISMOLOGY c c VOLUME V c c c c PROGRAM: DSPLOT10 c c c c COPYRIGHT 1985 R. B. Herrmann c c c c Department of Earth and Atmospheric Sciences c c Saint Louis University c c 221 North Grand Boulevard c c St. Louis, Missouri 63103 c c U. S. A. c c c c---------------------------------------------------------------------c parameter (LER=0, LIN=5, LOT=6) dimension x(2048),y(2048) integer*4 kk(2) character*4 sym(10) character*50 name(2) data sym/' ZDD',' RDD',' ZDS',' RDS',' TDS',' ZSS',' RSS', 1 ' TSS',' ZEP',' REP'/ c----- c call machine dependent initialization c----- call mchdep() c----- c parse command line c----- call gcmdln(name,jjplt) c----- c initilize open files c----- if(jjplt.gt.0)then c----- c file name or names from command line c----- jplt = jjplt do 102 i=1,jplt kk(i) = i open(kk(i),file=name(i),status='old', 1 form='formatted',access='sequential') rewind kk(i) 102 continue else c----- c single file input from standard input c----- kk(1) = LIN jplt = 1 endif c----- c initialize plot stream c----- 100 continue do 9997 k=1,jplt read(kk(k),10,end=9999)xr,yr,hs,nt,ti,dt,tau 9997 continue 10 format(3e16.9,i10/3e16.9) write(LOT,10)xr,hs do 101 i=1,nt 101 x(i)=ti + (i-1)*dt do 200 j=1,10 c guarantee array is initialized do 9996 k=1,jplt do 300 i=1,nt 300 y(i)=0.0 read(kk(k),11,end=9999)(y(i),i=1,nt) 11 format(5e16.9) call seispl(jplt,k,x,y,nt,xr,1,sym(j),hs,tau,dt) 9996 continue 200 continue go to 100 9999 continue if(jjplt.gt.1)then do 9990 k=1,jjplt close (kk(k) ) 9990 continue endif stop end subroutine seispl(jplt,kk,x,y,n,dist,id,sym,depth,tau,dt) parameter (LER=0,LIN=5, LOT=6) character*4 sym dimension x(1),y(1) dimension xx(2050),yy(2050) data x0/0.25/,y0/8.1/,icount/1/ if(jplt.eq.2)then if(kk.eq.1)then ypt=1.3 ysc = 0.55 ysn = 0.9 elseif(kk.eq.2)then ypt=.6 ysc=0.55 ysn = 0.15 endif elseif(jplt.eq.1)then ysc=0.9 ysn = 0.25 ypt = 0.95 endif ymin = 1.0e+38 xmin = 1.0e+38 xmax =-1.0e+38 ymax = -1.0e+38 do 100 i = 1,n if(y(i).gt.ymax) ymax = y(i) if(x(i).gt.xmax) xmax = x(i) if(x(i).lt.xmin) xmin = x(i) if(y(i).lt.ymin) ymin = y(i) 100 continue ymxx=ymax if(abs(ymin).gt.ymax) ymxx = abs(ymin) if(ymax.eq.ymin.or.ymxx.eq.0.0)ymxx=1.0 c----- c output max and min values of component c----- write(LOT,8)kk,sym,ymax,ymin,ymxx 8 format(i3,a,' ymax=',e10.3,' ymin=',e10.3,' ymxx=',e10.3) do 200 i=1,n xx(i)=4.*(x(i)-xmin)/(xmax-xmin) xx(i)=xx(i)+x0 yy(i)=ysc*y(i)/ymxx + ypt yy(i)=yy(i)+y0 200 continue if(jplt.eq.2 .and. kk.eq.1)return go to (11,12,13,14,15,16,17,18,19,20),icount 11 x0 = 4.25 y0 = 8.1 go to 30 12 x0 = 0.25 y0 = 6.2 go to 30 13 x0=4.25 y0=6.2 go to 30 14 x0=0.25 y0=4.3 go to 30 15 x0=4.25 y0=4.3 go to 30 16 x0=0.25 y0=2.4 go to 30 17 x0=4.25 y0=2.4 go to 30 18 x0=0.25 y0=0.5 go to 30 19 x0=4.25 t0=0.5 go to 30 20 x0=0.25 y0=8.1 30 icount=icount+1 if(icount.gt.10)icount=icount-10 return end subroutine gcmdln(names,jplt) integer*4 jplt character names(2)*(*) character*20 name integer numarg jplt = 0 nmarg = numarg() i = 0 11 i = i + 1 if(i.gt.nmarg)goto 13 call getarg(i,name) if(name(1:3).eq.'-f1')then i = i + 1 jplt = jplt + 1 call getarg(i,names(jplt)) elseif(name(1:3).eq.'-f2')then i = i + 1 jplt = jplt + 1 call getarg(i,names(jplt)) else call usage() endif goto 11 13 continue return end subroutine usage integer*4 LER, LIN, LOT parameter(LER=0, LIN=5, LOT=6) write(LER,*)'USAGE: dsplot10 [-f1 fname1] [-f2 fname2]' stop end EOF cat > TEST00 << EOF1 date hostname make all cat > hspec8.dat < rhwvinta.dat < hspec9.out < file10 ./dsplt10 < file10 rm hspec8.dat rhwvinta.dat hspec9.out file10 EOF1 chmod +x TEST00 TEST00 rm TEST00 Makefile *.f hspec9 rhfoc10 dsplt10 rm swe0010