C
C PROGRAM SEIS C ************** C C PROGRAM SEIS IS DESIGNED FOR A TWO POINT RAY TRACING AND C RAY SYNTHETIC SEISMOGRAM COMPUTATION IN A 2-D LATERALLY C INHOMOGENEOUS MEDIA WITH CURVED INTERFACES AND BLOCK STRUCTURES. C COMPLETE WAVE FIELD (ALL THREE COMPONENTS) GENERATED BY C SIMPLE 3-D POINT SOURCES IS COMPUTED. C C **************************************************************** C C CHARACTER*80 MTEXT,FILEIN,FILEOU,FILE1,FILE2 INTEGER CODE DIMENSION VEL(6),Y(4) COMMON/DIST/DST(100),NDST,IDIST,ASTART,ASTEP,AFIN,REPS,REPS1,RSTEP COMMON/COD/CODE(50),KREF,KC COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET COMMON/DEN/RHO1(19),RHO2(19),NRHO COMMON/SOUR/ROS,VPS,VSS COMMON/SH/KSH,NSH common/vsp/xvsp,icod,ivsp C C *** LTAPE=0 LIN=5 LOU=6 LU1=1 LU2=2 FILEIN='seis.dat' FILEOU='seis.out' FILE1='lu1.dat' FILE2='lu2.dat' WRITE(*,'(2A)') ' (SEIS) SPECIFY NAMES OF INPUT AND OUTPUT', 1' FILES LIN, LOU, LU1, LU2: ' READ(*,*) FILEIN,FILEOU,FILE1,FILE2 IF(FILE1.EQ.' ') LU1=0 IF(FILE2.EQ.' ') LU2=0 OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') IF(LU1.NE.0)OPEN(LU1,FILE=FILE1,FORM='FORMATTED') IF(LU2.NE.0)OPEN(LU2,FILE=FILE2,FORM='FORMATTED') C C *** C WRITE(LOU,777) 777 FORMAT(///,'***********************' 1,//,' PROGRAM S E I S ',//, 2'***********************',//) C C 1 READ(LIN,*)MTEXT WRITE(LOU,110)MTEXT READ(LIN,*)MPRINT,MM WRITE(LOU,100)MPRINT,MM C C SPECIFICATION OF THE MODEL C CALL MODEL(LIN,LOU,MM,MTEXT) IF(MM.LT.0)IND=19 IF(MM.LT.0)WRITE(LOU,112)IND IF(MM.LT.0)GO TO 99 C C SPECIFICATION OF THE SYNTHETIC SEISMOGRAMS C 2 CONTINUE ICONT=0 READ(LIN,*)ICONT,MEP,MOUT,MDIM,METHOD,IBP,IBS,iup,ius,IDP,IDS, 1IREAD,MREG,ITMAX,NLAY,LU1,LU2,LTAPE,KSH,Iloc IF(ICONT.EQ.0)GO TO 99 LTAPE=0 WRITE(LOU,102)ICONT,MEP,MOUT,MDIM,METHOD,IBP,IBS,iup,ius,IDP,IDS, 1IREAD,MREG,ITMAX,NLAY,LU1,LU2,LTAPE,KSH,Iloc IF(ICONT.EQ.(-1))GO TO 1 C C *** C *** C IF(LU1.NE.0)REWIND LU1 IF(LU2.NE.0)REWIND LU2 IF(NRHO.EQ.1)MREG=1 IF(ITMAX.EQ.0)ITMAX=20 IREC=1 IMET=0 ivsp=0 itpr=iloc IF(ITPR.EQ.0)ITPR=3 IF(ITPR.GT.100)IREC=ITPR-100 if(iloc.eq.1)ivsp=1 if(iloc.eq.1)itpr=22 if(ivsp.eq.1)mreg=1 if(ivsp.eq.1)irec=0 IF(IREC.NE.1)MREG=1 NLAYER=NINT-NLAY NDR=8 NDER=0 IF(MDIM.NE.0)NDER=1 IF(MDIM.NE.0)NDR=11 C IF(MEP.GT.0)GO TO 3 NDST=-MEP READ(LIN,*)(DST(I),I=1,NDST),xvsp WRITE(LOU,104)(DST(I),I=1,NDST),xvsp if(ndst.eq.1)rstep=1. if(ndst.eq.1)dst(2)=dst(1)+1. if(ndst.eq.1)go to 4 RSTEP=(DST(NDST)-DST(1))/FLOAT(NDST-1) GO TO 4 3 NDST=MEP READ(LIN,*)RMIN,RSTEP,xvsp WRITE(LOU,104)RMIN,RSTEP,xvsp DO 13 I=1,MEP 13 DST(I)=RMIN+(I-1)*RSTEP C 4 READ(LIN,*)XSOUR,ZSOUR,TSOUR,REPS,reps1,BRD(1),BRD(2) WRITE(LOU,104)XSOUR,ZSOUR,TSOUR,REPS,reps1,BRD(1),BRD(2) READ(LIN,*)DT,AMIN1,ASTEP1,AMAX1,AMIN2,ASTEP2,AMAX2,AC WRITE(LOU,104)DT,AMIN1,ASTEP1,AMAX1,AMIN2,ASTEP2,AMAX2,AC IF(DT.LT.0.0000001)DT=1. IF(AC.LT.0.0000001)AC=0.0001 TMAX=10000. IF(REPS.LT..00001)REPS=.05 IF(REPS1.LT..00001)REPS1=.05*REPS IND=-1 NC=NPNT(1) BL=X1(1,1) BR=X1(NC,1) BA=BRD(1) BB=BRD(2) IF(ABS(BB-BA).GT..00001)GO TO 5 BRD(1)=BL BRD(2)=BR 5 IF(BA.LT.BL)BRD(1)=BL IF(BB.GT.BR)BRD(2)=BR c c determination of the position of the source in the model c if(xsour.lt.brd(1).or.xsour.gt.brd(2))go to 205 intr=1 201 nl=npnt(intr)-1 do 202 i=1,nl j=i if(xsour.lt.x1(i+1,intr))go to 203 202 continue 203 a2=a1(j,intr) b2=b1(j,intr) c2=c1(j,intr) d2=d1(j,intr) x2=x1(j,intr) au=xsour-x2 zint=((d2*au+c2)*au+b2)*au+a2 if(abs(zsour).gt..0001)go to 204 isour=1 zsour=zint go to 206 204 if(zsour.lt.zint.and.intr.eq.1)go to 205 if(zsour.lt.zint)go to 206 if(abs(zsour-zint).lt..000001.and.intr.eq.nint)go to 206 if(intr.eq.nint)go to 205 isour=intr intr=intr+1 go to 201 205 ind=50 write(6,112)ind go to 99 206 lay=isour int1=isour ITYPE=-1 Y(1)=XSOUR Y(2)=ZSOUR Y(3)=AMIN1 Y(4)=0. CALL VELOC(Y,VEL) VSS=VEL(1) ITYPE=1 CALL VELOC(Y,VEL) VPS=VEL(1) ROS=RHO1(Isour)+RHO2(Isour)*VPS C C C GENERATE FILE LU2 FOR PLOTTING OF SYNTHETIC SEISMOGRAMS C IF(LU2.EQ.0)GO TO 25 IF(LTAPE.EQ.0)WRITE(LU2,110)MTEXT IF(LTAPE.EQ.1)WRITE(LU2)MTEXT IF(LTAPE.EQ.0)WRITE(LU2,100)NDST,KSH,itpr IF(LTAPE.EQ.1)WRITE(LU2)NDST,KSH,itpr IF(LTAPE.EQ.0)WRITE(LU2,104)XSOUR,ZSOUR,TSOUR,RSTEP,ROS,VPS,VSS IF(LTAPE.EQ.1)WRITE(LU2)XSOUR,ZSOUR,TSOUR,RSTEP,ROS,VPS,VSS IF(LTAPE.EQ.0)WRITE(LU2,104)(DST(I),I=1,NDST) IF(LTAPE.EQ.1)WRITE(LU2)(DST(I),I=1,NDST) 25 CONTINUE C C LOOP FOR ELEMENTARY WAVES C NAWX=0 IF(MOUT.EQ.0)WRITE(LOU,105) 20 CALL GENER(NAWX,IDP,IDS,IBP,IBS,IUP,IUS,IREAD,ISOUR,IREC,NLAYER, 1NLAY,IFIN,NCODE,lin,lou) NAWX=1 IF(IFIN.EQ.1)GO TO 49 IF(MOUT.NE.0)WRITE(LOU,107) WRITE(LOU,106)NCODE,KC,KREF,(CODE(I),I=1,KREF) icod=0 if(ivsp.eq.0)go to 35 do 34 i=1,kref icod=kref-i+1 if(icod.eq.1)go to 34 ic1=code(icod) ic2=code(icod-1) if((ic1-ic2).eq.0)go to 35 if((ic1*ic2).lt.0)go to 35 34 continue 35 continue NSH=0 DO 301 I=1,KREF IF(CODE(I).GT.0)GO TO 302 301 CONTINUE NSH=KSH 302 CONTINUE IF(KSH.EQ.1.AND.NSH.EQ.0)GO TO 20 IF(MOUT.NE.0)WRITE(LOU,108) IF(IFIN.EQ.(-1))GO TO 20 C C GENERATE FILE LU1 FOR PLOTTING OF RAY DIAGRAMS,TRAVEL TIMES AND C AMPLITUDES C IF(LU1.EQ.0)GO TO 21 if(ltape.eq.0)WRITE(LU1,100)ICONT,itpr if(ltape.eq.1)WRITE(LU1)ICONT,itpr if(ltape.eq.0)WRITE(LU1,100)NINT,(NPNT(I),I=1,NINT) if(ltape.eq.1)WRITE(LU1)NINT,(NPNT(I),I=1,NINT) DO 22 IIJ=1,NINT NI=NPNT(IIJ)-1 if(ltape.eq.0) 1WRITE(LU1,101)(A1(J,IIJ),B1(J,IIJ),C1(J,IIJ),D1(J,IIJ),X1(J,IIJ), 2III(J,IIJ),J=1,NI) if(ltape.eq.1) 1WRITE(LU1)(A1(J,IIJ),B1(J,IIJ),C1(J,IIJ),D1(J,IIJ),X1(J,IIJ), 2III(J,IIJ),J=1,NI) 22 continue if(ltape.eq.0)WRITE(LU1,104)XSOUR,ZSOUR,ROS,VPS,VSS if(ltape.eq.1)WRITE(LU1)XSOUR,ZSOUR,ROS,VPS,VSS 21 CONTINUE C C ASTART=AMIN1 ASTEP=ASTEP1 AFIN=AMAX1 IF(KC.LE.0)ASTART=AMIN2 IF(KC.LE.0)ASTEP=ASTEP2 IF(KC.LE.0)AFIN=AMAX2 CALL TRAMP(XSOUR,ZSOUR,TSOUR,DT,AC,TMAX,ITMAX,MOUT, 1MDIM,NCODE,lou,LU1,LU2,METHOD,ISOUR,LTAPE,ITPR) IF(IND.EQ.20.OR.IND.EQ.21)WRITE(LOU,111)IND,LAY IF(IND.EQ.20.OR.IND.EQ.21)GO TO 99 IF(IND.EQ.14)WRITE(LOU,112)IND GO TO 20 C C END OF LOOP FOR ELEMENTARY WAVES C C 49 IF(LU1.EQ.0)GO TO 50 JJ=0 if(ltape.eq.0)WRITE(LU1,100)JJ,JJ,JJ if(ltape.eq.1)WRITE(LU1)JJ,JJ,JJ REWIND LU1 50 IF(LU2.EQ.0)GO TO 2 REWIND LU2 C C *** C *** C c GO TO 2 C 100 FORMAT(26I3) 101 format(5e15.5,i5) 102 FORMAT(1H1,////,2X,26I3) 104 FORMAT(8F10.5) 105 FORMAT(//2X,'LISTING OF WAVE CODES'//2X,'INT.CODE',5X,'E X T E R N 1 A L C O D E') 106 FORMAT(I10,5X,31I3/18X,30I3) 107 FORMAT(//2X,'INT.CODE',5X,'E X T E R N A L C O D E') 108 FORMAT(//) 110 FORMAT(A) 111 FORMAT(/2X,'IND=',I5,2X,'LAY=',I5/) 112 format(/2x,'IND=',i5/) C C *** 99 CONTINUE C *** C STOP END C C *************************************************************** C SUBROUTINE TRAMP(XSOUR,ZSOUR,TSOUR,DT,AC,TMAX,ITMAX,MOUT, 1MDIM,NCODE,lou,LU1,LU2,METHOD,ISOUR,LTAPE,ITPR) C C TWO-POINT RAY TRACING C INTEGER CODE DIMENSION TIME(400),XCOOR(400),AMPX(400),AMPY(400),AMPZ(400), 1INDI(400),TAS(400),ANG(400),PHAX(400),PHAY(400),PHAZ(400),JC(50), 2QP(15) COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/DIST/DST(100),NDST,IDIST,ASTART,ASTEP,AFIN,REPS,REPS1,RSTEP COMMON/COD/CODE(50),KREF,KC COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),XX(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19), 1NQP(19),NQS(19),NABS COMMON/SH/KSH,NSH common/vsp/xvsp,icod,ivsp C TAST=0. PI=3.14159265 C AA=ASTART-ASTEP INDEX=0 INUM=0 IK1=0 ICLS=0 ICR=0 ISUC=0 INDS=ISOUR C C LOOP IN RAY PARAMETERS AA, FROM ASTART TO AFIN, WITH THE STEP C ASTEP C 1 AA=AA+ASTEP IF(ASTEP.GT.0..AND.AA.GT.AFIN)GO TO 99 IF(ASTEP.LT.0..AND.AA.LT.AFIN)GO TO 99 IND=INDS CALL RAY2(XSOUR,ZSOUR,TSOUR,AA,X,Z,T,FI,TMAX,DT,AC) IF(IND.EQ.14.OR.IND.EQ.20.OR.IND.EQ.21)RETURN if(itpr.eq.22)x=z IF(IND.EQ.ITPR)XAX=X IF(IND.EQ.ITPR)PNEW=AA IF(IND.EQ.9)XCR=X IF(MOUT.EQ.2)WRITE(LOU,100)IND,IND1,X,T,AA IF(INUM.NE.0)GO TO 2 AOLD=AA IOLD=IND IOLD1=IND1 XOLD=X TOLD=T INUM=1 GO TO 1 C C PARAMETERS FOR THE PRECEDING RAY: AA=AOLD, X=XOLD, IND=IOLD C PARAMETERS FOR THE NEW RAY: AA=ANEW, X=XNEW, IND=INEW C 2 INEW=IND INEW1=IND1 ANEW=AA XNEW=X TNEW=T IF(INEW.EQ.ITPR.AND.IOLD.EQ.ITPR)GO TO 21 IF(INEW.EQ.ITPR)GO TO 50 IF(IOLD.EQ.ITPR)GO TO 55 IF(INEW.EQ.9.AND.IOLD.NE.9.AND.IOLD.NE.8)GO TO 30 IF(INEW.NE.9.AND.INEW.NE.8.AND.IOLD.EQ.9)GO TO 35 GO TO 3 21 IF(IOLD1.NE.INEW1)IK1=2 IF(IK1.EQ.2)GO TO 55 GO TO 40 C C NO ITERATIONS, TAKE A NEW RAY IN THE LOOP C 3 IOLD=INEW IOLD1=INEW1 XOLD=XNEW AOLD=ANEW TOLD=TNEW GO TO 1 C C REGULAR CASE: IOLD=3, INEW=3 C 40 XXNEW=XNEW XXOLD=XOLD AANEW=ANEW AAOLD=AOLD TTNEW=TNEW TTOLD=TOLD IBNEW=0 IBOLD=0 41 IF(XXNEW.GT.XXOLD.AND.DST(2).GT.DST(1))GO TO 46 IF(XXNEW.LT.XXOLD.AND.DST(2).LT.DST(1))GO TO 46 C C REGULAR CASE: XXNEW.LE.XXOLD, ITREND=-1 (REVERSE BRANCH) C ITREND=-1 DO 42 I=1,NDST NDD=NDST-I+1 DD=DST(NDD) DX=XXOLD IF(IBOLD.EQ.1)DX=DX+REPS IF(DD.GE.DX)GO TO 42 DX=XXNEW IF(IBNEW.EQ.1)DX=DX-REPS IF(DD.LT.DX.AND.IK1.NE.0)GO TO 6 IF(DD.LT.DX)GO TO 3 II=NDD GO TO 43 42 CONTINUE IF(IK1.NE.0)GO TO 6 GO TO 3 C C REGULAR CASE: XXNEW.GT.XXOLD, ITREND=1 (REGULAR BRANCH) C 46 ITREND=1 DO 44 I=1,NDST DD=DST(I) DX=XXOLD IF(IBOLD.EQ.1)DX=DX-REPS IF(DD.LE.DX)GO TO 44 DX=XXNEW IF(IBNEW.EQ.1)DX=DX+REPS IF(DD.GT.DX.AND.IK1.NE.0)GO TO 6 IF(DD.GT.DX)GO TO 3 II=I GO TO 43 44 CONTINUE IF(IK1.NE.0)GO TO 6 GO TO 3 43 IF(METHOD.NE.3)GO TO 45 49 NPT=NPNT(1)-1 DO 47 I=1,NPT J=I IF(DD.LT.XX(I+1,1))GO TO 48 47 CONTINUE 48 AU=DD-XX(J,1) DDZ=((D1(J,1)*AU+C1(J,1))*AU+B1(J,1))*AU+A1(J,1) IF(ISUC.EQ.1)GO TO 60 45 P1=AAOLD P2=AANEW X1=XXOLD X2=XXNEW T1=TTOLD T2=TTNEW C C I T E R A T I O N S C 60 ITER=0 IF(ABS(X-DD).LT.REPS)GO TO 65 ISIGN=1 IPR1=0 IPR2=0 ISUC=0 GO TO 61 68 XAX=X PAX=PNEW 61 ITER=ITER+1 IF(ITER.GT.ITMAX)GO TO 80 IF(METHOD.EQ.3.AND.IND.EQ.itpr)GO TO 66 ISIGN=-ISIGN PNEW=0.5*(P1+P2) IF(ISIGN.LT.0.OR.METHOD.GT.1)GO TO 69 AUX=X2-X1 IF(ABS(AUX).LT..000001)GO TO 69 AUX=P1+(DD-X1)*(P2-P1)/AUX IF((AUX-P1)*(AUX-P2).GT.0.)GO TO 69 PNEW=AUX GO TO 69 66 if(itpr.eq.22)then dx=0. dz=dd-x end if if(itpr.ne.22)then DZ=DDZ-Z DX=DD-X end if NNN=N IREF1=IREF 64 N=KINT(IREF1) IF(N.LT.0)IREF1=IREF1-1 IF(N.LT.0)GO TO 64 V=AY(6,N) TX=AY(4,N)*V TZ=AY(5,N)*V DN=DX*TZ-DZ*TX DSS=-(DX*TX+DZ*TZ) N=0 KMAH=0 CALL JPAR(qp,afact,KMAH) N=NNN QS=Qp(3)+V*qP(4)*DSS IF(ABS(QS).LT..0001)then pnew=0.5*(p1+p2) go to 69 end if V0=AY(6,1) SN=(DN*V0)/QS PNEW=PNEW-SN IF(p1.lt.p2.and.(pnew.lt.p1.or.pnew.gt.p2))then pnew=0.5*(p1+p2) go to 69 end if IF(p1.gt.p2.and.(pnew.gt.p1.or.pnew.lt.p2))then pnew=0.5*(p1+p2) go to 69 end if 69 IND=INDS CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC) IF(IND.EQ.20.OR.IND.EQ.21)RETURN if(itpr.eq.22)x=z IF(MOUT.EQ.2)WRITE(LOU,101)IND,IND1,ITER,DD,X,T,PNEW IF(ICLS.NE.0)GO TO 70 IF(IND.NE.ITPR)P2=PNEW IF(IND.NE.ITPR)GO TO 61 IF(ABS(X-DD).LT.REPS)GO TO 65 IF(X1.LT.X2.AND.DD.GT.X)GO TO 63 IF(X1.GT.X2.AND.DD.LT.X)GO TO 63 62 IF(ABS(X-X1).LT..000001)GO TO 67 P2=PNEW X2=X T2=T IPR2=1 GO TO 68 63 IF(ABS(X-X2).LT..000001)GO TO 67 P1=PNEW X1=X T1=T IPR1=1 GO TO 68 67 IF(ABS(PNEW-PAX).GT..000001)ITER=ITMAX AX1=X1-DD AX2=X2-DD IF((IPR1*IPR2).EQ.0)ITER=ITMAX X=X1 PNEW=P1 IF(ABS(AX1).GT.ABS(AX2))PNEW=P2 IF(ABS(AX1).GT.ABS(AX2))X=X2 IF(ITER.EQ.ITMAX)GO TO 61 ICLS=1 GO TO 69 70 ICLS=0 GO TO 65 C C SUCCESSFUL ITERATIONS C 65 INDEX=INDEX+1 XAX=X IF(LU1.EQ.0)GO TO 800 TIME(INDEX)=T XCOOR(INDEX)=X INDI(INDEX)=IND ANG(INDEX)=PNEW 800 CONTINUE IF(ITPR.NE.22)AUX=AY(4,N) IF(ITPR.EQ.22)AUX=AY(5,N) TAUX=T T=Taux+AUX*(DD-X) IF(MOUT.EQ.2)WRITE(LOU,102)IND,IND1,ITER,II,INDEX,DD,X,T,PNEW IF(MOUT.EQ.2.AND.MDIM.EQ.0)WRITE(LOU,114) if(itpr.eq.22)z=x if(itpr.eq.22)x=xvsp IF(MOUT.EQ.1.AND.MDIM.EQ.0)WRITE(LOU,113)DD,X,Z,T,PNEW,IND, 1IND1,ITER,II if(itpr.eq.22)x=z if(lu1.eq.0)go to 802 if(ltape.eq.0)WRITE(LU1,108)N,IND if(ltape.eq.0)WRITE(LU1,109)(AY(2,I),AY(3,I),I=1,N) if(ltape.eq.1)WRITE(LU1)N,IND if(ltape.eq.1)WRITE(LU1)(AY(2,I),AY(3,I),I=1,N) 802 continue PH=0. IF(MDIM.EQ.0)GO TO 80 IF(MDIM.EQ.1)SPR=1. IF(MDIM.EQ.1)GO TO 10 NNN=N KMAH=0 N=0 CALL JPAR(QP,AFACT,KMAH) PH=PH-1.57079632*KMAH SPR=ABS(QP(3)/AY(6,1)) N=NNN 10 CALL AMPL(SPR,AX,AAY,AZ,PHX,PHY,PHZ) PHX=PHX+PH PHY=PHY+PH PHZ=PHZ+PH TAST=QP(5) IF(MDIM.LT.3)GO TO 12 AUX=QP(6)/AY(6,1) AUX=SQRT(1./ABS(AUX)) SPR=SPR/AUX AX=AX*AUX AAY=AAY*AUX AZ=AZ*AUX 12 CONTINUE IF(PHX.GE.(-PI))GO TO 13 PHX=PHX+2.*PI GO TO 12 13 IF(PHX.LE.PI)GO TO 14 PHX=PHX-2.*PI GO TO 13 14 IF(PHZ.GE.(-PI))GO TO 15 PHZ=PHZ+2.*PI GO TO 14 15 IF(PHZ.LE.PI)GO TO 16 PHZ=PHZ-2.*PI GO TO 15 16 CONTINUE IF(LU1.EQ.0)GO TO 801 TAS(INDEX)=TAST PHAX(INDEX)=PHX PHAY(INDEX)=PHY PHAZ(INDEX)=PHZ AMPZ(INDEX)=AZ AMPY(INDEX)=AAY AMPX(INDEX)=AX 801 CONTINUE if(itpr.eq.22)z=x if(itpr.eq.22)x=xvsp IF(MOUT.EQ.1)WRITE(LOU,105)DD,X,Z,T,TAST,PNEW,SPR, 1IND,IND1,ITER,II,KMAH,AX,PHX,AAY,PHY,AZ,PHZ if(itpr.eq.22)x=z IF(MOUT.EQ.2)WRITE(LOU,110)AX,PHX,AAY,PHY,AZ,PHZ,SPR,TAST,KMAH NCD=NCODE IF(CODE(1).LT.0)NCD=-NCODE C IF(LU2.NE.0.AND.LTAPE.EQ.0) 1WRITE(LU2,116)NCD,II,T,AX,AAY,AZ,PHX,PHY,PHZ,PNEW,TAST IF(LU2.NE.0.AND.LTAPE.EQ.1) 1WRITE(LU2)NCD,II,T,AX,AAY,AZ,PHX,PHY,PHZ,PNEW,TAST C C C 80 P1=PNEW X1=X T1=TAUX IF(ITER.GT.ITMAX)P1=AAOLD IF(ITER.GT.ITMAX)X1=XXOLD IF(ITER.GT.ITMAX)T1=TTOLD P2=AANEW X2=XXNEW T2=TTNEW IF(ITREND.EQ.(-1))GO TO 85 II=II+1 IF(II.GT.NDST.AND.IK1.NE.0)GO TO 6 IF(II.GT.NDST)GO TO 3 DD=DST(II) DX=XXNEW IF(IBNEW.EQ.1)DX=DX+REPS IF(DD.GT.DX.AND.IK1.NE.0)GO TO 6 IF(DD.GT.DX)GO TO 3 ISUC=1 IF(METHOD.EQ.3)GO TO 49 GO TO 60 85 II=II-1 IF(II.LT.1.AND.IK1.NE.0)GO TO 6 IF(II.LT.1)GO TO 3 DD=DST(II) DX=XXNEW IF(IBNEW.EQ.1)DX=DX-REPS IF(DD.LT.DX.AND.IK1.NE.0)GO TO 6 IF(DD.LT.DX)GO TO 3 ISUC=1 IF(METHOD.EQ.3)GO TO 49 GO TO 60 C C 6 IF(IK1.EQ.1)GO TO 7 IK1=1 P1=ANEW P2=AANEW IOLD1=INEW1 GO TO 59 7 IK1=0 XXOLD=XXNEW AAOLD=AANEW TTOLD=TTNEW IBOLD=IBNEW XXNEW=XNEW AANEW=ANEW TTNEW=TNEW IBNEW=0 GO TO 41 C C E N D O F I T E R A T I O N S C C C BOUNDARY RAYS. CASE IOLD.NE.ITPR, INEW=ITPR C 50 XXNEW=XNEW TTNEW=TNEW AANEW=ANEW IBNEW=0 P1=AOLD P2=ANEW 54 IRES=0 ITER=0 51 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC) IF(IND.EQ.20.OR.IND.EQ.21)RETURN if(itpr.eq.22)x=z IF(MOUT.EQ.2)WRITE(LOU,103)IND,IND1,ITER,X,T,PNEW IF(IND.EQ.ITPR)GO TO 52 P1=PNEW GO TO 53 52 XXOLD=X AAOLD=PNEW TTOLD=T IBOLD=1 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 XAX=X ZAX=Z IAX=IND IAX1=IND1 T1=T P2=PNEW 53 IF(ITER.LT.ITMAX)GO TO 51 IF(MOUT.EQ.1.AND.IRES.EQ.1) 1WRITE(LOU,107)XAX,ZAX,T1,P2,IAX,IAX1,IRES if(itpr.eq.22)z=x if(itpr.eq.22)x=xvsp IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,107)X,Z,T,PNEW,IND,IND1,IRES if(itpr.eq.22)x=z IF(IRES.EQ.1)GO TO 41 GO TO 3 C C BOUNDARY RAYS. CASE IOLD=ITPR, INEW.NE.ITPR C OR CASE IOLD=ITPR, INEW=ITPR BUT IOLD1.NE.INEW1 C 55 XXOLD=XOLD AAOLD=AOLD TTOLD=TOLD IBOLD=0 P1=AOLD P2=ANEW 59 IRES=0 ITER=0 56 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC) IF(IND.EQ.20.OR.IND.EQ.21)RETURN if(itpr.eq.22)x=z IF(MOUT.EQ.2)WRITE(LOU,103)IND,IND1,ITER,X,T,PNEW IF(IND.EQ.ITPR.AND.IBOLD.EQ.1)GO TO 57 IF(IND.EQ.ITPR.AND.IND1.EQ.IOLD1)GO TO 57 P2=PNEW GO TO 58 57 XXNEW=X AANEW=PNEW TTNEW=T IBNEW=1 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 XAX=X ZAX=Z IAX=IND IAX1=IND1 T2=T P1=PNEW 58 IF(ITER.LT.ITMAX)GO TO 56 IF(MOUT.EQ.1.AND.IRES.EQ.1) 1WRITE(LOU,107)XAX,ZAX,T2,P1,IAX,IAX1,IRES if(itpr.eq.22)z=x if(itpr.eq.22)x=xvsp IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,107)X,Z,T,PNEW,IND,IND1,IRES if(itpr.eq.22)x=z IF(IRES.EQ.1.AND.IK1.EQ.1)GO TO 7 IF(IRES.EQ.1)GO TO 41 GO TO 3 C C CRITICAL RAYS. CASE IOLD.NE.9, IOLD.NE.ITPR, INEW=9 C 30 ITER=0 XCR=XNEW P1=AOLD P2=ANEW IRES=0 31 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC) IF(IND.EQ.20.OR.IND.EQ.21)RETURN if(itpr.eq.22)x=z IF(MOUT.EQ.2)WRITE(LOU,104)IND,IND1,ITER,X,T,PNEW IF(IND.EQ.9)GO TO 32 IF(IND.EQ.ITPR)GO TO 33 P1=PNEW GO TO 34 32 IF(IND1.NE.INEW1)P1=PNEW IF(IND1.NE.INEW1)GO TO 34 P2=PNEW IF(ABS(X-XCR).LT..01.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 89 XCR=X GO TO 34 89 DO 86 K=1,KREF 86 JC(K)=CODE(K) IRF3=IREF+3 DO 87 K=IRF3,KREF 87 CODE(K-2)=JC(K) KREF1=KREF KREF=KREF-2 ICR=1 ITER=ITMAX-1 GO TO 31 33 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 XAX=X ZAX=Z IAX=IND IAX1=IND1 T2=T P1=PNEW PAP=PNEW 34 IF(ITER.LT.ITMAX)GO TO 31 IF(MOUT.EQ.1.AND.IRES.EQ.1) 1WRITE(LOU,111)XAX,ZAX,T2,PAP,IAX,IAX1,IRES if(itpr.eq.22)z=x if(itpr.eq.22)x=xvsp IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,111)X,Z,T,PNEW,IND,IND1,IRES if(itpr.eq.22)x=z IF(IRES.EQ.0)GO TO 3 P2=PAP XXNEW=XAX AANEW=P2 TTNEW=T2 IBNEW=1 P1=AOLD IF(ICR.EQ.0)GO TO 54 ICR=0 KREF=KREF1 DO 88 K=1,KREF 88 CODE(K)=JC(K) GO TO 54 C C CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.ITPR C 35 ITER=0 XCR=0 P1=AOLD P2=ANEW IRES=0 36 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC) IF(IND.EQ.20.OR.IND.EQ.21)RETURN if(itpr.eq.22)x=z IF(MOUT.EQ.2)WRITE(LOU,104)IND,IND1,ITER,X,T,PNEW IF(IND.EQ.9)GO TO 37 IF(IND.EQ.ITPR)GO TO 38 P2=PNEW GO TO 39 37 IF(IND1.NE.IOLD1)P2=PNEW IF(IND1.NE.IOLD1)GO TO 39 P1=PNEW IF(ABS(X-XCR).LT..01.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 94 XCR=X GO TO 39 94 DO 91 K=1,KREF 91 JC(K)=CODE(K) IRF3=IREF+3 DO 92 K=IRF3,KREF 92 CODE(K-2)=JC(K) KREF1=KREF KREF=KREF-2 ICR=1 ITER=ITMAX-1 GO TO 36 38 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 XAX=X ZAX=Z IAX=IND IAX1=IND1 P2=PNEW PAP=PNEW T1=T 39 IF(ITER.LT.ITMAX)GO TO 36 IF(MOUT.EQ.1.AND.IRES.EQ.1) 1WRITE(LOU,111)XAX,ZAX,T1,PAP,IAX,IAX1,IRES if(itpr.eq.22)z=x if(itpr.eq.22)x=xvsp IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,111)X,Z,T,PNEW,IND,IND1,IRES if(itpr.eq.22)x=z IF(IRES.EQ.0)GO TO 3 P1=PAP XXOLD=XAX AAOLD=P1 TTOLD=T1 IBOLD=1 P2=ANEW IF(ICR.EQ.0)GO TO 59 ICR=0 KREF=KREF1 DO 93 K=1,KREF 93 CODE(K)=JC(K) GO TO 59 C C 100 FORMAT(2I3,2F10.5,F15.10) 101 FORMAT(1X,'ITERATIONS',5X,3I3,3F10.5,F15.10) 102 FORMAT(/,1X,'SUCCESSFUL ITERATION',3X,5I3,3F10.5,F15.10) 103 FORMAT(2X,'BOUNDARY RAY',5X,3I3,2F10.5,F15.10) 104 FORMAT(2X,'CRITICAL RAY',5X,3I3,2F10.5,F15.10) 105 FORMAT(4F10.5,2F10.6,F10.3,5I3/13X,3(E15.5,F10.5)) 106 FORMAT(2I5,5E15.10) 107 FORMAT(10X,3F10.5,F15.10,3I3,5X,'BOUNDARY RAY') 108 format(2i5) 109 format(2e15.5) 110 FORMAT(5X,3(E15.5,F10.5),2F15.7,I5/) 111 FORMAT(10X,3F10.5,F15.10,3I3,5X,'CRITICAL RAY') 113 FORMAT(4F10.5,F15.10,4I3) 114 FORMAT(/) 115 FORMAT(2I3,F10.5,2E12.6,4F10.6) 116 FORMAT(2I3,F10.5,3E12.6,5F10.6) 117 format(i5,10e15.5) C C 99 N=0 NAUX=0 if(lu1.eq.0)return if(ltape.eq.0)WRITE(LU1,108)N,NAUX if(ltape.eq.1)WRITE(LU1)N,NAUX IF(CODE(1).LT.0)INDEX=-INDEX IF(ltape.eq.0)WRITE(LU1,108)INDEX IF(ltape.eq.1)WRITE(LU1)INDEX IF(INDEX.EQ.0)RETURN IF(INDEX.LT.0)INDEX=-INDEX IF(ltape.eq.0)WRITE(LU1,117)(INDI(I),XCOOR(I),TIME(I),TAS(I), 1ANG(I),AMPX(I),AMPY(I),AMPZ(I),PHAX(I),PHAY(I),PHAZ(I),I=1,INDEX) IF(ltape.eq.1)WRITE(LU1)(INDI(I),XCOOR(I),TIME(I),TAS(I),ANG(I), 1AMPX(I),AMPY(I),AMPZ(I),PHAX(I),PHAY(I),PHAZ(I),I=1,INDEX) RETURN END C C ************************************************************** C SUBROUTINE GENER(N,IDP,IDS,IBP,IBS,IUP,IUS,IREAD,IS,IR,NS,NL,IFIN, 1NCODE,lin,lou) C C THE ROUTINE GENER IS DESIGNED TO GENERATE NUMERICAL CODES OF C ELEMENTARY SEISMIC BODY WAVES C DIMENSION JCA(50) COMMON/COD/JC(50),KCA,KC C C IFIN=0 if(is.le.0.or.is.gt.ns)go to 101 if(ir.lt.0.or.ir.gt.ns)go to 101 irc=ir IF(N.GT.0)GO TO 1 NCODE=1 INEW=1 INDEX=1 GO TO 51 1 NCODE=NCODE+1 IF(INDEX.GE.11)GO TO 80 IF(INEW.EQ.1)GO TO 50 JAUX=0 DO 20 I=1,KCA IF(JCA(I).EQ.JAUX)GO TO 21 IAUX=I 20 JAUX=JCA(I) 21 IF(INDEX.LE.6.AND.JAUX.GE.(NS-1))GO TO 50 IF(INDEX.EQ.4.AND.NL.EQ.0.AND.JAUX.GE.(NS-2))GO TO 50 IF(INDEX.EQ.6.AND.NL.EQ.0.AND.JAUX.GE.(NS-2))GO TO 50 IF(INDEX.GT.6.AND.JAUX.EQ.1)GO TO 50 IAUX1=IAUX+1 IAUX2=IAUX+2 DO 22 I=IAUX1,KCA II=KCA+IAUX1-I JC(II+2)=JC(II) 22 JCA(II+2)=JCA(II) KCA=KCA+2 IF(INDEX.GE.7)GO TO 31 JC(IAUX1)=JAUX+1 JC(IAUX2)=JAUX+1 JCA(IAUX1)=JC(IAUX1) JCA(IAUX2)=JC(IAUX2) IF(INDEX.GE.5)JC(IAUX1)=-JC(IAUX1) IF(INDEX.EQ.4.OR.INDEX.EQ.5)JC(IAUX2)=-JC(IAUX2) RETURN 31 JC(IAUX1)=JAUX-1 JC(IAUX2)=JAUX-1 JCA(IAUX1)=JC(IAUX1) JCA(IAUX2)=JC(IAUX2) IF(INDEX.GE.9)JC(IAUX1)=-JC(IAUX1) IF(INDEX.EQ.8.OR.INDEX.EQ.9)JC(IAUX2)=-JC(IAUX2) RETURN C 50 INDEX=INDEX+1 INEW=1 IF(INDEX.GE.11)GO TO 80 51 IF(INDEX.EQ.1.AND.IDP.EQ.0)GO TO 50 IF(INDEX.EQ.2.AND.IDS.EQ.0)GO TO 50 IF(INDEX.EQ.3.AND.IBP.EQ.0)GO TO 50 IF(INDEX.EQ.4.AND.IBP.NE.2)GO TO 50 IF(INDEX.EQ.5.AND.IBS.EQ.0)GO TO 50 IF(INDEX.EQ.6.AND.IBS.NE.2)GO TO 50 IF(INDEX.EQ.7.AND.IUP.EQ.0)GO TO 50 IF(INDEX.EQ.8.AND.IUP.NE.2)GO TO 50 IF(INDEX.EQ.9.AND.IUS.EQ.0)GO TO 50 IF(INDEX.EQ.10.AND.IUS.NE.2)GO TO 50 if(index.ge.7)irc=ir-1 IMIN=IS IMAX=IRc IF(IS.GT.IRc)IMIN=IRc IF(IS.GT.IRc)IMAX=IS ISUM=IMAX-IMIN+2 IF(INDEX.LE.2)GO TO 70 INEW=0 IF(INDEX.GT.6)GO TO 60 ISM=IMAX-IS+1 DO 55 I=1,ISM JCA(I)=IS+I-1 JC(I)=JCA(I) 55 IF(INDEX.GT.4)JC(I)=-JCA(I) ISM2=ISM+1 DO 56 I=ISM2,ISUM JCA(I)=IMAX+ISM2-I JC(I)=JCA(I) 56 IF(INDEX.EQ.4.OR.INDEX.EQ.5)JC(I)=-JCA(I) KCA=ISUM KC=1 RETURN 60 ISM=IS-IMIN+1 DO 65 I=1,ISM JCA(I)=IS+1-I JC(I)=JCA(I) 65 IF(INDEX.GT.8)JC(I)=-JCA(I) ISM2=ISM+1 DO 66 I=ISM2,ISUM JCA(I)=IMIN+I-ISM2 JC(I)=JCA(I) 66 IF(INDEX.EQ.8.OR.INDEX.EQ.9)JC(I)=-JCA(I) KCA=ISUM KC=-1 INEW=0 RETURN 70 KCA=IMAX-IMIN+1 KC=1 IF(IS.GE.IR)KC=-1 DO 71 I=1,KCA JC(I)=IS+I-1 IF(IS.GT.IR)JC(I)=IS-I+1 JCA(I)=JC(I) IF(INDEX.EQ.2)JC(I)=-JC(I) 71 CONTINUE RETURN 80 IF(IREAD.EQ.0)GO TO 101 KCA=0 READ(LIN,*)KC,KCA,(JC(I),I=1,KCA) IF(KCA.EQ.0)GOTO 101 WRITE(LOU,100)KC,KCA,(JC(I),I=1,KCA) IF(KCA.GT.50)GO TO 99 DO 81 I=1,KCA 81 JCA(I)=IABS(JC(I)) JAUX=JCA(1) IF(KC.EQ.0)RETURN IF(ir.gt.0.and.(JCA(KCA).GT.IR.OR.JCA(KCA).LT.(IR-1)))GO TO 99 IF(KCA.EQ.1)RETURN ID=1 IF(KC.LT.0)ID=-1 DO 103 I=2,KCA JJ=JCA(I) IF(JJ.LE.0.OR.JJ.GE.NS)GO TO 99 IF(JJ.EQ.JAUX)GO TO 105 IF(ID.EQ.1.AND.JJ.EQ.(JAUX+1))GO TO 103 IF(ID.EQ.(-1).AND.JJ.EQ.(JAUX-1))GO TO 103 GO TO 99 105 ID=-ID 103 JAUX=JJ RETURN 99 IFIN=-1 RETURN 101 IFIN=1 RETURN 100 FORMAT(24I3) END C C======================================================================= C INCLUDE 'ray.for' C ray.for C C Interpolation method: C Include just one of the following files 'mod*.for': C (a) Isosurface interpolation: INCLUDE 'modis.for' C modis.for C (b) Bicubic B-spline interpolation: * INCLUDE 'modsp.for' C modsp.for C C======================================================================= C