subroutine input(lu,nlay,ifree) implicit real*8(a-h,o-z) c-LK complex*16 px,py,omega c c........|.........|.........|.........|.........|.........|.........|.. c........|.........|.........|.........|.........|.........|.........|.. c c Rmatrix and associated programs may be used freely by University c researchers and teachers. Commercial parties interested in using c the program should contact C. J. Thomson. c c........|.........|.........|.........|.........|.........|.........|.. c........|.........|.........|.........|.........|.........|.........|.. c c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c-LK: Input logical unit number 10 changed to LU. c-LK: Dummy arguments PX,PY,OMEGA of subroutine INPUT discarded and c-LK: replaced with new dummy argument LU. c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c read(lu,*)nlay,ifree c-LK read(lu,*)omega,px,py do 99 ilay=1,nlay read(lu,*)isoflg(ilay),thickn(ilay),rho(ilay) if(isoflg(ilay).eq.0)then do 100 l=1,3 do 101 k=1,3 read(lu,*)((a(i,j,k,l,ilay),i=1,3),j=1,3) 101 continue 100 continue else read(lu,*)a(3,3,3,3,ilay),a(2,3,2,3,ilay) iflat=0 if(iflat.ne.0)then c-LK print *,' ' c-LK print *,'Earth flattening transformation being called.' c-LK print *,' ' call flatten(ilay,thickn(ilay), . a(3,3,3,3,ilay),a(2,3,2,3,ilay),rho(ilay)) else c-LK print *,' ' c-LK print *,'Earth flattening transformation not being called.' c-LK print *,' ' endif c next lines for conversion from velocity to elasticity a(3,3,3,3,ilay)=(a(3,3,3,3,ilay)**2) a(2,3,2,3,ilay)=(a(2,3,2,3,ilay)**2) endif 99 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine flatten(ilay,thickn,vp,vs,rho) implicit real*8 (a-h,o-z) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. save zlast,rad,rad2,a data a/6371.d3/ c if(ilay.eq.1)then rad=a rad2=a+thickn/2.d0 thickl=thickn zlast=0.d0 endif rad=rad-thickn rad2=rad2-(thickn+thickl)/2.d0 thickl=thickn z=a*dlog(a/rad) thickn=z-zlast zlast=z vp=vp*a/rad2 vs=vs*a/rad2 rho=rho*(rad2/a)**5 c write(11,*)ilay,thickn,vp,vs,rho return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine rmatrix(px,py,omega,nlay,ifree,tus,rus,tds,rds,iter) c c Program to compute r/t matrices for a stack of homogeneous c layers. Uses reflection matrix algorithm of B. L. N. Kennett. c General anisotropy permitted. c c C. J. Thomson, August--November 1993. c c c........|.........|.........|.........|.........|.........|.........|.. c........|.........|.........|.........|.........|.........|.........|.. c c Rmatrix and associated programs may be used freely by University c researchers and teachers. Commercial parties interested in using c the program should contact C. J. Thomson. c c........|.........|.........|.........|.........|.........|.........|.. c........|.........|.........|.........|.........|.........|.........|.. c c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8(a-h,o-z) complex*16 px,py,omega c eigenvalue, eigenvector and inverse eigenvector arrays complex*16 eval1(6), evec1(6,6), eveci1(6,6) complex*16 eval2(6), evec2(6,6), eveci2(6,6) c partitions of wave propagator complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3) c r/t coefficient matrices - single interface and stack of layers complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3) c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) c complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c c start at bottom layer ilay=nlay c only compute e-vecs, etc, if on first pass for this slowness c as they are frequency independent if(iter.eq.1)then call evset(px,py,omega,ilay,eval1,evec1) call estore(eval1,evec1,evals(1,ilay),evecs(1,1,ilay)) endif c invert eigenvalue matrix call xeveci(nlay,ilay,evecs(1,1,ilay),eveci1) c call write2(ilay,evecs(1,1,ilay),eveci1) c set j2 symmetry constant for bottom half space call j2sym(nlay,ilay,evecs(1,1,ilay),j2cstn) c next layer and loop point 999 ilay=ilay-1 c check if on first pass for this slowness if(iter.eq.1)then call evset(px,py,omega,ilay,eval2,evec2) call estore(eval2,evec2,evals(1,ilay),evecs(1,1,ilay)) c wave propagator for interface call qmatrix(eveci1,evec2,q11,q12,q21,q22) call qstore(q11,q12,q21,q22,q11i(1,1,ilay),q12i(1,1,ilay), . q21i(1,1,ilay),q22i(1,1,ilay)) c r/t matrices for interface from wave propagator call retra(q11,q12,q21,q22,tu,ru,td,rd) c call write3(ilay,tu,ru,td,rd,evec1,evec2) call rtstore(tu,ru,td,rd,tui(1,1,ilay),rui(1,1,ilay), . tdi(1,1,ilay),rdi(1,1,ilay)) else call rtget(tu,ru,td,rd,tui(1,1,ilay),rui(1,1,ilay), . tdi(1,1,ilay),rdi(1,1,ilay)) endif c addition rule for previous and present r/t matrices call addit(nlay,ilay,tus,rus,tds,rds,tu,ru,td,rd) call rtstore(tus,rus,tds,rds,tusi(1,1,ilay),rusi(1,1,ilay), . tdsi(1,1,ilay),rdsi(1,1,ilay)) if(ilay.gt.1)then c layer up to next interface c-LK print *,'entering layer, ilay=',ilay call layer(ilay,omega,thickn(ilay),evals(1,ilay), . tus,rus,tds,rds) if(iter.eq.1)then c transfer data call trnsf6(eval1,eval2,evec1,evec2) c invert eigenvalue matrix call xeveci(nlay,ilay,evec1,eveci1) c call write2(ilay,evec1,eveci1) endif c loop go to 999 endif c top layer of stack -- if no free surface above just write out results c-LK call write4(tus,rus,tds,rds) c now test some symmetries of stack results... note requires that the c simple inverse constant of proportionality for the upper half- c space be found first... call xeveci(nlay,ilay,evecs(1,1,ilay),eveci2) c-LK call j1symtst(nlay,ilay,tus,rus,tds,rds) c set j2 symmetry constant for top half space call j2sym(nlay,ilay,evecs(1,1,ilay),j2cst1) c and test j2 symmetry c-LK call j2symtst(nlay,ilay,tus,rus,tds,rds) if(ifree.eq.1)then c free surface at top of first layer -- so include phase shift from c first interface up to free surface c-LK print *,'entering layer, ilay=',ilay call layer(ilay,omega,thickn(ilay),evals(1,ilay), . tus,rus,tds,rds) endif return end c c........|.........|.........|.........|.........|.........|.........|.. c