c****************************************************************** c Modified version of Ray Buland's SETBRN program which c runs on both PC's and SUN's c c Note: On PC's, SETBRN is too big to run under MSDOS c to run it under windows using MS Fortran 5.1, compile and c link it as follows (I assume that you installed your compiler c in the \FORTRAN subdirectory) c c fl -c -Gt -AH -MW setbrn.for c fl -c -Gt -AH -MW tau.for c link setbrn+tau,,,/NOD:llibfor7 llibfew,\fortran\binb\fl.def; c c Don't try using the programmer's workbench (PWB): its NMAKE facility c is all screwed up when it comes to making Window's executables c c 3/94 Recompiled using MS Fortran Powerstation: (32 bit compiler which c generates executables that run under MSDOS on a 386 or better). c c Should also work with other 32 bit compilers (Watcom, Microway, etc) c c jh mar 1 99 ------------- Version 7.0 check ---------------- c no changes c c program setbrn c c Herewith the new version of setbrn.f with separated table and header c and correct index assignments for direct access (hopefully) c include 'limits.inc' save character*8 code,phcd character*20 modnam character*80 fnam double precision zm,pm,pb,pu,taup,xp,taul,px,xt,xl,pux,pt,taut, 1 coef,xa double precision tmp(nsl1,2),xm(nsl1,2),deg,dtol,zmax,zoc,zic, 1 z0 dimension ndx2(nsr0,2) common/umodc/zm(nsr0,2),pm(nsr0,2),ndex(nsr0,2),mt(2) common/brkptb/lcb(2),lbb(2),lbrk(nbr2,2) common/brkptc/code(nbr1,2) common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) common/segc/fcs(jseg,3),nafl(jseg,3),indx(jseg,2),kndx(jseg,2), 1 nseg common/brnb/px(jbrn,2),xt(jbrn,2),jndx(jbrn,2),mndx(jbrn,2),nbrn common/brnc/phcd(jbrn) common/msqc/pux(jbrn,2),km(2),midx(jbrn,2) common/outc/pt(jout),taut(jout),coef(5,jout),xa(jout),nl data nin,nout,xmin,dtol/1,2,200.,1d-6/ c c print version c include 'version.inc' out_version_date='July 23, 2001' if (version_new) out_version_date=version_date call print_ver deg=180d0/3.1415927d0 c c write(6,*) "rec length for dasign:" c read(5,*) ndasr c fnam='remodl.hed'//char(0) call assign(nin,-1,fnam) read(nin)ndasr,modnam,zmax,zoc,zic,kb,(pb(i),i=1,kb(2)), 1 mt,lt,lbb,lcb,xn,pn,tn read(nin)((lbrk(i,nph),i=1,lbb(nph)),(code(i,nph), 1 i=1,lcb(nph)),(zm(i,nph),pm(i,nph),ndex(i,nph),i=1,mt(nph)), 2 (lvz(i,nph),taul(i,nph),xl(i,nph),i=1,lt(nph)),nph=1,2) call iasp91_filename(modnam) call retrns(nin) print *,'ndasr =',ndasr,' modnam = ',modnam c call dasign(nin,-1,'remodl.tbl'//char(0),ndasr) nrec=0 do 1 nph=1,2 n1=kb(nph) ind=0 do 2 k=1,n1 2 xm(k,nph)=0d0 3 nrec=nrec+1 read(nin,rec=nrec)z0,n,(tmp(k,1),k=1,n),(tmp(k,2),k=1,n) if(ind.gt.0) go to 4 if(dabs(z0-zoc).le.dtol) go to 4 j=1 do 5 k=2,n xm(k,nph)=dmax1(xm(k,nph),dabs(tmp(j,2)-tmp(k,2))) 5 j=k if(n+1.eq.n1) xm(n1,nph)=tmp(n,2) go to 3 4 ind=ind+1 do 6 k=1,n taup(k,ind,nph)=tmp(k,1) 6 xp(k,ind,nph)=tmp(k,2) if(ind.lt.3) go to 3 1 continue c xmin=xn*xmin c fnam='setbrn1.lis'//char(0) call assign(10,2,fnam) write(10,*)'kb mt lt lbb lcb',kb,mt,lt,lbb,lcb write(10,*)'xn pn tn xmin',xn,pn,tn,xmin cn=1./xn write(10,209)(i,(lbrk(i,j),code(i,j),j=1,2),i=1,lbb(1)) 209 format(/(1x,2i5,2x,a,i5,2x,a)) write(10,210)(i,lbrk(i,2),code(i,2),i=lbb(1)+1,lbb(2)) 210 format(1x,i5,15x,i5,2x,a) write(10,200,iostat=ios)(i,(zm(i,j),pm(i,j),ndex(i,j),j=1,2), 1 i=1,mt(1)) 200 format(/(1x,i5,2f12.6,i5,2x,2f12.6,i5)) write(10,201,iostat=ios)(i,zm(i,2),pm(i,2),ndex(i,2), 1 i=mt(1)+1,mt(2)) 201 format(1x,i5,31x,2f12.6,i5) write(10,217)((nph,i,lvz(i,nph),taul(i,nph),deg*xl(i,nph), 1 i=1,lt(nph)),nph=1,2) 217 format(/(1x,3i5,f12.6,f12.2)) write(10,202)(i,pb(i),cn*xm(i,1),cn*xm(i,2),i=1,kb(1)) 202 format(/(5x,i5,f12.6,2f12.2)) write(10,203)(i,pb(i),cn*xm(i,2),i=kb(1)+1,kb(2)) 203 format(5x,i5,f12.6,12x,f12.2) call retrns(10) fnam='setbrn2.lis'//char(0) call assign(10,2,fnam) c do 8 nph=1,2 n1=kb(nph) do 9 i=2,n1 xm(i,nph)=xm(i-1,nph)+xm(i,nph) pu(i,nph)=pb(i) 9 kuse(i,nph)=-1 do 8 ind=3,2,-1 jnd=ind-1 do 8 i=1,n1 taup(i,ind,nph)=taup(i,ind,nph)-taup(i,jnd,nph) 8 xp(i,ind,nph)=xp(i,ind,nph)-xp(i,jnd,nph) do 10 nph=1,2 10 call pdecx(kb(nph),nph,xm,2.) c write(10,*)'ku',ku write(10,204)(i,(pu(i,nph),cn*xm(i,nph),cn*(xm(i+1,nph)- 1 xm(i,nph)),nph=1,2),i=1,ku(1)) 204 format(/(5x,i5,2(f12.6,2f12.2))) write(10,205)(i,pu(i,2),cn*xm(i,2),cn*(xm(i+1,2)-xm(i,2)), 1 i=ku(1)+1,ku(2)) 205 format(5x,i5,36x,f12.6,2f12.2) do 207 nph=1,2 207 write(10,206)(i,pb(i),(taup(i,j,nph),j=1,3),(deg*xp(i,j,nph), 1 j=1,3),i=1,kb(nph)) 206 format(/(1x,i5,4f10.6,3f10.2)) c call layout c write(10,214)(i,pb(i),(kuse(i,j),j=1,2),i=1,kb(2)) c214 format(/(5x,i5,f12.6,2i5)) do 11 nph=1,2 n1=kb(nph) k=0 do 12 i=1,n1 if(kuse(i,nph).lt.0) go to 12 k=k+1 pu(k,nph)=pb(i) 12 continue 11 ku(nph)=k call kseq call mseq write(10,215)(i,(pu(i,j),j=1,2),i=1,ku(1)) 215 format(/(5x,i5,2f12.6)) write(10,216)(i,pu(i,2),i=ku(1)+1,ku(2)) 216 format(5x,i5,12x,f12.6) write(10,208)(i,(nafl(i,j),j=1,3),(indx(i,j),j=1,2),(kndx(i,j), 1 j=1,2),(fcs(i,j),j=1,3),i=1,nseg) 208 format(/(1x,8i6,3f6.1)) write(10,211)(i,(jndx(i,j),j=1,2),(mndx(i,j),j=1,2),(px(i,j), 1 j=1,2),(deg*xt(i,j),j=1,2),phcd(i),i=1,nbrn) 211 format(/(1x,i3,4i5,2f12.6,2f10.2,2x,a)) write(10,218)(i,(midx(i,j),j=1,2),(pux(i,j),j=1,2), 1 i=1,max0(km(1),km(2))) 218 format(/(1x,i3,2i5,2f12.6)) write(10,212,iostat=ios)(i,pt(i),taut(i),deg*xa(i), 1 cn*(xa(i)-xa(i+1)),(coef(j,i),j=1,5),i=1,nl) 212 format(/(1x,i4,0p2f12.6,2f10.2,1p5d10.2)) call retrns(10) c fnam='setbrn3.lis'//char(0) call assign(10,2,fnam) do 20 nph=1,2 mt(nph)=mt(nph)-3 ku(nph)=ku(nph)-1 20 km(nph)=km(nph)-1 c icor=33 - originally 32 records used as header in setbrn c and 2 records used as header in remodl. icor=3 do 14 nph=1,2 m1=mt(nph) icor=icor-3 do 14 i=2,m1 14 ndx2(i,nph)=ndex(i,nph)+icor len1=ku(2)+km(2) len0=8*len1 len2=5*nl write(10,*)'nseg nbrn mt ku km len len1',nseg,nbrn,mt,ku,km,len0, 1 len1 write(10,*) nasgr = len0 write(6,*) 'reclength for direct access', nasgr c++ c write(6,*) 'enter model name' c read(5,*) cmodel nb=index(modnam,' ')-1 if(nb.le.0) nb=len(modnam) c cnam1 = cmodel(1:nb)//'.tbl' c cnam2 = cmodel(1:nb)//'.hed' write(6,*) 'header file :',modnam(1:nb)//'.HED' write(6,*) 'table file :',modnam(1:nb)//'.TBL' fnam=modnam(1:nb)//'.HED'//char(0) call assign(nout,-2,fnam) c++ write(nout) nasgr,nl,len2,xn,pn,tn,mt,nseg,nbrn,ku,km,fcs do i=1,jseg write(nout)(nafl(i,j),j=1,3) write(nout)(indx(i,j),j=1,2) write(nout)(kndx(i,j),j=1,2) enddo do i=1,nsr0 write(nout)(pm(i,j),j=1,2) write(nout)(zm(i,j),j=1,2) end do do i=1,nsr0 write(nout)(ndx2(i,j),j=1,2) enddo do i=1,nsl1 write(nout)(pu(i,j),j=1,2) end do do i=1,jbrn write(nout)(pux(i,j),j=1,2) end do write(nout) phcd do i=1,jbrn write(nout)(px(i,j),j=1,2) write(nout)(xt(i,j),j=1,2) write(nout)(jndx(i,j),j=1,2) end do write(nout) pt,taut do i=1,jout write(nout)(coef(j,i),j=1,5) c write(*,'(5f12.4)')(coef(j,i),j=1,5) end do call retrns(nout) c call dasign(nout,-2,modnam(1:nb)//'.TBL'//char(0),nasgr) nrec = 0 do 16 nph=1,2 m1=mt(nph) n1=ku(nph) k1=km(nph) write(10,*)'nph m1 n1 k1',nph,m1,n1,k1 do 16 m=2,m1 if(ndex(m,nph).eq.ndex(m-1,nph)) go to 16 read(nin,rec=ndex(m,nph))z0,n,(tmp(k,1),k=1,n),(tmp(k,2),k=1,n) write(10,*)'m nph ndex n',m,nph,ndex(m,nph),n k=0 l=1 do 17 i=1,n if(kuse(i,nph).lt.0) go to 17 if(dabs(pux(l,nph)-pb(i)).gt.dtol) go to 18 tmp(l,2)=tmp(i,2) l=l+1 18 k=k+1 tmp(k,1)=tmp(i,1) 17 continue write(10,*)'k l nrec',k,l-1,nrec+1,ndx2(m,nph),sngl(tmp(1,1)) if(k.ge.n1) go to 19 k=k+1 do 21 i=k,n1 21 tmp(i,1)=0d0 19 if(l.gt.k1) go to 23 do 22 i=l,k1 22 tmp(i,2)=0d0 23 nrec=nrec+1 write(nout,rec=nrec)(tmp(i,1),i=1,n1),(tmp(i,2),i=1,k1) 16 continue c call retrns(10) call retrns(nin) call retrns(nout) call vexit(0) end c subroutine pdecx(n1,nph,xm,fac) save include 'limits.inc' double precision pb,pu,taup,xp,taul,xl double precision xm(nsl1,2),ptol,pa,pax,plim common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) data ptol/.03/ c call collct(1,n1,xm(1,nph),fac*xmin) k=0 plim=.7d0*pu(n1,nph) do 1 i=1,n1 if(xm(i,nph).lt.0d0) go to 1 if(pu(i,nph).lt.plim) go to 2 if(pu(i,nph)-pu(k,nph).le.ptol) go to 2 pa=pu(k,nph)+.75d0*(pu(i,nph)-pu(k,nph)) pax=1d10 m=0 do 3 j=i1,i if(dabs(pu(j,nph)-pa).ge.pax) go to 3 m=j pax=dabs(pu(j,nph)-pa) 3 continue if(m.eq.i1.or.m.eq.i) go to 2 k=k+1 pu(k,nph)=pu(m,nph) xm(k,nph)=0d0 kuse(m,nph)=1 2 k=k+1 i1=i pu(k,nph)=pu(i,nph) xm(k,nph)=xm(i,nph) kuse(i,nph)=1 1 continue ku(nph)=k return end c subroutine collct(i1,i2,x,xmn) c c $$$$$ calls varn $$$$$ c save double precision x(i2) c is=i1+1 ie=i2-1 if(ie.lt.is) return k1=i1 var=0. m=0 do 1 i=is,ie dx1=dabs(x(k1)-x(i))-xmn dx2=dabs(x(k1)-x(i+1))-xmn if(abs(dx2).ge.abs(dx1)) go to 2 x(i)=-x(i) go to 1 2 if(k1.le.i1) kb=i k1=i var=var+dx1*dx1 m=m+1 1 continue dx1=dabs(x(k1)-x(i2))-xmn var=var+dx1*dx1 m=m+1 7 if(m.le.1) return k1=i1 k2=kb ks=kb+1 nch=0 write(10,*)'ks,i2',ks,i2 do 8 i=ks,i2 if(x(i).lt.0d0) go to 8 k0=k1 k1=k2 k2=i var1=varn(x,k0,k1,k2,k1-1,xmn,var,m,m1) var2=varn(x,k0,k1,k2,k1+1,xmn,var,m,m2) if(amin1(var1/m1,var2/m2).ge.var/m) go to 6 nch=nch+1 x(k1)=-x(k1) if(var1/m1-var2/m2)3,4,5 4 if(m1-m2)3,3,5 3 k1=k1-1 x(k1)=dabs(x(k1)) var=var1 m=m1 go to 6 5 k1=k1+1 x(k1)=dabs(x(k1)) var=var2 m=m2 6 if(k0.eq.i1) kb=k1 8 continue if(nch.gt.0) go to 7 return end c function varn(x,k0,k1,k2,kt,xmn,var,m,mn) c c $$$$$ calls only library routines $$$$$ c save double precision x(k2) c dx1=dabs(x(k0)-x(k1))-xmn dx2=dabs(x(k1)-x(k2))-xmn varn=var-dx1*dx1-dx2*dx2 if(kt.le.k0.or.kt.ge.k2) go to 1 dx1=dabs(x(k0)-dabs(x(kt)))-xmn dx2=dabs(dabs(x(kt))-x(k2))-xmn varn=varn+dx1*dx1+dx2*dx2 mn=m return 1 dx1=dabs(x(k0)-dabs(x(k2)))-xmn varn=varn+dx1*dx1 mn=m-1 return end c subroutine layout c c Layout contains the program for the desired travel-time segments c implemented as calls to the mk_br entry points. Each call does one c segment (which may have many branches). c save include 'limits.inc' double precision pb,pu,taup,xp,taul,xl,px,xt,pt,taut,coef,xa double precision dir(3),cref(3),sref(3) common/brkptb/lcb(2),lbb(2),lbrk(nbr2,2) common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) common/segc/fcs(jseg,3),nafl(jseg,3),indx(jseg,2),kndx(jseg,2), 1 nseg common/brnb/px(jbrn,2),xt(jbrn,2),jndx(jbrn,2),mndx(jbrn,2),nbrn common/outc/pt(jout),taut(jout),coef(5,jout),xa(jout),nl data dir,cref,sref/1d0,1d0,1d0,1d0,2d0,2d0,2d0,2d0,2d0/ c c Initialize variables. nseg=0 nbrn=0 nl=0 do 1 j=1,3 do 1 i=1,jseg 1 fcs(i,j)=0. do 2 i=1,jout 2 taut(i)=0d0 do 3 j=1,2 do 3 i=1,jbrn 3 xt(i,j)=0d0 c c Do all of the segments. c c P (up-going branch). print *,'Layout: do Pup' call mkubr(ku(1), +1) c P, Pdiff, and PKP. print *,'Layout: do P and PKP' call mkdbr(1,lbb(1),-1,3,1,1,dir) c PKiKP. print *,'Layout: do PKiKP' call mkrbr(2, -1,2,1,1,dir) c pP. print *,'Layout: do pP' call mkdbr(1,lbb(1),+1,3,1,1,dir) c sP. print *,'Layout: do sP' call mkdbr(1,lbb(1),+2,3,1,1,dir) c pPKiKP. print *,'Layout: do pPKiKP' call mkrbr(2, +1,2,1,1,dir) c sPKiKP. print *,'Layout: do sPKiKP' call mkrbr(2, +2,2,1,1,dir) c PcP. print *,'Layout: do PcP' call mkrbr(3, -1,1,1,1,dir) c ScP. print *,'Layout: do ScP' call mkrbr(3, -2,1,2,1,dir) c SKP. print *,'Layout: do SKP' call mkdbr(1,3, -2,3,2,1,dir) c SKiKP. print *,'Layout: do SKiKP' call mkrbr(2, -2,2,2,1,dir) c PKKP. print *,'Layout: do PKKP' call mkdbr(1,3, -1,3,1,1,cref) c SKKP. print *,'Layout: do SKKP' call mkdbr(1,3, -2,3,2,1,cref) c PP and P'P'. print *,'Layout: do PP, P''P''' call mkdbr(1,lbb(1),-1,3,1,1,sref) c S (up-going branch). print *,'Layout: do Sup' call mkubr(ku(2), +2) c S, Sdiff, and SKS. print *,'Layout: do S and SKS' call mkdbr(1,lbb(2),-2,3,2,2,dir) c pS print *,'Layout: do pS' call mkdbr(1,lbb(1),+1,3,2,2,dir) c sS print *,'Layout: do sS' call mkdbr(1,lbb(2),+2,3,2,2,dir) c ScS print *,'Layout: do ScS' call mkrbr(4, -2,1,2,2,dir) c PcS print *,'Layout: do PcS' call mkrbr(3, -1,1,1,2,dir) c PKS print *,'Layout: do PKS' call mkdbr(1,3, -1,3,1,2,dir) c PKKS print *,'Layout: do PKKS' call mkdbr(1,3, -1,3,1,2,cref) c SKKS print *,'Layout: do SKKS' call mkdbr(1,3, -2,3,2,2,cref) c SS and S'S'. print *,'Layout: do SS and S''S''' call mkdbr(1,lbb(2),-2,3,2,2,sref) c SP print *,'Layout: do SP' call mkcbr(4,lbb(1),-2,1,2,1,sref) c PS print *,'Layout: do PS' call mkcbr(4,lbb(1),-1,1,1,2,sref) return end c subroutine mkdbr(l1,l2,isgn,lyr,nph,kph,fac) c c Mkdbr sets up a simple refracted wave segment. L1 and l2 point to the c lbrk array of slowness break point pointers. Note that the P and S c break point arrays don't necessarily line up layer by layer. This is c not generally a problem as most phases need only worry about the c pointer to the surface slowness for one wave type and a pointer c somewhere in the core (which is constrained to be the same for both c P and S). Isgn is positive if the wave starts out going up and c negative if the wave starts out going down. Iabs(isng) is 1 if the c wave starts out as a P wave and 2 if the wave starts out as an S wave. c Lyr gives the number of major layers (mantle, outer core, and inner c core) that the wave penetrates. Nph and kph give the wave type (1 for c P and 2 for S) on the down-going and up-going legs of the ray path c respectively. Fac is a three element array giving the number of c repeats of the ray path in each major layer. This scheme incorporates c turning rays (e.g., P and S), turning rays reflected, but not c converted at the surface (e.g., PP and SS), up-going rays reflected c and/or converted at the surface into turning rays (e.g., pP and sP), c turning rays converted during transmission through an interface (e.g., c SKP and PKS), and rays which turn multiple times while reflecting from c the bottom side of a layer (e.g., PKKP or SKKP). Mkdbr does not c include up-going rays (to the receiver), rays reflected from the top c side of a discontinuity, or rays which are reflected and converted at c the free surface. See mkubr, mkrbr, and mkcbr respectively for c routines which handle these types of rays. c save include 'limits.inc' character*8 code,phcd,ks,phcd1,phcd2 double precision pb,pu,taup,xp,taul,xl,px,xt,pt,taut,coef,xa double precision fac(3) common/brkptb/lcb(2),lbb(2),lbrk(nbr2,2) common/brkptc/code(nbr1,2) common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) common/segc/fcs(jseg,3),nafl(jseg,3),indx(jseg,2),kndx(jseg,2), 1 nseg common/brnb/px(jbrn,2),xt(jbrn,2),jndx(jbrn,2),mndx(jbrn,2),nbrn common/brnc/phcd(jbrn) common/outc/pt(jout),taut(jout),coef(5,jout),xa(jout),nl data ks/'KKKKKKKK'/ c c Remember the programming as part of the final phase construction is c done in depcor. nseg=nseg+1 nafl(nseg,1)=isgn nafl(nseg,2)=nph nafl(nseg,3)=kph indx(nseg,1)=nl+1 kndx(nseg,1)=1 c Using l1 and l2 to get the breakpoints has some shortcommings, c particularly for converted phases. It would be more general to c have separate indicies for the breakpoints and the layers covered. if(l1.gt.1) kndx(nseg,1)=lbrk(l1-1,nph) kndx(nseg,2)=min0(min0(lbrk(l2,nph),lbrk(l2,kph)), 1 lbrk(l2,iabs(isgn))) print *,'Mkdbr: l1 l2 isgn lyr nph kph =',l1,l2,isgn,lyr,nph,kph print *,'Mkdbr: nseg kndx indx =',nseg,kndx(nseg,1),kndx(nseg,2), 1 indx(nseg,1) xfc=0. do 14 m=1,lyr fcs(nseg,m)=fac(m) 14 xfc=amax1(xfc,fcs(nseg,m)) c c Set up the required slownesses, taus and distances. c j=kndx(nseg,1) lz1=1 lz2=1 c Loop over the layers of interest. do 1 i=l1,l2 c Be sure that the phase cuts off at the right place. l=min0(lbrk(i,nph),kndx(nseg,2)) c Skip all total internal reflections. if(code(i,nph)(1:1).eq.'r'.or.j.ge.l) go to 1 c Set the starting branch pointer. nbrn=nbrn+1 nt=nl+1 jndx(nbrn,1)=nt c Copy in the desired slownesses. do 2 k=j,l nl=nl+1 pt(nl)=pb(k) c Add up the tau contributions. do 2 m=1,lyr 2 taut(nl)=taut(nl)+fac(m)*(taup(k,m,nph)+taup(k,m,kph)) c Take care of branch end pointers and slownesses. mndx(nbrn,1)=j mndx(nbrn,2)=l px(nbrn,1)=pb(j) px(nbrn,2)=pb(l) c Add up distance contributions for the branch end points only. do 3 m=1,lyr xt(nbrn,1)=xt(nbrn,1)+fac(m)*(xp(j,m,nph)+xp(j,m,kph)) 3 xt(nbrn,2)=xt(nbrn,2)+fac(m)*(xp(l,m,nph)+xp(l,m,kph)) c Take care of the contribution due to low velocity zones for the c down-going leg(s). if(j.ne.lvz(lz1,nph)) go to 9 do 11 m=1,lyr taut(nt)=taut(nt)-fac(m)*taup(j,m,nph) 11 xt(nbrn,1)=xt(nbrn,1)-fac(m)*xp(j,m,nph) taut(nt)=taut(nt)+fac(1)*taul(lz1,nph) xt(nbrn,1)=xt(nbrn,1)+fac(1)*xl(lz1,nph) lz1=lz1+1 c Take care of the contributions due to low velocity zones for the c up-going leg(s). 9 if(j.ne.lvz(lz2,kph)) go to 10 do 12 m=1,lyr taut(nt)=taut(nt)-fac(m)*taup(j,m,kph) 12 xt(nbrn,1)=xt(nbrn,1)-fac(m)*xp(j,m,kph) taut(nt)=taut(nt)+fac(1)*taul(lz2,kph) xt(nbrn,1)=xt(nbrn,1)+fac(1)*xl(lz2,kph) lz2=lz2+1 c Decimate the slownesses if the branch is oversampled in distance. 10 call pdect(jndx(nbrn,1),nl,j,iabs(isgn),xfc) c Set up the interpolation. call tauspl(jndx(nbrn,1),nl,pt,coef) c Remember the final branch end slowness value. jndx(nbrn,2)=nl c c Take care of the branch name. First, set up a default. phcd1=code(i,kph)(3:) phcd(nbrn)=code(i,nph)(2:2)//phcd1(1:7) if(idint(fac(1)+.5d0).gt.1) go to 5 if(idint(fac(2)+.5d0).le.1) go to 4 c Re-do the name if the ray is reflected from the underside of the c core-mantle boundary. ind=idint(fac(2)-.5d0) phcd1=code(i,kph)(3:) phcd(nbrn)=code(i,nph)(2:2)//ks(1:ind)//phcd1(1:7-ind) go to 4 c Re-do the name if the ray is reflected from the surface. 5 phcd1=code(i,kph)(2:) if(code(i,nph)(3:3).eq.' ') phcd(nbrn)=code(i,nph)(2:2)// 1 phcd1(1:7) phcd1=code(i,kph)(2:) if(code(i,nph)(3:3).ne.' '.and.code(i,nph)(3:3).ne.'K') 1 phcd(nbrn)=code(i,nph)(2:3)//phcd1(1:6) phcd1=code(i,kph)(2:2) phcd2=code(i,kph)(5:) if(code(i,nph)(3:3).eq.'K') phcd(nbrn)=code(i,nph)(2:2)//''''// 1 phcd1(1:1)//''''//phcd2(1:4) c Take care . 4 ind=max0(index(phcd(nbrn),'KSab'),index(phcd(nbrn),'S''ab')) if(phcd(nbrn)(1:1).eq.'S'.and.ind.gt.0) phcd(nbrn)(ind+2:ind+3)= 1 'ac' phcd1=phcd(nbrn) if(isgn.eq.1) phcd(nbrn)='p'//phcd1(1:7) if(isgn.eq.2) phcd(nbrn)='s'//phcd1(1:7) 1 j=l indx(nseg,2)=nl return c c Mkubr handles up-going P and S. L1 and isgn are as for mkdbr (except c that l1 actually plays the role of l2 with the beginning break point c assumed to be zero). The other arguments are not needed. c entry mkubr(l1,isgn) nseg=nseg+1 nafl(nseg,1)=isgn nafl(nseg,2)=0 nafl(nseg,3)=0 indx(nseg,1)=nl+1 kndx(nseg,1)=1 l=kb(iabs(isgn)) kndx(nseg,2)=l print *,'Mkubr: l1 isgn =',l1,isgn print *,'Mkubr: nseg kndx indx =',nseg,kndx(nseg,1),kndx(nseg,2), 1 indx(nseg,1) nbrn=nbrn+1 jndx(nbrn,1)=nl+1 do 6 k=1,l1 nl=nl+1 pt(nl)=pu(k,isgn) 6 xa(nl)=0d0 mndx(nbrn,1)=1 mndx(nbrn,2)=l px(nbrn,1)=pb(1) px(nbrn,2)=pb(l) call tauspl(jndx(nbrn,1),nl,pt,coef) jndx(nbrn,2)=nl phcd(nbrn)=code(1,iabs(isgn))(2:2) indx(nseg,2)=nl return c c Mkrbr handles reflected phases possibly with a conversion such as c PcP, PcS, and PkiKP. Arguments are as for mkdbr (except that l1 c actually plays the role of l2 with the beginning break point assumed c to be zero). c entry mkrbr(l1,isgn,lyr,nph,kph,fac) nseg=nseg+1 nafl(nseg,1)=isgn nafl(nseg,2)=nph nafl(nseg,3)=kph indx(nseg,1)=nl+1 kndx(nseg,1)=1 l=min0(lbrk(l1,nph),lbrk(l1,kph)) kndx(nseg,2)=l print *,'Mkrbr: l1 isgn lyr nph kph =',l1,isgn,lyr,nph,kph print *,'Mkrbr: nseg kndx indx =',nseg,kndx(nseg,1),kndx(nseg,2), 1 indx(nseg,1) xfc=0. do 15 m=1,lyr fcs(nseg,m)=fac(m) 15 xfc=amax1(xfc,fcs(nseg,m)) if(lyr.ge.2) xfc=2. c nbrn=nbrn+1 jndx(nbrn,1)=nl+1 do 7 k=1,l nl=nl+1 pt(nl)=pb(k) do 7 m=1,lyr 7 taut(nl)=taut(nl)+fac(m)*(taup(k,m,nph)+taup(k,m,kph)) mndx(nbrn,1)=1 mndx(nbrn,2)=l px(nbrn,1)=pb(1) px(nbrn,2)=pb(l) do 8 m=1,lyr 8 xt(nbrn,2)=xt(nbrn,2)+fac(m)*(xp(l,m,nph)+xp(l,m,kph)) c write(10,*)'jndx1',jndx(nbrn,1) call pdect(jndx(nbrn,1),nl,1,iabs(isgn),xfc) call tauspl(jndx(nbrn,1),nl,pt,coef) c write(10,*)'jndx2',jndx(nbrn,1) jndx(nbrn,2)=nl c changed BRL phcd1=code(l1,kph)(2:2) if(lyr.eq.1) phcd(nbrn)=code(l1,nph)(2:2)//'c'//phcd1(1:6) phcd1=code(l1,kph)(3:) if(lyr.eq.2) phcd(nbrn)=code(l1,nph)(2:2)//phcd1(1:7) phcd1=phcd(nbrn)(1:7) if(isgn.eq.1) phcd(nbrn)='p'//phcd1(1:7) if(isgn.eq.2) phcd(nbrn)='s'//phcd1(1:7) indx(nseg,2)=nl return c c Mkcbr handles phases reflected and converted at the surface such as c PS and SP. Arguments are as for mkdbr. c entry mkcbr(l1,l2,isgn,lyr,nph,kph,fac) if(nph.gt.0.and.kph.gt.0.and.nph.ne.kph) go to 29 print *,'Mkcbr: bad call - nph kph =',nph,kph call vexit(1) 29 nseg=nseg+1 nafl(nseg,1)=isgn nafl(nseg,2)=nph nafl(nseg,3)=kph indx(nseg,1)=nl+1 kndx(nseg,1)=1 if(l1.gt.1) kndx(nseg,1)=min0(lbrk(l1,nph),lbrk(l1,kph)) kndx(nseg,2)=min0(min0(lbrk(l2,nph),lbrk(l2,kph)), 1 lbrk(l2,iabs(isgn))) print *,'Mkcbr: l1 l2 isgn lyr nph kph =',l1,l2,isgn,lyr,nph,kph print *,'Mkcbr: nseg kndx indx =',nseg,kndx(nseg,1),kndx(nseg,2), 1 indx(nseg,1) xfc=0. do 16 m=1,lyr fcs(nseg,m)=fac(m) 16 xfc=amax1(xfc,fcs(nseg,m)) c j=kndx(nseg,1) lz1=1 lz2=1 ik=l1 c print *,'Mkcbr: start loop' do 17 in=l1,l2 31 l=min0(lbrk(in,nph),kndx(nseg,2)) if(code(in,nph)(1:1).eq.'r'.or.j.ge.l) go to 17 l=min0(lbrk(ik,kph),kndx(nseg,2)) if(code(ik,kph)(1:1).ne.'r'.and.j.lt.l.or.ik.ge.l2) go to 28 j=max0(j,l) ik=ik+1 go to 31 c 28 if(lbrk(in,nph).le.lbrk(ik,kph)) go to 26 l=min0(lbrk(ik,kph),kndx(nseg,2)) print *,'kph: kph ik j l code =',kph,ik,j,l,' ',code(ik,kph) isw=2 go to 27 26 l=min0(lbrk(in,nph),kndx(nseg,2)) print *,'nph: nph in j l code =',nph,in,j,l,' ',code(in,nph) isw=1 c 27 nbrn=nbrn+1 nt=nl+1 jndx(nbrn,1)=nt do 18 k=j,l nl=nl+1 pt(nl)=pb(k) do 18 m=1,lyr 18 taut(nl)=taut(nl)+fac(m)*(taup(k,m,nph)+taup(k,m,kph)) mndx(nbrn,1)=j mndx(nbrn,2)=l px(nbrn,1)=pb(j) px(nbrn,2)=pb(l) do 19 m=1,lyr xt(nbrn,1)=xt(nbrn,1)+fac(m)*(xp(j,m,nph)+xp(j,m,kph)) 19 xt(nbrn,2)=xt(nbrn,2)+fac(m)*(xp(l,m,nph)+xp(l,m,kph)) if(j.ne.lvz(lz1,nph)) go to 20 do 21 m=1,lyr taut(nt)=taut(nt)-fac(m)*taup(j,m,nph) 21 xt(nbrn,1)=xt(nbrn,1)-fac(m)*xp(j,m,nph) taut(nt)=taut(nt)+fac(1)*taul(lz1,nph) xt(nbrn,1)=xt(nbrn,1)+fac(1)*xl(lz1,nph) lz1=lz1+1 20 if(j.ne.lvz(lz2,kph)) go to 22 do 23 m=1,lyr taut(nt)=taut(nt)-fac(m)*taup(j,m,kph) 23 xt(nbrn,1)=xt(nbrn,1)-fac(m)*xp(j,m,kph) taut(nt)=taut(nt)+fac(1)*taul(lz2,kph) xt(nbrn,1)=xt(nbrn,1)+fac(1)*xl(lz2,kph) lz2=lz2+1 22 call pdect(jndx(nbrn,1),nl,j,iabs(isgn),xfc) call tauspl(jndx(nbrn,1),nl,pt,coef) jndx(nbrn,2)=nl c c changed BRL phcd1=code(ik,kph)(2:) if(code(in,nph)(3:3).eq.' ') phcd(nbrn)=code(in,nph)(2:2)// 1 phcd1(1:7) if(code(in,nph)(3:3).ne.' '.and.code(in,nph)(3:3).ne.'K') 1 phcd(nbrn)=code(in,nph)(2:3)//phcd1(1:6) phcd1=code(ik,kph)(5:) if(code(in,nph)(3:3).eq.'K') phcd(nbrn)=code(in,nph)(2:2)//''''// 1 code(ik,kph)(2:2)//''''//phcd1(1:4) phcd1=phcd(nbrn) if(isgn.eq.1) phcd(nbrn)='p'//phcd1(1:7) if(isgn.eq.2) phcd(nbrn)='s'//phcd1(1:7) print *,'phcd: in ik phcd =',in,ik,' ',phcd(nbrn) if(isw.le.1) go to 17 ik=ik+1 j=max0(j,l) go to 31 17 j=max0(j,l) indx(nseg,2)=nl return end c subroutine pdect(i1,i2,j1,nph,fac) save include 'limits.inc' double precision pb,pu,taup,xp,taul,xl,px,xt,pt,taut,coef,xa double precision h1,h2,hh dimension ib(2,2) common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) common/brnb/px(jbrn,2),xt(jbrn,2),jndx(jbrn,2),mndx(jbrn,2),nbrn common/outc/pt(jout),taut(jout),coef(5,jout),xa(jout),nl c xmn=fac*xmin isg=1 do 1 i=1,2 ib(i,1)=i1 1 ib(i,2)=i2 ii=i1+1 ie=i2-1 xa(i1)=xt(nbrn,1) do 2 i=ii,ie h1=pt(i-1)-pt(i) h2=pt(i+1)-pt(i) hh=h1*h2*(h1-h2) h1=h1*h1 h2=-h2*h2 2 xa(i)=-(h2*taut(i-1)-(h2+h1)*taut(i)+h1*taut(i+1))/hh xa(i2)=xt(nbrn,2) do 3 i=ii,ie if((xa(i+1)-xa(i))*(xa(i)-xa(i-1)).gt.0d0) go to 3 isg=2 ib(1,2)=i-2 ib(2,1)=i+2 3 continue do 4 it=1,isg write(10,*)'collct1',it,isg,ib(it,1),ib(it,2) 4 call collct(ib(it,1),ib(it,2),xa,xmn) k=i1-1 j=j1 do 5 i=i1,i2 if(xa(i).lt.0d0) go to 5 k=k+1 pt(k)=pt(i) taut(k)=taut(i) xa(k)=xa(i) kuse(j,nph)=1 5 j=j+1 if(k.eq.nl) return ii=k+1 do 6 i=ii,nl 6 taut(i)=0d0 nl=k return end c subroutine kseq c c Kseq makes a correspondence between model slownesses in array pb and c the subset of the same slownesses used for sampling tau which are c stored in pu (separate sets for P and S). The net result is to c translate the kndx pointers to critical slowness values (bounding c branches actually implemented) from pointing into pb to pointing c into pu. c c save include 'limits.inc' double precision pb,pu,taup,xp,taul,xl dimension kl(2),kk(jseg,2,2) common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) common/segc/fcs(jseg,3),nafl(jseg,3),indx(jseg,2),kndx(jseg,2), 1 nseg data kl/0,0/ c c Compile a sorted list of unique kndx values in the first column of c kk. c do 1 i=1,nseg nph=iabs(nafl(i,1)) k=kl(nph) do 2 j=1,2 if(k.le.0) go to 3 do 4 m=1,k if(kk(m,nph,1)-kndx(i,j))4,2,5 4 continue 3 k=k+1 kk(k,nph,1)=kndx(i,j) go to 2 5 do 6 l=k,m,-1 6 kk(l+1,nph,1)=kk(l,nph,1) k=k+1 kk(m,nph,1)=kndx(i,j) 2 continue 1 kl(nph)=k c c Make the correspondence between pb and pu for each kndx and save it c in the second column of kk. c do 7 nph=1,2 n1=ku(nph) k=1 ki=kk(k,nph,1) do 8 i=1,n1 if(pu(i,nph)-pb(ki))8,9,10 9 kk(k,nph,2)=i if(k.ge.kl(nph)) go to 7 k=k+1 ki=kk(k,nph,1) 8 continue 10 write(*,100)ki,pb(ki),nph 100 format(1x,'Kseq: pb(',i3,') =',f7.4,' not found in pu(*,',i1, 1 ').') call vexit(1) 7 continue c c Replace each kndx pb pointer with the corresponding pu pointer. c do 11 i=1,nseg nph=iabs(nafl(i,1)) k=kl(nph) do 11 j=1,2 do 12 m=1,k if(kk(m,nph,1)-kndx(i,j))12,11,13 12 continue 13 write(*,101)kndx(i,j) 101 format(1x,'Kseq: kndx value',i4,' not translated.') call vexit(1) 11 kndx(i,j)=kk(m,nph,2) return end c subroutine mseq c partial reordering of tables save include 'limits.inc' double precision pb,pu,taup,xp,taul,xl,px,xt,pux common/intc/pb(nsl1),pu(nsl1,2),taup(nsl1,3,2),xp(nsl1,3,2), 1 taul(nlvz0,2),xl(nlvz0,2),xmin,kb(2),ku(2),lt(2),lvz(nlvz0,2), 2 kuse(nsl1,2) common/segc/fcs(jseg,3),nafl(jseg,3),indx(jseg,2),kndx(jseg,2), 1 nseg common/brnb/px(jbrn,2),xt(jbrn,2),jndx(jbrn,2),mndx(jbrn,2),nbrn common/msqc/pux(jbrn,2),km(2),midx(jbrn,2) data km/0,0/ c is=1 do 1 i=1,nbrn 8 if(jndx(i,2).le.indx(is,2)) go to 7 is=is+1 go to 8 7 nph=iabs(nafl(is,1)) k=km(nph) do 2 j=1,2 if(k.le.0) go to 3 do 4 m=1,k if(midx(m,nph)-mndx(i,j))4,2,5 4 continue 3 k=k+1 midx(k,nph)=mndx(i,j) pux(k,nph)=px(i,j) go to 2 5 do 6 l=k,m,-1 midx(l+1,nph)=midx(l,nph) 6 pux(l+1,nph)=pux(l,nph) k=k+1 midx(m,nph)=mndx(i,j) pux(m,nph)=px(i,j) 2 continue 1 km(nph)=k return end