C *** FOR BUGS OR COMMENTS, PLEASE E-MAIL to "hisada@cc.kogakuin.ac.jp", C OR MAIL TO "Yoshiaki Hisada, Dept. of Architecture, C Kogakuin Univeristy, Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan" c c *** This code compute phase velocities with small amplitude of the secular c functions, which are used for computing Green's functions. c c *** Note that some of them are not the phase velocities of the surface waves, c but the body waves. c c *** This code is open only for the academic use. c c (originally made in 1992, and modified in Feb, 1996 from "normal.f") C PROGRAM phs3sQ C IMPLICIT REAL*8(A-H,O-Z) PARAMETER(MOM=700) COMMON /OUTPUT/ISEC,IVCT C WRITE(*,1000) 1000 FORMAT(//' "phs3.f" computes the normal mode solutions'/ * ' (the phase and group velocities, and eigen vectors of'/ * ' the Love and Rayleigh waves).') WRITE(*,*)' For trouble, please check "fort.10".' WRITE(*,*)' This code is open for the academic use.' WRITE(*,*) *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".' WRITE(*,*) *' (Originally made in 1992, and modified in 1995)' C CALL INDAT(MOM) C CALL PHASE(MOM) CALL ENGUU(MOM) C WRITE(*,*) WRITE(*,*)'Data "mdl_dat, mdr_dat, phsl_dat, and phsr_dat"' WRITE(*,*)'are used for "grflt4.f" as data of poles' WRITE(*,*) WRITE(*,*)'You can plot the dispersion curves using "phsl.m"' WRITE(*,*)'and "phsr.m" for "M-file" of "Matlab"' WRITE(*,*) write(*,*)'The followings are output data ofExcel format:' WRITE(*,*)' Ldisper.csv and Rdisper.csv: the dipersion curves' c if(ISEC.eq.1)then WRITE(*,*)' sec-Lov.csv: Secular Functions of Love Wave' WRITE(*,*)' sec-Ray.csv: Secular Functions of Rayleigh Wave' WRITE(*,*)' Ldis.csv: Eigen Displacement Vector of Love Wave' WRITE(*,*)' Lstr.csv: Eigen Stress Vector of Love Wave' WRITE(*,*)' Rdis.csv: Eigen Displacement Vector of Rayleigh Wave' WRITE(*,*)' Rstr.csv: Eigen Stress Vector of Rayleigh Wave' write(*,*)' Lmedres.csv: Medium Response of Love Wave' write(*,*)' Rmedres.csv: Medium Response of Rayleigh Wave' end if C STOP END C C C *** DATA INPUT *** C SUBROUTINE INDAT(MOM) C IMPLICIT REAL*8(A-H,O-Z) character infile*20 COMMON /FREQ/DOM,IOM,NOM * /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL * /MINMAX/VSMIN,VSMAX,CCMIN,CCMAX,CMIN,CMAX,TOL,NITER,NCC * /OUTPUT/ISEC,IVCT /NINTG/ZINT,NZINT C C *** PARAMETERS IN COMMOM STATEMENT *** C NML =20 C C *** OMEGA DATA *** C write(*,*) write(*,*)'Enter Input File Name for "grflt12sx-in.csv"' read(*,'(a20)')infile c infile='grflt04bx-in.csv' OPEN (UNIT=7,FILE=infile,STATUS='OLD') C read(7,*) read(7,*) read(7,*)dt,nt write(*,*) write(*,*)'Delta Time (s) =',dt,', Number of Time =',nt a=nt 1 a=a/2 if(a.gt.1.d0)then goto 1 else if(a.ne.1.d0)then write(*,*) write(*,*)'NT is not a power of 2 !!' write(10,*)'NT is not a power of 2 !!' stop end if end if read(7,*) read(7,*)period c write(*,*)'Minimum Period (s) =',period duration=dt*nt pi2=datan(1.d0)*8.d0 dom=pi2/duration iom=0 nom=int(duration/period) c write(*,*) write(*,*)'Maximum frequency (Hz) =',1.d0/period write(*,*)'Check Secular Functions & Eigen Vectors? (N=0 or Y=1)' read(*,*)ichk if(ichk.eq.1)then nom=1 iom=1 write(*,*)'Enter the frequency for checking (Hz)' read(*,*)fchk dom=fchk*pi2 end if c WRITE(10,*)'FROM indat:' WRITE(10,*) WRITE(10,*)'dt,nt,period',dt,nt,period WRITE(10,*)'IOM,NOM,DOM=',IOM,NOM,DOM C IF(NOM.GT.MOM)THEN WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!' STOP END IF C C *** MEDIUM DATA *** C read(7,*) read(7,*) READ(7,*)NL read(7,*) write(*,*)'Number of Layers=',NL write(*,*)'Layer, Density, Vp, Qp, Fp, Vs, Qs, Fs, Thickness' DO 100 IL=1,NL READ(7,*)i,DNS(IL),VP(IL),QP,FP,VS(IL),QS,FS,THK(IL) write(*,*)i,DNS(IL),VP(IL),QP,FP,VS(IL),QS,FS,THK(IL) 100 CONTINUE C WRITE(10,*)'NL=',NL DO 150 IL=1,NL WRITE(10,*)IL,DNS(IL),VP(IL),VS(IL),THK(IL) 150 CONTINUE C IF(NL.GT.NML)THEN WRITE(10,*)' NL>NML !!!!!!!!!!!!!' STOP END IF C CLOSE(7) C C *** MIN & MAX VS AND CC *** C VSMIN=VS(NL) VSMAX=0.0D0 VPMAX=0.0D0 MINIL=1 MAXIL=1 DO 200 IL=1,NL IF(VS(IL).LE.VSMIN)THEN VSMIN=VS(IL) MINIL=IL END IF IF(VS(IL).GE.VSMAX)THEN VSMAX=VS(IL) MAXIL=IL END IF VPMAX=MAX(VPMAX,VP(IL)) 200 CONTINUE C IF(VSMAX.NE.VS(NL))THEN WRITE(*,*)' FROM SUBROUTINE INDAT :' WRITE(*,*)' VS(NL).NE.VSMAX, THUS STOP !!!' STOP END IF IF(VPMAX.NE.VP(NL))THEN WRITE(*,*)' FROM SUBROUTINE INDAT :' WRITE(*,*)' VP(NL).NE.VPMAX, THUS STOP !!!' STOP END IF C WRITE(*,1602) READ(*,*)TOL c 1602 FORMAT(/' Enter Tolerance for Secular Function (ex., 0.01):') C WRITE(*,1603) READ(*,*)NITER 1603 FORMAT(/' Enter Max Number of Iterations (ex., 200):') C CCMIN=FRHF(VP(MINIL),VS(MINIL)) CCMAX=FRHF(VP(MAXIL),VS(MAXIL)) C WRITE(*,1600)VSMAX,VPMAX WRITE(*,1610)CCMIN,CCMAX 1600 FORMAT(/' VSMAX, VPMAX =',2F9.3) 1610 FORMAT( ' CCMIN, CCMAX =',2F9.3) c NCC=1 IF(NL.EQ.1)GOTO 25 c WRITE(*,1620) c READ(*,*)IVC c IF(IVC.EQ.1)THEN c WRITE(*,1630) c READ(*,*)CMIN,CMAX c ELSE CMIN=CCMIN-CCMIN/50.D0 CMAX=VSMAX+VSMAX/100.D0 c END IF c 1620 FORMAT( ' Change CMIN and CMAX? (yes=1):') c 1630 FORMAT( ' Enter CMIN and CMAX:') c WRITE(*,1640)CMIN,CMAX 1640 FORMAT(/' Min and Max Velocities (Cmin and Cmax: m/s) =',2F9.3) c WRITE(*,1650) 1650 FORMAT *( ' Enter Number of Partitions from Cmin to Cmax (ex.,1000):') READ(*,*)NCC DCC=(CMAX-CMIN)/NCC WRITE(*,1660)DCC 1660 FORMAT( ' Increment of Phase Velocity (DCC: m/s) =',F10.5) c ISEC=0 NZINT=1 write(*,*)NOM,IOM IF((NOM-IOM).ne.0)return IF(IOM.EQ.0)RETURN c c WRITE(*,1680) c 1680 FORMAT( ' Output Secular Functions? (yes=1) :') c READ(*,*)ISEC ISEC=1 IF(ISEC.EQ.1)THEN WRITE(*,*) WRITE(*,*)' You can plot the secular functions using' WRITE(*,*)' "secl.m and secr.m" for M-files of "Matlab"' END IF C 25 NZINT=1 IF((NOM-IOM).ne.0)return c WRITE(*,1800) c READ(*,*)IVCT IVCT=1 IF(IVCT.NE.1)RETURN WRITE(*,1810) READ(*,*)ZINT IF(ZINT.LT.THK(1))THEN WRITE(*,*) WRITE(*,*)' Depth is smaller than the thichness of 1st layer !' END IF WRITE(*,1820) READ(*,*)NZINT WRITE(*,*) IF(IVCT.EQ.1)THEN WRITE(*,*) WRITE(*,*)' You can plot the eigen vectors using "vld.m, vrd.m' WRITE(*,*)' vls.m, and vrs.m" for "M-file" of "Matlab"' END IF c 1800 FORMAT(/' Output Eigen Vectors? (yes=1) :') 1810 FORMAT(/' Enter Depth to Output the Eigen Vectors (m):') 1820 FORMAT( ' Enter Partition Number along the Depth (ex., 200):') C RETURN END C C C *** PHASE VELOCITY OF THE RAYLEIGH WAVE IN THE HOMOGENEOUS HALF-SPACE *** C FUNCTION FRHF(VP,VS) IMPLICIT REAL*8(A-H,O-Z) C DZ=0.0001d0 C SP=VS/VP IF(SP.GE.0.0D0) Z= 0.91263D0 IF(SP.GE.0.2D0) Z= 0.90767D0 IF(SP.GE.0.4D0) Z= 0.88899D0 IF(SP.GE.0.6D0) Z= 0.83575D0 IF(SP.GE.0.7D0) Z= 0.77062D0 IF(SP.GE.0.8D0) Z= 0.63338D0 C F0=FRHF1(SP,Z) 100 Z=Z-DZ F1=FRHF1(SP,Z) AF=F0/F1 F00=F0 F0 =F1 IF(AF.GE.0.D0)GOTO 100 FZ=Z+DZ*F0/(F0-F00) C FRHF=VS*DSQRT(FZ) C RETURN END C C FUNCTION FRHF1(SP,Z) IMPLICIT REAL*8(A-H,O-Z) C Z2=Z*Z Z3=Z2*Z C F1=16.D0*(Z-1.D0)*SP*SP FRHF1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0 C RETURN END C c C *** COMPUTATION OF PHASE VELOCITY BY THE BI-SECTION METHOD *** C SUBROUTINE PHASE(MOM) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 GRU0L,GTDL,GRDL,GTDR,GRDR,SCL,SCR,AA(2,2) * ,SC1,SC2,SC3,CMS,CMP,ER(4,4),DDS,DDP COMMON /FREQ/DOM,IOM,NOM * /MINMAX/VSMIN,VSMAX,CCMIN,CCMAX,CMIN,CMAX,TOL,NITER,NCC * /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL * /GRTLV/GRU0L,GTDL(20),GRDL(20) * /GRTRY/GTDR(2,2,20),GRDR(2,2,20) * /CCMD/CCL(0:701,0:230),CCR(0:701,0:230),MDL(0:701),MDR(0:701) * /OUTPUT/ISEC,IVCT C C *** CHECK OF DIMENSION *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM PHASE: NL.GT.20 !!!!!!!!!!!!!!!' WRITE(10,*)' FROM PHASE: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C NC=701 IF(MOM.GT.NC)THEN WRITE(*,*)' FROM PHASE: MOM.GT.701 !!!!!!!!!!!!!!!' WRITE(10,*)' FROM PHASE: MOM.GT.701 !!!!!!!!!!!!!!!' STOP END IF C C *** PHASE VELOCITY *** C OPEN(20,FILE='om_dat',STATUS='unknown') C DCC=(CMAX-CMIN)/NCC write(*,*) write(10,*) DO 100 IO=IOM,NOM C OM=DOM*DFLOAT(IO) if(io.eq.0)om=0.001d0 write(20,1100)om write(*,*) write(*,1110)om,om/6.283185307d0 write(10,1110)om,om/6.283185307d0 1100 FORMAT(e15.7) 1110 FORMAT(' omega and freq.=',f8.3,f8.4) c c IF(IO.EQ.0)THEN c MDL(IO)=0 c CCL(IO,0)=VSMAX c GOTO 20 c END IF C IF(ISEC.EQ.1)THEN OPEN(11,FILE='sclc_dat',STATUS='unknown') OPEN(12,FILE='sclr_dat',STATUS='unknown') OPEN(13,FILE='scli_dat',STATUS='unknown') OPEN(14,FILE='scrc_dat',STATUS='unknown') OPEN(15,FILE='scrr_dat',STATUS='unknown') OPEN(16,FILE='scri_dat',STATUS='unknown') open(17,file='sec-Lov.csv',STATUS='unknown') open(18,file='sec-Ray.csv',STATUS='unknown') write(17,*)'Phase Velocity (m/s), Real Part, Imag Part *,Absolute Value' write(18,*)'Phase Velocity (m/s), Real Part, Imag Part *,Absolute Value' END IF C IF(NL.EQ.1)GOTO 20 C C ** SEARCH OF LOVE POLES BY THE BI-SECTION METHOD ** C MD=0 C C1=CMIN WN=OM/C1 CALL GNRTLV(OM,WN) SC1=1.D0-GRU0L*GRDL(1) AS1=CDABS(SC1) C IF(ISEC.EQ.1)then WRITE(11,1100)C1 WRITE(12,1100)dreal(sc1) WRITE(13,1100)dimag(sc1) write(17,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1) end if 1101 FORMAT(e15.7,3(',',e15.7)) C C2=C1+DCC WN=OM/C2 CALL GNRTLV(OM,WN) SC2=1.D0-GRU0L*GRDL(1) AS2=CDABS(SC2) C 10 C3=C2+DCC if(C3.gt.(CMAX+DCC))goto 15 WN=OM/C3 CALL GNRTLV(OM,WN) SC3=1.D0-GRU0L*GRDL(1) AS3=CDABS(SC3) C C1N=C1 C2N=C2 C3N=C3 AS2N=AS2 c c write(2,9999)c1n,as1,c2n,as2,c3n,as3 c 9999 format(3(f10.4,f10.6)) IF((AS2.LT.AS1).AND.(AS2.LT.AS3))THEN C c write(2,*)0,c1n,as1 c write(2,*)0,c2n,as2 c write(2,*)0,c3n,as3 c do 110 iter=1,NITER c C12=(C1N+C2N)/2.D0 WN=OM/C12 CALL GNRTLV(OM,WN) SCL=1.D0-GRU0L*GRDL(1) AS12=CDABS(SCL) c write(2,*)iter,c12,as12 c if(as12.eq.as2n)goto 12 c C23=(C2N+C3N)/2.D0 WN=OM/C23 CALL GNRTLV(OM,WN) SCL=1.D0-GRU0L*GRDL(1) AS23=CDABS(SCL) c write(2,*)iter,c23,as23 c write(2,*)100,c12,as12 c write(2,*)200,c2n,as2n c write(2,*)300,c23,as23 c IF(AS12.LT.AS2N)THEN C3N =C2N C2N =C12 AS2N=AS12 ELSE IF(AS23.LT.AS2N)THEN C1N =C2N C2N =C23 AS2N=AS23 ELSE C1N =C12 C3N =C23 END IF END IF c IF(AS2N.LT.TOL)THEN CCL(IO,MD)=C2N MD=MD+1 goto 14 END IF C 110 CONTINUE C 12 write(*,*)'Love Wave: Did not converge!' write(10,*)'Love Wave: Did not converge!' write(*,1200)MD,iter,C2N,AS2N write(10,1200)MD,iter,C2N,AS2N 1200 format(' Mode,Iter,Phase Vel,Amp=',i3,i5,f9.2,f9.3) c 14 C1=C2N WN=OM/C1 CALL GNRTLV(OM,WN) SC1=1.D0-GRU0L*GRDL(1) AS1=CDABS(SC1) c IF(ISEC.EQ.1)THEN WRITE(11,1100)C1 WRITE(12,1100)dreal(sc1) WRITE(13,1100)dimag(sc1) write(17,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1) END IF C C2=C1+DCC WN=OM/C2 CALL GNRTLV(OM,WN) SC2=1.D0-GRU0L*GRDL(1) AS2=CDABS(SC2) C ELSE c IF(ISEC.EQ.1)THEN WRITE(11,1100)C2 WRITE(12,1100)dreal(sc2) WRITE(13,1100)dimag(sc2) write(17,1101)C2,dreal(sc2),dimag(sc2),cdabs(sc2) END IF C C1=C2 SC1=SC2 AS1=AS2 C C2=C3 SC2=SC3 AS2=AS3 C END IF goto 10 C 15 MDL(IO)=MD-1 NC=230 IF(MDL(IO).GT.NC)THEN WRITE(*,*)' FROM PHASE: LOVE MODE.GT.230 !!!!!!!!!!!!!!' STOP END IF C C *** SEARCH OF RAYLEIGH POLES BY THE BI-SECTION METHOD *** C 20 MD=0 C write(*,*)'1;R start,om',om IF(NL.EQ.1)THEN CCR(IO,MD)=CCMIN GO TO 100 END IF c IF(IO.EQ.0)THEN MDR(IO)=0 CCR(IO,0)=CCMAX GOTO 100 END IF C SM2=DNS(1)*VS(1)*VS(1)*DNS(1)*VS(1)*VS(1) C C1=CMIN IF(dabs(C1-VS(1)).lt.1.0d-7)C1=C1+C1/1000.D0 IF(dabs(C1-VP(1)).lt.1.0d-7)C1=C1+C1/1000.D0 WN=OM/C1 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1)+ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1)+ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1)+ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1)+ER(4,4)*DDS*GRDR(2,2,1) SC1=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS1=CDABS(SC1) AMAX=AS1 C IF(ISEC.EQ.1)then WRITE(14,1100)C1 WRITE(15,1100)dreal(sc1) WRITE(16,1100)dimag(sc1) write(18,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1) end if C C2=C1+DCC IF(dabs(C2-VS(1)).lt.1.0d-7)C2=C2+C2/1000.D0 IF(dabs(C2-VP(1)).lt.1.0d-7)C2=C2+C2/1000.D0 WN=OM/C2 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1)+ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1)+ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1)+ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1)+ER(4,4)*DDS*GRDR(2,2,1) SC2=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS2=CDABS(SC2) IF(AS2.GT.AMAX)AMAX=AS2 C 30 C3=C2+DCC if(C3.gt.(CMAX+DCC))goto 35 IF(dabs(C3-VS(1)).lt.1.0d-7)C3=C3+C3/1000.D0 IF(dabs(C3-VP(1)).lt.1.0d-7)C3=C3+C3/1000.D0 WN=OM/C3 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) * +ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) * +ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) * +ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) * +ER(4,4)*DDS*GRDR(2,2,1) SC3=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS3=CDABS(SC3) IF(AS3.GT.AMAX)AMAX=AS3 C C1N=C1 C2N=C2 C3N=C3 AS2N=AS2 c c write(2,9999)c1n,as1,c2n,as2,c3n,as3 c 9999 format(3(f10.4,f10.6)) C IF((AS2.LT.AS1).AND.(AS2.LT.AS3))THEN c c write(2,*)0,c1n,as1 c write(2,*)0,c2n,as2 c write(2,*)0,c3n,as3 c do 310 iter=1,NITER C C12=(C1N+C2N)/2.D0 IF(dabs(C12-VS(1)).lt.1.0d-7)C12=C12+C12/1000.D0 IF(dabs(C12-VP(1)).lt.1.0d-7)C12=C12+C12/1000.D0 WN=OM/C12 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) * +ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) * +ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) * +ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) * +ER(4,4)*DDS*GRDR(2,2,1) SCR=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS12=CDABS(SCR) c write(2,*)iter,c12,as12 c if(as12.eq.as2n)goto 32 C C23=(C2N+C3N)/2.D0 IF(dabs(C23-VS(1)).lt.1.0d-7)C23=C23+C23/1000.D0 IF(dabs(C23-VP(1)).lt.1.0d-7)C23=C23+C23/1000.D0 WN=OM/C23 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) * +ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) * +ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) * +ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) * +ER(4,4)*DDS*GRDR(2,2,1) SCR=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS23=CDABS(SCR) c c write(2,*)iter,c23,as23 c write(2,*)100,c12,as12 c write(2,*)200,c2n,as2n c write(2,*)300,c23,as23 C IF(AS12.LT.AS2N)THEN C3N =C2N C2N =C12 AS2N=AS12 ELSE IF(AS23.LT.AS2N)THEN C1N =C2N C2N =C23 AS2N=AS23 ELSE C1N =C12 C3N =C23 END IF END IF C IF(AS2N.LT.TOL)THEN CCR(IO,MD)=C2N MD=MD+1 GOTO 34 END IF C 310 continue C 32 write(*,*)'Rayleigh Wave: Did not converge!' write(10,*)'Rayleigh Wave: Did not converge!' write(*,3200)MD,iter,C2N,AS2N write(10,3200)MD,iter,C2N,AS2N 3200 format(' Mode,Iter,Phase Vel,Amp=',i3,i5,f9.2,f9.3) C 34 C1=C2N WN=OM/C1 IF(C1.EQ.VS(1))C1=C1+C1/1000.D0 IF(C1.EQ.VP(1))C1=C1+C1/1000.D0 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) * +ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) * +ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) * +ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) * +ER(4,4)*DDS*GRDR(2,2,1) SC1=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS1=CDABS(SC1) IF(ISEC.EQ.1)THEN WRITE(14,1100)C1 WRITE(15,1100)dreal(sc1) WRITE(16,1100)dimag(sc1) write(18,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1) END IF C C2=C1+DCC WN=OM/C2 IF(C2.EQ.VS(1))C2=C2+C2/1000.D0 IF(C2.EQ.VP(1))C2=C2+C2/1000.D0 CALL GNRTRY(OM,WN) CALL MTXRY(OM,WN,1,CMS,CMP,ER) DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) * +ER(3,4)*DDS*GRDR(2,1,1) AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) * +ER(3,4)*DDS*GRDR(2,2,1) AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) * +ER(4,4)*DDS*GRDR(2,1,1) AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) * +ER(4,4)*DDS*GRDR(2,2,1) SC2=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) AS2=CDABS(SC2) c ELSE C IF(ISEC.EQ.1)then WRITE(14,1100)C2 WRITE(15,1100)dreal(sc2) WRITE(16,1100)dimag(sc2) write(18,1101)C2,dreal(sc2),dimag(sc2),cdabs(sc2) end if C C1=C2 SC1=SC2 AS1=AS2 C C2=C3 SC2=SC3 AS2=AS3 C END IF C GOTO 30 C 35 MDR(IO)=MD-1 NC=230 IF(MDR(IO).GT.NC)THEN WRITE(*,*)' FROM PHASE: RAYLEIGH MODE.GT.230 !!!!!!!!!!!!!!!' STOP END IF C write(*,*) * ' The highest modes for L- and R-waves',mdl(io),mdr(io) C 100 CONTINUE C CLOSE(11) CLOSE(12) CLOSE(13) CLOSE(14) CLOSE(15) CLOSE(16) CLOSE(17) CLOSE(18) CLOSE(20) c OPEN(22,FILE='mdl_dat',STATUS='unknown') OPEN(23,FILE='mdr_dat',STATUS='unknown') DO 300 IO=IOM,NOM WRITE(22,1000)MDL(IO) WRITE(23,1000)MDR(IO) 300 CONTINUE 1000 FORMAT(I5) CLOSE(22) CLOSE(23) c OPEN(25,FILE='phsl1_dat',STATUS='unknown') OPEN(26,FILE='phsr1_dat',STATUS='unknown') DO 210 IO=IOM,NOM DO 200 MD=0,MDL(IO) WRITE(25,1100)CCL(IO,MD) 200 CONTINUE DO 205 MD=0,MDR(IO) WRITE(26,1100)CCR(IO,MD) 205 CONTINUE 210 CONTINUE CLOSE(25) CLOSE(26) c RETURN END C C C *** COMPUTATION OF EIGEN-VECTORS, ENERGY INTEGRALS AND GROUP VELOCITIES *** C SUBROUTINE ENGUU(MOM) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 GRU0L,GTDL,GRDL,GTDR,GRDR * ,EL(2,2),CDS,CUS,DDS,DUS,CMS,CDPU,VR20 * ,ER(4,4),CDP,CUP,DDP,DUP,CMP,DR1,DR2,DR14,DR23 * ,A1,A2,A3,A4,A5,A6,B1,B2,B3,B4,B5,B6 dimension res(0:230) COMMON /FREQ/DOM,IOM,NOM * /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL * /OUTPUT/ISEC,IVCT /NINTG/ZINT,NZINT * /GRTLV/GRU0L,GTDL(20),GRDL(20) * /GRTRY/GTDR(2,2,20),GRDR(2,2,20) * /CCMD/CCL(0:701,0:230),CCR(0:701,0:230),MDL(0:701),MDR(0:701) * /VCTL/VL1(702,0:230),VL2(702,0:230) * /VCLR/VR1(702,0:230),VR2(702,0:230),VR3(702,0:230),VR4(702,0:230) * /UEGL/UUL(0:701,0:230),E1L(0:701,0:230) * /UEGR/UUR(0:701,0:230),E1R(0:701,0:230),HoV(0:701,0:230) * /MINMAX/VSMIN,VSMAX,CCMIN,CCMAX,CMIN,CMAX,TOL,NITER,NCC C C *** CHECK OF DIMENSION *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM ENGUU: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C NC=702 IF((NZINT+2).GT.NC)THEN WRITE(*,*)' FROM ENGUU: NZINT+2.GT.702 !!!!!!!!!!!!!!!' STOP END IF C NC=701 IF(MOM.GT.NC)THEN WRITE(*,*)' FROM ENGUU: MOM.GT.701 !!!!!!!!!!!!!!!' STOP END IF C C *** EIGEN VECTORS FOR LOVE WAVE *** C IF(IVCT.EQ.1)THEN OPEN(9,FILE='vzn_dat',STATUS='unknown') write(9,1000)zint write(9,1100)nzint write(9,1100)nl do 90 l=1,nl write(9,1000)thk(l) 90 continue close(9) OPEN(10,FILE='vl1_dat',STATUS='unknown') OPEN(11,FILE='vl2_dat',STATUS='unknown') OPEN(20,FILE='vr1_dat',STATUS='unknown') OPEN(21,FILE='vr2_dat',STATUS='unknown') OPEN(22,FILE='vr3_dat',STATUS='unknown') OPEN(23,FILE='vr4_dat',STATUS='unknown') OPEN(24,FILE='Ldis.csv',STATUS='unknown') OPEN(25,FILE='Lstr.csv',STATUS='unknown') OPEN(26,FILE='Rdis.csv',STATUS='unknown') OPEN(27,FILE='Rstr.csv',STATUS='unknown') WRITE(24,2024)(MD,MD=MDL(NOM),0,-1) WRITE(25,2024)(MD,MD=MDL(NOM),0,-1) WRITE(26,2025)(MD,MD,MD=MDR(NOM),0,-1) WRITE(27,2025)(MD,MD,MD=MDR(NOM),0,-1) END IF 2024 FORMAT('depth (km)',400(',',I3,' mode')) 2025 FORMAT('depth (km)',400(',',I3,' mode (h),',I3,' mode (v)')) 1000 FORMAT(e15.7) 1100 format(i5) C DZ=ZINT/NZINT DO 100 IO=IOM,NOM C OM=DOM*DFLOAT(IO) IF(IO.EQ.0)OM=DOM/1000.d0 OM2=OM*OM C IF(NL.EQ.1)GOTO 20 c IF(IO.EQ.0)THEN UUL(IO,0)=VSMAX c GOTO 20 END IF C DO 110 MD=0,MDL(IO) C WN=OM/CCL(IO,MD) WN2=WN*WN CALL GNRTLV(OM,WN) C IF(IVCT.NE.1)GOTO 10 C LAY=1 UDEP=0.D0 BDEP=THK(1) C VS2=VS(1)*VS(1) WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS SM=DNS(1)*VS2 EL(2,1)=-SM*CMS EL(2,2)= SM*CMS C CDS=(0.5D0,0.D0) CUS=GRDL(1)*0.5D0 C DO 120 IZ=1,NZINT+1 C Z=DZ*(IZ-1) IF(Z.GT.BDEP)THEN C LAY=LAY+1 UDEP=BDEP IF(LAY.EQ.NL)THEN BDEP=BDEP+ZINT ELSE BDEP=BDEP+THK(LAY) END IF VS2=VS(LAY)*VS(LAY) WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS SM=DNS(LAY)*VS2 EL(2,1)=-SM*CMS EL(2,2)= SM*CMS C CDS=GTDL(LAY-1)*CDS CUS=GRDL(LAY)*CDS C END IF C DDS=CDEXP(-CMS*(Z-UDEP)) IF(LAY.EQ.NL)THEN VL1(IZ,MD)=dreal(DDS*CDS) VL2(IZ,MD)=dreal(EL(2,1)*DDS*CDS) ELSE DUS=CDEXP(-CMS*(BDEP-Z)) VL1(IZ,MD)=dreal(DDS*CDS+DUS*CUS) VL2(IZ,MD)=dreal(EL(2,1)*DDS*CDS+EL(2,2)*DUS*CUS) END IF C 120 CONTINUE C DO 140 I=1,NZINT+1 WRITE(10,1000)VL1(I,MD) WRITE(11,1000)VL2(I,MD) 140 CONTINUE C C *** ENERGY INTEGRALS AND GROUP VELOCITIES FOR LOVE WAVE *** C 10 E1L(IO,MD)=0.D0 E2L =0.D0 C DO 160 LAY=1,NL C IF(LAY.EQ.1)THEN CDS=(0.5D0,0.D0) CUS=GRDL(1)*0.5D0 ELSE CDS=GTDL(LAY-1)*CDS CUS=GRDL(LAY)*CDS END IF VS2=VS(LAY)*VS(LAY) WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS C IF(LAY.EQ.NL)THEN DR1=CDS*CDS/(2.D0*CMS) ELSE DDS=CDEXP(-2.D0*CMS*THK(LAY)) DR1=(CDS*CDS+CUS*CUS)*(1.D0-DDS)/(2.D0*CMS) DDS=CDEXP(-CMS*THK(LAY)) DR1=DR1+2.D0*CDS*CUS*THK(LAY)*DDS END IF C SM=DNS(LAY)*VS2 E1L(IO,MD)=dreal(E1L(IO,MD)+DNS(LAY)*DR1) E2L =dreal(E2L +SM*DR1) C 160 CONTINUE C E1L(IO,MD)=E1L(IO,MD)/2.D0 E2L =E2L /2.D0 UUL(IO,MD)=E2L/(E1L(IO,MD)*CCL(IO,MD)) C 110 CONTINUE c do MD=0,MDL(IO) VDIS=0.001d0 VSTR=0.001d0 do IZ=1,NZINT+1 if(dabs(VL1(IZ,MD)).gt.VDIS)VDIS=dabs(VL1(IZ,MD)) if(dabs(VL2(IZ,MD)).gt.VSTR)VSTR=dabs(VL2(IZ,MD)) end do do IZ=1,NZINT+1 VL1(IZ,MD)=VL1(IZ,MD)/VDIS VL2(IZ,MD)=VL2(IZ,MD)/VSTR end do end do c IF(IVCT.EQ.1)THEN do IZ=1,NZINT+1 Z=DZ*(IZ-1) write(24,2500)z/1000.d0,(VL1(IZ,MD)/2.d0+MD,MD=MDL(IO),0,-1) write(25,2500)z/1000.d0,(VL2(IZ,MD)/2.d0+MD,MD=MDL(IO),0,-1) end do end if C C *** EIGEN VECTORS FOR RAYLEIGH WAVE *** C 20 IF(IO.EQ.0)THEN UUR(IO,0)=CCMAX c goto 100 END IF c DO 210 MD=0,MDR(IO) C WN=OM/CCR(IO,MD) CALL GNRTRY(OM,WN) C LAY=1 UDEP=0.D0 CALL MTXRY(OM,WN,1,CMS,CMP,ER) C C * HOMOGENEOUS HALF-SPACE SOLUTION * C IF(NL.EQ.1)THEN C B1=-ER(3,1)/ER(3,2) B2= ER(2,1)+ER(2,2)*B1 CDP= 1.D0/B2 CDS= B1*CDP C C NNZ=NZINT+1 IF(IVCT.NE.1)NNZ=1 DO 215 IZ=1,NNZ C Z=DZ*(IZ-1) DDS=CDEXP(-CMS*Z) DDP=CDEXP(-CMP*Z) A1=DDP*CDP A2=DDS*CDS VR1(IZ,MD)=dreal(ER(1,1)*A1+ER(1,2)*A2) VR2(IZ,MD)=dreal(ER(2,1)*A1+ER(2,2)*A2) VR3(IZ,MD)=dreal(ER(3,1)*A1+ER(3,2)*A2) VR4(IZ,MD)=dreal(ER(4,1)*A1+ER(4,2)*A2) C 215 CONTINUE C C * LAYERED HALF-SPACE SOLUTION * C ELSE C DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) A1=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1)+ER(3,4)*DDS*GRDR(2,1,1) A2=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1)+ER(3,4)*DDS*GRDR(2,2,1) B1=-A1/A2 C CUP=GRDR(1,1,1)+GRDR(1,2,1)*B1 CUS=GRDR(2,1,1)+GRDR(2,2,1)*B1 A3=DDP*CUP A4=DDS*CUS VR20=ER(2,1)+ER(2,2)*B1+ER(2,3)*A3+ER(2,4)*A4 C CDP=1.D0/VR20 CDS=B1*CDP CUP=GRDR(1,1,1)*CDP+GRDR(1,2,1)*CDS CUS=GRDR(2,1,1)*CDP+GRDR(2,2,1)*CDS C BDEP=THK(1) NNZ=NZINT+1 IF(IVCT.NE.1)NNZ=1 DO 220 IZ=1,NNZ C Z=DZ*(IZ-1) IF(Z.GT.BDEP)THEN C LAY=LAY+1 UDEP=BDEP IF(LAY.EQ.NL)THEN BDEP=BDEP+ZINT ELSE BDEP=BDEP+THK(LAY) END IF CALL MTXRY(OM,WN,LAY,CMS,CMP,ER) C CDPU=CDP CDP=GTDR(1,1,LAY-1)*CDPU+GTDR(1,2,LAY-1)*CDS CDS=GTDR(2,1,LAY-1)*CDPU+GTDR(2,2,LAY-1)*CDS CUP=GRDR(1,1,LAY)*CDP+GRDR(1,2,LAY)*CDS CUS=GRDR(2,1,LAY)*CDP+GRDR(2,2,LAY)*CDS C END IF C DDS=CDEXP(-CMS*(Z-UDEP)) DDP=CDEXP(-CMP*(Z-UDEP)) A1=DDP*CDP A2=DDS*CDS IF(LAY.EQ.NL)THEN VR1(IZ,MD)=dreal(ER(1,1)*A1+ER(1,2)*A2) VR2(IZ,MD)=dreal(ER(2,1)*A1+ER(2,2)*A2) VR3(IZ,MD)=dreal(ER(3,1)*A1+ER(3,2)*A2) VR4(IZ,MD)=dreal(ER(4,1)*A1+ER(4,2)*A2) ELSE DUS=CDEXP(-CMS*(BDEP-Z)) DUP=CDEXP(-CMP*(BDEP-Z)) A3=DUP*CUP A4=DUS*CUS VR1(IZ,MD)=dreal(ER(1,1)*A1+ER(1,2)*A2+ER(1,3)*A3+ER(1,4)*A4) VR2(IZ,MD)=dreal(ER(2,1)*A1+ER(2,2)*A2+ER(2,3)*A3+ER(2,4)*A4) VR3(IZ,MD)=dreal(ER(3,1)*A1+ER(3,2)*A2+ER(3,3)*A3+ER(3,4)*A4) VR4(IZ,MD)=dreal(ER(4,1)*A1+ER(4,2)*A2+ER(4,3)*A3+ER(4,4)*A4) END IF C 220 CONTINUE C END IF C IF(IVCT.NE.1)goto 30 DO 240 I=1,NZINT+1 WRITE(20,1000)VR1(I,MD) WRITE(21,1000)VR2(I,MD) WRITE(22,1000)VR3(I,MD) WRITE(23,1000)VR4(I,MD) 240 CONTINUE C C *** ENERGY INTEGRAL, GROUP VELOCITY and H over V FOR RAYLEIGH WAVE *** C 30 E1R(IO,MD)=0.D0 E2R =0.D0 E3R =0.D0 C DO 260 LAY=1,NL C CALL MTXRY(OM,WN,LAY,CMS,CMP,ER) C IF(NL.EQ.1)THEN B1=-ER(3,1)/ER(3,2) B2= ER(2,1)+ER(2,2)*B1 CDP= 1.D0/B2 CDS= B1*CDP ELSE IF(LAY.EQ.1)THEN DDS=CDEXP(-CMS*THK(1)) DDP=CDEXP(-CMP*THK(1)) A1=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) * +ER(3,4)*DDS*GRDR(2,1,1) A2=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) * +ER(3,4)*DDS*GRDR(2,2,1) B1=-A1/A2 CUP=GRDR(1,1,1)+GRDR(1,2,1)*B1 CUS=GRDR(2,1,1)+GRDR(2,2,1)*B1 A3=DDP*CUP A4=DDS*CUS VR20=ER(2,1)+ER(2,2)*B1+ER(2,3)*A3+ER(2,4)*A4 CDP=1.D0/VR20 CDS=B1*CDP ELSE CDPU=CDP CDP=GTDR(1,1,LAY-1)*CDPU+GTDR(1,2,LAY-1)*CDS CDS=GTDR(2,1,LAY-1)*CDPU+GTDR(2,2,LAY-1)*CDS END IF CUP=GRDR(1,1,LAY)*CDP+GRDR(1,2,LAY)*CDS CUS=GRDR(2,1,LAY)*CDP+GRDR(2,2,LAY)*CDS END IF C IF(LAY.EQ.NL)THEN A1=0.5D0/CMP A2=0.5D0/CMS A3=1.D0/(CMP+CMS) C B1=ER(1,1)*ER(1,1)*CDP*CDP B2=ER(1,2)*ER(1,2)*CDS*CDS B3=2.D0*ER(1,1)*ER(1,2)*CDP*CDS DR1=B1*A1+B2*A2+B3*A3 B1=ER(2,1)*ER(2,1)*CDP*CDP B2=ER(2,2)*ER(2,2)*CDS*CDS B3=2.D0*ER(2,1)*ER(2,2)*CDP*CDS DR2=B1*A1+B2*A2+B3*A3 B1=ER(1,1)*ER(4,1)*CDP*CDP B2=ER(1,2)*ER(4,2)*CDS*CDS B3=(ER(1,1)*ER(4,2)+ER(1,2)*ER(4,1))*CDP*CDS DR14=B1*A1+B2*A2+B3*A3 B1=ER(2,1)*ER(3,1)*CDP*CDP B2=ER(2,2)*ER(3,2)*CDS*CDS B3=(ER(2,1)*ER(3,2)+ER(2,2)*ER(3,1))*CDP*CDS DR23=B1*A1+B2*A2+B3*A3 ELSE DDP=CDEXP(-2.D0*CMP*THK(LAY)) A1=0.5D0*(1.D0-DDP)/CMP DDS=CDEXP(-2.D0*CMS*THK(LAY)) A2=0.5D0*(1.D0-DDS)/CMS DDS=CDEXP(-(CMP+CMS)*THK(LAY)) A3=(1.D0-DDS)/(CMP+CMS) DDP=CDEXP(-CMP*THK(LAY)) A4=THK(LAY)*DDP DDS=CDEXP(-CMS*THK(LAY)) A5=(DDS-DDP)/(CMP-CMS) A6=THK(LAY)*DDS C B1=ER(1,1)*ER(1,1)*CDP*CDP+ER(1,3)*ER(1,3)*CUP*CUP B2=ER(1,2)*ER(1,2)*CDS*CDS+ER(1,4)*ER(1,4)*CUS*CUS B3=ER(1,1)*ER(1,2)*CDP*CDS+ER(1,3)*ER(1,4)*CUP*CUS B4=ER(1,1)*ER(1,3)*CDP*CUP B5=ER(1,1)*ER(1,4)*CDP*CUS+ER(1,2)*ER(1,3)*CUP*CDS B6=ER(1,2)*ER(1,4)*CDS*CUS DR1=B1*A1+B2*A2+2.D0*(B3*A3+B4*A4+B5*A5+B6*A6) C B1=ER(2,1)*ER(2,1)*CDP*CDP+ER(2,3)*ER(2,3)*CUP*CUP B2=ER(2,2)*ER(2,2)*CDS*CDS+ER(2,4)*ER(2,4)*CUS*CUS B3=ER(2,1)*ER(2,2)*CDP*CDS+ER(2,3)*ER(2,4)*CUP*CUS B4=ER(2,1)*ER(2,3)*CDP*CUP B5=ER(2,1)*ER(2,4)*CDP*CUS+ER(2,2)*ER(2,3)*CUP*CDS B6=ER(2,2)*ER(2,4)*CDS*CUS DR2=B1*A1+B2*A2+2.D0*(B3*A3+B4*A4+B5*A5+B6*A6) C B1=ER(1,1)*ER(4,1)*CDP*CDP+ER(1,3)*ER(4,3)*CUP*CUP B2=ER(1,2)*ER(4,2)*CDS*CDS+ER(1,4)*ER(4,4)*CUS*CUS B3=(ER(1,1)*ER(4,2)+ER(1,2)*ER(4,1))*CDP*CDS * +(ER(1,3)*ER(4,4)+ER(1,4)*ER(4,3))*CUP*CUS B4=(ER(1,1)*ER(4,3)+ER(1,3)*ER(4,1))*CDP*CUP B5=(ER(1,1)*ER(4,4)+ER(1,4)*ER(4,1))*CDP*CUS * +(ER(1,2)*ER(4,3)+ER(1,3)*ER(4,2))*CUP*CDS B6=(ER(1,2)*ER(4,4)+ER(1,4)*ER(4,2))*CDS*CUS DR14=B1*A1+B2*A2+B3*A3+B4*A4+B5*A5+B6*A6 C B1=ER(2,1)*ER(3,1)*CDP*CDP+ER(2,3)*ER(3,3)*CUP*CUP B2=ER(2,2)*ER(3,2)*CDS*CDS+ER(2,4)*ER(3,4)*CUS*CUS B3=(ER(2,1)*ER(3,2)+ER(2,2)*ER(3,1))*CDP*CDS * +(ER(2,3)*ER(3,4)+ER(2,4)*ER(3,3))*CUP*CUS B4=(ER(2,1)*ER(3,3)+ER(2,3)*ER(3,1))*CDP*CUP B5=(ER(2,1)*ER(3,4)+ER(2,4)*ER(3,1))*CDP*CUS * +(ER(2,2)*ER(3,3)+ER(2,3)*ER(3,2))*CUP*CDS B6=(ER(2,2)*ER(3,4)+ER(2,4)*ER(3,2))*CDS*CUS DR23=B1*A1+B2*A2+B3*A3+B4*A4+B5*A5+B6*A6 END IF C SM=DNS(LAY)*VS(LAY)*VS(LAY) RAM=DNS(LAY)*VP(LAY)*VP(LAY)-2.D0*SM A=RAM+2.D0*SM C E1R(IO,MD)=dreal(E1R(IO,MD)+DNS(LAY)*(DR1+DR2)) E2R =dreal(E2R +A*DR1+SM*DR2 ) E3R =dreal(E3R +RAM*(DR14-WN*RAM*DR1)/A-SM*WN*DR2-DR23) C 260 CONTINUE C E1R(IO,MD)=E1R(IO,MD)/2.D0 E2R =E2R /2.D0 UUR(IO,MD)=(E2R+E3R/(2.D0*WN))/(E1R(IO,MD)*CCR(IO,MD)) HoV(IO,MD)=dabs(VR1(1,MD)/VR2(1,MD)) c write(33,*)'io,md',io,md,VR1(1,MD),VR2(1,MD) C 210 CONTINUE c do MD=0,MDR(IO) VDIS=0.001d0 VSTR=0.001d0 do IZ=1,NZINT+1 if(dabs(VR1(IZ,MD)).gt.VDIS)VDIS=dabs(VR1(IZ,MD)) if(dabs(VR2(IZ,MD)).gt.VDIS)VDIS=dabs(VR2(IZ,MD)) if(dabs(VR3(IZ,MD)).gt.VSTR)VSTR=dabs(VR3(IZ,MD)) if(dabs(VR4(IZ,MD)).gt.VSTR)VSTR=dabs(VR4(IZ,MD)) end do do IZ=1,NZINT+1 VR1(IZ,MD)=VR1(IZ,MD)/VDIS VR2(IZ,MD)=VR2(IZ,MD)/VDIS VR3(IZ,MD)=VR3(IZ,MD)/VSTR VR4(IZ,MD)=VR4(IZ,MD)/VSTR end do end do c IF(IVCT.EQ.1)THEN do IZ=1,NZINT+1 Z=DZ*(IZ-1) write(26,2500)z/1000.d0,(VR1(IZ,MD)/2.d0+MD * ,VR2(IZ,MD)/2.d0+MD,MD=MDR(IO),0,-1) write(27,2500)z/1000.d0,(VR3(IZ,MD)/2.d0+MD * ,VR4(IZ,MD)/2.d0+MD,MD=MDR(IO),0,-1) end do END IF C 100 CONTINUE C IF(IVCT.EQ.1)THEN CLOSE(10) CLOSE(11) CLOSE(20) CLOSE(21) CLOSE(22) CLOSE(23) CLOSE(24) CLOSE(25) CLOSE(26) CLOSE(27) END IF c OPEN(32,FILE='grpl_dat',STATUS='unknown') OPEN(33,FILE='grpr_dat',STATUS='unknown') C DO 200 MD=0,MDL(NOM) DO 200 IO=IOM,NOM WRITE(32,2100)UUL(IO,MD) 200 CONTINUE C DO 300 MD=0,MDR(NOM) DO 300 IO=IOM,NOM WRITE(33,2100)UUR(IO,MD) 300 CONTINUE 2100 format(e15.7) C CLOSE(32) CLOSE(33) 2500 FORMAT(f8.4,500(',',e12.4)) c c *** Output of Dispersion Curves and Medium Response from Aki and Richards, c * MedRes-L (A & K, p302, (7.106)) -> |Gyy|=1/(4*k*c*U*I1)=1/(4*om*U*I1) c * MedRes-R (A & K, p303, (7.110)) -> |Gzz|=1/(4*k*c*U*I1)=1/(4*om*U*I1) c * "H over H Spectra" of the Rayleigh waves c OPEN(44,FILE='Ldisper.csv',STATUS='unknown') OPEN(45,FILE='Rdisper.csv',STATUS='unknown') OPEN(46,FILE='Lmedres.csv',STATUS='unknown') OPEN(47,FILE='Rmedres.csv',STATUS='unknown') OPEN(48,FILE='HoverV.csv',STATUS='unknown') WRITE(44,2600)(MD,MD,MD=0,MDL(NOM)) WRITE(45,2600)(MD,MD,MD=0,MDR(NOM)) WRITE(46,2650)(MD,MD=0,MDL(NOM)) WRITE(47,2650)(MD,MD=0,MDR(NOM)) WRITE(48,2650)(MD,MD=0,MDR(NOM)) do IO=IOM,NOM om=DOM*DFLOAT(IO) Hz=DOM*DFLOAT(IO)/6.283185307d0 WRITE(44,2620)Hz,(CCL(IO,MD),UUL(IO,MD),MD=0,MDL(NOM)) WRITE(45,2620)Hz,(CCR(IO,MD),UUR(IO,MD),MD=0,MDR(NOM)) do MD=0,MDL(NOM) if(io.eq.0)om=dom/1000.d0 res(MD)=0.d0 if (UUL(IO,MD).gt.0.00001d0)then res(MD)=0.25d0/(om*UUL(IO,MD)*E1L(IO,MD)) end if end do WRITE(46,2670)Hz,(res(MD),MD=0,MDL(NOM)) do MD=0,MDR(NOM) if(io.eq.0)om=dom/1000.d0 res(MD)=0.d0 if (UUR(IO,MD).gt.0.00001d0)then res(MD)=0.25d0/(om*UUR(IO,MD)*E1R(IO,MD)) end if end do WRITE(47,2670)Hz,(res(MD),MD=0,MDR(NOM)) WRITE(48,2670)Hz,(HoV(IO,MD),MD=0,MDR(NOM)) end do 2600 FORMAT('Freq (Hz)',500(',C',I3,',U',I3)) 2620 FORMAT(f8.4,500(',',e12.4)) 2650 FORMAT('Freq (Hz)',500(',M',I3)) 2670 FORMAT(f8.4,500(',',e12.4)) CLOSE(44) CLOSE(45) CLOSE(46) CLOSE(47) CLOSE(48) C RETURN END C C C *** SH GENERALIZED R/T MATIX *** C SUBROUTINE GNRTLV(OM,WN) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 GRU0,GTD,GRD,TD,TU,RD,RU COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL COMMON /GRTLV/GRU0,GTD(20),GRD(20) * /MRTLV/TD(20),TU(20),RD(20),RU(0:20) C C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM GRTLV: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C C *** MODIFIED R/T MATRIX *** C CALL MDRTLV(OM,WN) C C *** GENERALIZED R/T MATRIX *** C GRU0=RU(0) GTD(NL-1)=TD(NL-1) GRD(NL-1)=RD(NL-1) C DO 200 I=NL-2,1,-1 if(cdabs(1.D0-RU(I)*GRD(I+1)).lt.1.d-15)then write(10,*)'from GNRTLV:om,wn',om,wn write(10,*)' 1.D0-RU(I)*GRD(I+1) < 1.d-15 !!' write(10,*)1.D0-RU(I)*GRD(I+1) write(10,*)' Try a different number of Partitions !!' end if GTD(I)=TD(I)/(1.D0-RU(I)*GRD(I+1)) GRD(I)=RD(I)+TU(I)*GRD(I+1)*GTD(I) 200 CONTINUE C RETURN END C C C C *** SH MODIFIED R/T MATIX *** C SUBROUTINE MDRTLV(OM,WN) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 CMS,EA(2,2),EB(2,2),TD,TU,RD,RU,DDA,DUB,EAB COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL COMMON /MRTLV/TD(20),TU(20),RD(20),RU(0:20) C C *** CHECK FOR THE DIMENSION OF TD, TU, RD, AND RU *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM MDRTSH: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C C *** MODIFIED R/T MATRIX *** C WN2=WN*WN OM2=OM*OM C VS2=VS(1)*VS(1) WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS C SM=DNS(1)*VS2 EA(1,1)=(1.D0,0.D0) EA(1,2)=(1.D0,0.D0) EA(2,1)=-SM*CMS EA(2,2)= SM*CMS C DUB=CDEXP(-CMS*THK(1)) DDA=DUB C RU(0)=DUB C DO 100 I=1, NL-2 C VS2=VS(I+1)*VS(I+1) WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS C SM=DNS(I+1)*VS2 EB(1,1)=(1.D0,0.D0) EB(1,2)=(1.D0,0.D0) EB(2,1)=-SM*CMS EB(2,2)= SM*CMS C DUB=CDEXP(-CMS*THK(I+1)) C EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) TD(I)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB RU(I)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB RD(I)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB TU(I)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB C EA(1,1)=EB(1,1) EA(1,2)=EB(1,2) EA(2,1)=EB(2,1) EA(2,2)=EB(2,2) C DDA=DUB C 100 CONTINUE C VS2=VS(NL)*VS(NL) WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS C SM=DNS(NL)*VS2 EB(1,1)=(1.D0,0.D0) EB(1,2)=(1.D0,0.D0) EB(2,1)=-SM*CMS EB(2,2)= SM*CMS C EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2) C TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB C RETURN END C C C C *** P-SV GENERALIZED R/T MATIX *** C SUBROUTINE GNRTRY(OM,WN) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 GTD,GRD,TD,TU,RD,RU,E1(2,2),E2(2,2),DT COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL COMMON /MRTRY/TD(2,2,20),TU(2,2,20),RD(2,2,20),RU(2,2,20) * /GRTRY/GTD(2,2,20),GRD(2,2,20) C C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM GRTPRY: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C C *** MODIFIED R/T MATRIX *** C IF(NL.EQ.1)RETURN CALL MDRTRY(OM,WN) C C *** GENERALIZED R/T MATRIX *** C GTD(1,1,NL-1)=TD(1,1,NL-1) GTD(1,2,NL-1)=TD(1,2,NL-1) GTD(2,1,NL-1)=TD(2,1,NL-1) GTD(2,2,NL-1)=TD(2,2,NL-1) GRD(1,1,NL-1)=RD(1,1,NL-1) GRD(1,2,NL-1)=RD(1,2,NL-1) GRD(2,1,NL-1)=RD(2,1,NL-1) GRD(2,2,NL-1)=RD(2,2,NL-1) C DO 200 I=NL-2,1,-1 C E1(1,1)=1.D0-RU(1,1,I)*GRD(1,1,I+1)-RU(1,2,I)*GRD(2,1,I+1) E1(1,2)= -RU(1,1,I)*GRD(1,2,I+1)-RU(1,2,I)*GRD(2,2,I+1) E1(2,1)= -RU(2,1,I)*GRD(1,1,I+1)-RU(2,2,I)*GRD(2,1,I+1) E1(2,2)=1.D0-RU(2,1,I)*GRD(1,2,I+1)-RU(2,2,I)*GRD(2,2,I+1) DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1) E2(1,1)= E1(2,2)/DT E2(1,2)=-E1(1,2)/DT E2(2,1)=-E1(2,1)/DT E2(2,2)= E1(1,1)/DT GTD(1,1,I)=E2(1,1)*TD(1,1,I)+E2(1,2)*TD(2,1,I) GTD(1,2,I)=E2(1,1)*TD(1,2,I)+E2(1,2)*TD(2,2,I) GTD(2,1,I)=E2(2,1)*TD(1,1,I)+E2(2,2)*TD(2,1,I) GTD(2,2,I)=E2(2,1)*TD(1,2,I)+E2(2,2)*TD(2,2,I) C E1(1,1)=TU(1,1,I)*GRD(1,1,I+1)+TU(1,2,I)*GRD(2,1,I+1) E1(1,2)=TU(1,1,I)*GRD(1,2,I+1)+TU(1,2,I)*GRD(2,2,I+1) E1(2,1)=TU(2,1,I)*GRD(1,1,I+1)+TU(2,2,I)*GRD(2,1,I+1) E1(2,2)=TU(2,1,I)*GRD(1,2,I+1)+TU(2,2,I)*GRD(2,2,I+1) E2(1,1)=E1(1,1)*GTD(1,1,I)+E1(1,2)*GTD(2,1,I) E2(1,2)=E1(1,1)*GTD(1,2,I)+E1(1,2)*GTD(2,2,I) E2(2,1)=E1(2,1)*GTD(1,1,I)+E1(2,2)*GTD(2,1,I) E2(2,2)=E1(2,1)*GTD(1,2,I)+E1(2,2)*GTD(2,2,I) GRD(1,1,I)=RD(1,1,I)+E2(1,1) GRD(1,2,I)=RD(1,2,I)+E2(1,2) GRD(2,1,I)=RD(2,1,I)+E2(2,1) GRD(2,2,I)=RD(2,2,I)+E2(2,2) C 200 CONTINUE C RETURN END C C C C *** P-SV MODIFIED R/T MATIX *** C SUBROUTINE MDRTRY(OM,WN) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 CMS,CMP,EA(4,4),EB(4,4),DDA(2,2),DUB(2,2) COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL C C *** CHECK FOR THE DIMENSION OF COMMON *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM MDRTPSV: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C C *** MODIFIED R/T MATRIX *** C CALL MTXRY(OM,WN,1,CMS,CMP,EA) CALL DDURY(THK(1),CMS,CMP,DDA) C DO 100 I=1, NL-1 C CALL MTXRY(OM,WN,I+1,CMS,CMP,EB) CALL DDURY(THK(I+1),CMS,CMP,DUB) C CALL MTRRY(EA,EB,DDA,DUB,I) C DO 110 J=1,4 DO 120 K=1,4 EA(J,K)=EB(J,K) 120 CONTINUE 110 CONTINUE C DDA(1,1)=DUB(1,1) DDA(1,2)=DUB(1,2) DDA(2,1)=DUB(2,1) DDA(2,2)=DUB(2,2) C 100 CONTINUE C RETURN END C C C C *** P-SV MODIFIED R/T MATRIX FOR Ith LAYER *** C SUBROUTINE MTRRY(EA,EB,DDA,DUB,L) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 EA(4,4),EB(4,4),DDA(2,2),DUB(2,2),TD,TU,RD,RU * ,E1(4,4),E2(4,4),E3(4,4) COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL * /MRTRY/TD(2,2,20),TU(2,2,20),RD(2,2,20),RU(2,2,20) C C *** CHECK FOR THE DIMENSION OF COMMON *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM MTRRY: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C C *** MATRIX *** C DO 100 I=1,2 DO 110 J=1,2 E1(I,J)= EB(I,J) E1(I,J+2)=-EA(I,J+2) E1(I+2,J)= EB(I+2,J) E1(I+2,J+2)=-EA(I+2,J+2) E2(I,J)= EA(I,J) E2(I,J+2)=-EB(I,J+2) E2(I+2,J)= EA(I+2,J) E2(I+2,J+2)=-EB(I+2,J+2) 110 CONTINUE 100 CONTINUE C CALL DCINVA(E1) C DO 200 I=1,4 DO 210 J=1,4 E3(I,J)=(0.D0,0.D0) DO 220 K=1,4 E3(I,J)=E3(I,J)+E1(I,K)*E2(K,J) 220 CONTINUE 210 CONTINUE 200 CONTINUE C TD(1,1,L)=E3(1,1)*DDA(1,1)+E3(1,2)*DDA(2,1) TD(2,1,L)=E3(2,1)*DDA(1,1)+E3(2,2)*DDA(2,1) TD(1,2,L)=E3(1,1)*DDA(1,2)+E3(1,2)*DDA(2,2) TD(2,2,L)=E3(2,1)*DDA(1,2)+E3(2,2)*DDA(2,2) RD(1,1,L)=E3(3,1)*DDA(1,1)+E3(3,2)*DDA(2,1) RD(2,1,L)=E3(4,1)*DDA(1,1)+E3(4,2)*DDA(2,1) RD(1,2,L)=E3(3,1)*DDA(1,2)+E3(3,2)*DDA(2,2) RD(2,2,L)=E3(4,1)*DDA(1,2)+E3(4,2)*DDA(2,2) C IF(L.EQ.NL)RETURN C RU(1,1,L)=E3(1,3)*DUB(1,1)+E3(1,4)*DUB(2,1) RU(2,1,L)=E3(2,3)*DUB(1,1)+E3(2,4)*DUB(2,1) RU(1,2,L)=E3(1,3)*DUB(1,2)+E3(1,4)*DUB(2,2) RU(2,2,L)=E3(2,3)*DUB(1,2)+E3(2,4)*DUB(2,2) TU(1,1,L)=E3(3,3)*DUB(1,1)+E3(3,4)*DUB(2,1) TU(2,1,L)=E3(4,3)*DUB(1,1)+E3(4,4)*DUB(2,1) TU(1,2,L)=E3(3,3)*DUB(1,2)+E3(3,4)*DUB(2,2) TU(2,2,L)=E3(4,3)*DUB(1,2)+E3(4,4)*DUB(2,2) C RETURN END C C C *** Inverse of 4x4 Matrix *** C SUBROUTINE DCINVA(E) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 E(4,4),ee(4,4),aa * ,a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3 C do 100 i=1,4 do 100 j=1,4 ee(i,j)=e(i,j) 100 continue C DO 200 I=1,4 C IF(I.EQ.1)THEN J=2 K=3 M=4 ELSE IF(I.EQ.2)THEN J=3 K=4 M=1 ELSE IF(I.EQ.3)THEN J=4 K=1 M=2 ELSE J=1 K=2 M=3 END IF END IF END IF C A1=EE(I,4)*EE(J,1)-EE(I,1)*EE(J,4) A2=EE(I,4)*EE(J,2)-EE(I,2)*EE(J,4) A3=EE(I,4)*EE(J,3)-EE(I,3)*EE(J,4) B1=EE(I,3)*EE(J,2)-EE(I,2)*EE(J,3) B2=EE(I,1)*EE(J,3)-EE(I,3)*EE(J,1) B3=EE(I,2)*EE(J,1)-EE(I,1)*EE(J,2) C1=EE(K,3)*EE(M,2)-EE(K,2)*EE(M,3) C2=EE(K,1)*EE(M,3)-EE(K,3)*EE(M,1) C3=EE(K,2)*EE(M,1)-EE(K,1)*EE(M,2) D1=EE(K,4)*EE(M,1)-EE(K,1)*EE(M,4) D2=EE(K,4)*EE(M,2)-EE(K,2)*EE(M,4) D3=EE(K,4)*EE(M,3)-EE(K,3)*EE(M,4) C AA=A1*C1+A2*C2+A3*C3+B1*D1+B2*D2+B3*D3 C E(1,I)=(-EE(J,2)*D3+EE(J,3)*D2-EE(J,4)*C1)/AA E(2,I)=( EE(J,1)*D3-EE(J,3)*D1-EE(J,4)*C2)/AA E(3,I)=(-EE(J,1)*D2+EE(J,2)*D1-EE(J,4)*C3)/AA E(4,I)=( EE(J,1)*C1+EE(J,2)*C2+EE(J,3)*C3)/AA C 200 continue C return end c C C *** P-SV ELEMENT MATIX *** C SUBROUTINE MTXRY(OM,WN,L,CMS,CMP,EE) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 CMS,CMP,EE(4,4),CMSG,CMPG COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL C C *** CHECK FOR THE DIMENSION OF COMMOM *** C NC=20 IF(NL.GT.NC)THEN WRITE(*,*)' FROM MTXRY: NL.GT.20 !!!!!!!!!!!!!!!' STOP END IF C C *** ELEMENT MATRIX *** C VS2=VS(L)*VS(L) VP2=VP(L)*VP(L) SM=DNS(L)*VS2 C IF(OM.EQ.0.D0)THEN C X0=(1.D0+VS2/VP2)/(1.D0-VS2/VP2) AA=WN*SM C EE(1,1)=(1.D0,0.D0) EE(1,2)=(1.D0,0.D0) EE(1,3)=(1.D0,0.D0) EE(1,4)=(1.D0,0.D0) EE(2,1)=-(X0-1.D0) EE(2,2)=(1.D0,0.D0) EE(2,3)= X0-1.D0 EE(2,4)=(-1.D0,0.D0) EE(3,1)= (X0-3.D0)*AA EE(3,2)= -2.D0*AA EE(3,3)=-(X0-3.D0)*AA EE(3,4)= 2.D0*AA EE(4,1)= (X0-1.D0)*AA EE(4,2)= -2.D0*AA EE(4,3)= (X0-1.D0)*AA EE(4,4)=-2.D0*AA C ELSE C WN2=WN*WN OM2=OM*OM C WS2=OM2/VS2 CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) IF(DREAL(CMS).LT.0.D0)CMS=-CMS CMSG=CMS/WN C WP2=OM2/VP2 CMP=CDSQRT(DCMPLX(WN2-WP2,0.D0)) IF(DREAL(CMP).LT.0.D0)CMP=-CMP CMPG=CMP/WN C X0=2.D0*WN-WS2/WN C EE(1,1)=(-1.D0,0.D0) EE(1,2)= CMSG EE(1,3)=(-1.D0,0.D0) EE(1,4)= CMSG EE(2,1)=-CMPG EE(2,2)=(1.D0,0.D0) EE(2,3)= CMPG EE(2,4)=(-1.D0,0.D0) EE(3,1)= 2.D0*SM*CMP EE(3,2)=-SM*X0 EE(3,3)=-2.D0*SM*CMP EE(3,4)= SM*X0 EE(4,1)= SM*X0 EE(4,2)=-2.D0*SM*CMS EE(4,3)= SM*X0 EE(4,4)=-2.D0*SM*CMS C END IF C RETURN END C C C *** P-SV UP/DOWN-GOING MATRIX *** C SUBROUTINE DDURY(THK,CMS,CMP,DU) C IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 DU(2,2),CMS,CMP C DU(1,1)= CDEXP(-CMP*THK) DU(1,2)= (0.D0,0.D0) DU(2,1)= (0.D0,0.D0) DU(2,2)= CDEXP(-CMS*THK) C RETURN END