C
C PROGRAM S Y N T A N C ********************** C C PROGRAM SYNTAN IS DESIGNED FOR THE COMPUTATION OF RAY C SYNTHETIC SEISMOGRAMS FROM THE IMPULSE SYNTHETIC SEISMOGRAMS C STORED IN THE FILE LU2 GENERATED BY THE PROGRAM ANRAY C C ************************************************************** C CHARACTER*80 MTEXT,FILEIN,FILEOU,FILE1,FILE2,FILE3 COMPLEX CAUX,CAX,CAY,CAZ,CAX1,CAY1,CAZ1,AUX,CSOUR,IMAG DIMENSION IDIST(2000),TIME(2000),A(2000),PH(2000), 1SEIS(3001),TAS(2000),DIST(100),SMAX(100),ISEL(100),XX(100), 2IEP(100),JS(20),AUX(3),p(3),pol(3),pol1(3) COMMON/SOUR/ROS C C************************************************** C LIN=5 LOU=6 LU2=1 LU4=2 LU5=3 FILEIN='syntan.dat' FILEOU='syntan.out' FILE1='lu2.dat' FILE2='lu4.dat' FILE3='lu5.dat' WRITE(*,'(2A)') ' (SYNTAN) SPECIFY NAMES OF INPUT AND OUTPUT', 1' FILES LIN, LOU, LU2, LU4, LU5: ' READ(*,*) FILEIN,FILEOU,FILE1,FILE2,FILE3 IF(FILE1.EQ.' ') GO TO 99 IF(FILE2.EQ.' ') LU4=0 IF(FILE3.EQ.' ') LU5=0 OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') OPEN(LU2,FILE=FILE1,FORM='FORMATTED',STATUS='OLD') IF(LU4.NE.0)OPEN(LU4,FILE=FILE2,FORM='FORMATTED') IF(LU5.NE.0)OPEN(LU5,FILE=FILE3,FORM='UNFORMATTED') C C************************************************** C IMAG=(0.,1.) PI=3.14159265 SHF=0. 1 READ(LIN,*)MPRINT WRITE(LOU,100)MPRINT C C *** C REWIND LU2 IF(LU4.NE.0)REWIND LU5 IF(LU4.NE.0)REWIND LU4 2 READ(LU2,101)MTEXT READ(LU2,100)NDST,KSH,ILOC READ(LU2,102)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS IF(NDST.EQ.1)THEN READ(LU2,102)XREC,YREC,XPRF,XATAN ELSE READ(LU2,102)(DIST(K),K=1,NDST) END IF C WRITE(LOU,101)MTEXT WRITE(LOU,100)NDST,KSH,ILOC WRITE(LOU,102)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS IF(NDST.EQ.1)WRITE(LOU,102)XREC,YREC,XPRF,XATAN IF(NDST.NE.1)WRITE(LOU,102)(DIST(K),K=1,NDST) IF(NDST.EQ.1)DIST(1)=XPRF READ(LIN,*)MCOMP,MRED,MSELEC,MEPIC,MSHIFT,KABS,MSOUR WRITE(LOU,100)MCOMP,MRED,MSELEC,MEPIC,MSHIFT,KABS,MSOUR READ(LIN,*)TMIN,DT,TMAX,FREQ,GAMA,PSI,VRED,SHIFT,AROT WRITE(LOU,102)TMIN,DT,TMAX,FREQ,GAMA,PSI,VRED,SHIFT,AROT CSRT=COS(AROT) SNRT=SIN(AROT) IF(MEPIC.EQ.0)GO TO 11 READ(LIN,*)NEPIC,(IEP(I),I=1,NEPIC) WRITE(LOU,100)NEPIC,(IEP(I),I=1,NEPIC) 11 IF(MSELEC.EQ.0)GO TO 12 READ(LIN,*)NSELEC,(ISEL(I),I=1,NSELEC) WRITE(LOU,100)NSELEC,(ISEL(I),I=1,NSELEC) 12 CONTINUE IF(KABS.EQ.0)GO TO 5 READ(LIN,*) FREF,QRED WRITE(LOU,102)FREF,QRED IF(FREF.LT..00001)FREF=1. IF(QRED.LT..00001)QRED=1. RF=FREQ/FREF C 5 IF(MSOUR.NE.0)CALL SOURCE(LIN,LOU,0,0,MSOUR,P,POL,AMSOUR,PHSOUR) C NCD=0 I=0 3 I=I+1 IF(I.GT.2000)WRITE(LOU,104) IF(I.GT.2000)GO TO 99 READ(LU2,116,END=7)NC,IDIST(I),TIME(I),CAX,CAY,CAZ,TAST IF(NC.LT.0)READ(LU2,117)CAX1,CAY1,CAZ1 READ(LU2,118)(P(K),K=1,3) READ(LU2,118)(POL(K),K=1,3) IF(NC.LT.0)READ(LU2,118)(POL1(K),K=1,3) CAUX=CAX*CSRT+CAY*SNRT CAY=-CAX*SNRT+CAY*CSRT CAX=CAUX CAUX=CAX1*CSRT+CAY1*SNRT CAY1=-CAX1*SNRT+CAY1*CSRT CAX1=CAUX NC1=IABS(NC) NC1=NC1+NCD IF(MSELEC.EQ.0)GO TO 14 DO 13 J=1,NSELEC IF(ISEL(J).EQ.NC1)GO TO 14 13 CONTINUE GO TO 3 14 IF(MEPIC.EQ.0)GO TO 16 DO 15 J=1,NEPIC IF(IEP(J).EQ.IDIST(I))GO TO 16 15 CONTINUE GO TO 3 16 CONTINUE C IF(MSOUR.NE.0)THEN CALL SOURCE(LIN,LOU,1,1,MSOUR,P,POL,AMDIR,PHDIR) CSOUR=AMDIR*CEXP(IMAG*PHDIR) END IF IF(MSOUR.EQ.0)CSOUR=(1.,0.) AUX(1)=CAX*CSOUR AUX(2)=CAY*CSOUR AUX(3)=CAZ*CSOUR IF(NC.LT.0)THEN IF(MSOUR.NE.0)THEN CALL SOURCE(LIN,LOU,1,1,MSOUR,P,POL1,AMDIR,PHDIR) CSOUR=AMDIR*CEXP(IMAG*PHDIR) END IF IF(MSOUR.EQ.0)CSOUR=(1.,0.) AUX(1)=AUX(1)+CAX1*CSOUR AUX(2)=AUX(2)+CAY1*CSOUR AUX(3)=AUX(3)+CAZ1*CSOUR END IF IF(MCOMP.EQ.0)K=3 IF(MCOMP.EQ.1)K=1 IF(MCOMP.EQ.2)K=2 RW=REAL(AUX(K)) CW=AIMAG(AUX(K)) A(I)=SQRT(RW*RW+CW*CW) IF(ABS(RW).LT..000001.AND.ABS(CW).LT..000001)THEN PH(I)=0. ELSE PH(I)=ATAN2(CW,RW) END IF IF(MPRINT.GE.1)WRITE(LOU,114) 1NC,IDIST(I),TIME(I),A(I),PH(I),P,TAST IF(KABS.NE.0)TAS(I)=TAST IF(KABS.EQ.2)PH(I)=PH(I)-2.*FREQ*TAST*ALOG(RF) GO TO 3 C 7 ISUM=I-1 IF(VRED.LT.0.0001)VRED=8. IF(FREQ.LT.0.0001)FREQ=4. IF(GAMA.LT.0.0001)GAMA=4.0 IF(MSHIFT.EQ.0)SHF=0. IF(MSHIFT.EQ.1)SHF=.241506*GAMA/FREQ IF(MSHIFT.EQ.2)SHF=SHIFT OMEGA=2.*PI*FREQ OMRF=2.*PI*FREF OG=OMEGA/GAMA OGG=OG*OG TSH=0.45*GAMA/FREQ C C CONSTRUCTION OF SYNTHETIC SEISMOGRAMS C IF(NDST.GT.100)WRITE(LOU,105) IF(NDST.GT.100)GO TO 99 NPTS=(TMAX-TMIN)/DT+1.5 IF(NPTS.GT.3001)WRITE(LOU,106) IF(NPTS.GT.3001)GO TO 99 C C MMD=NDST IF(MEPIC.NE.0)MMD=NEPIC IF(LU4.EQ.0)GO TO 6 WRITE(LU5)NPTS,TMIN,DT,NDST,DIST(1),RSTEP WRITE(LU4,101)MTEXT WRITE(LU4,108)MMD,MRED,MCOMP,ILOC,VRED,RSTEP,XSOUR,YSOUR,DT C C LOOP FOR RECEIVER POSITIONS C 6 CONTINUE DO 30 I=1,MMD JJ=I IF(MEPIC.NE.0)JJ=IEP(I) XX(I)=DIST(JJ) XXD=ABS(XX(I)-XSOUR) XXV=0. IF(MRED.NE.0)XXV=XXD/VRED AMAX=0. DO 20 J1=1,NPTS 20 SEIS(J1)=0. DO 21 K=1,ISUM IF(IDIST(K).NE.JJ)GO TO 21 TR=TIME(K)-XXV+SHF IF((TR+TSH).LT.TMIN)GO TO 21 IF((TR-TSH).GT.TMAX)GO TO 21 AMPL=A(K) PHASE=PH(K) IF(AMPL.LT.0.00005*AMAX)GO TO 21 IF(AMPL.GT.AMAX)AMAX=AMPL IF1=IFIX((TR-TSH-TMIN)/DT+0.1) IF2=IFIX((TR+TSH-TMIN)/DT+0.1) IFF1=1 IF(IF1.GT.0)IFF1=IF1 IF1=IFF1 IFF2=NPTS IF(IF2.LT.NPTS)IFF2=IF2 IF2=IFF2 IF(KABS.EQ.0)GO TO 27 TAST=TAS(K)*QRED OMAS=OMEGA*(1.-OMEGA*TAST/(GAMA*GAMA)) IF(OMAS.LE.0.)WRITE(LOU,112) IF(OMAS.LE.0.)GO TO 32 AABS1=OMEGA*TAST/GAMA AABS1=.25*(AABS1*AABS1) AABS2=.5*OMAS*TAST IF(KABS.NE.2)GO TO 27 AABS4=OMAS/OMRF AABS5=ALOG(AABS4) AABS3=(TAST/PI)*(1.+AABS5) 27 DO 24 KK=IF1,IF2 TT=TMIN+DT*FLOAT(KK-1) TD=TT-TR IF(KABS.EQ.2)TD=TD+AABS3 AA=PSI-PHASE BB=-OGG*TD*TD IF(KABS.EQ.0)AA=AA+OMEGA*TD IF(KABS.NE.0)AA=AA+OMAS*TD IF(KABS.EQ.2)AA=AA-OMAS*TAST/PI IF(KABS.NE.0)BB=BB-AABS1-AABS2 AA=AMPL*COS(AA)*EXP(BB) 24 SEIS(KK)=SEIS(KK)+AA 21 CONTINUE SMAX(I)=0. DO 25 KK=1,NPTS IF(SMAX(I).GT.ABS(SEIS(KK)))GO TO 25 SMAX(I)=ABS(SEIS(KK)) 25 CONTINUE C C IF(LU4.NE.0)WRITE(LU5)(SEIS(J),J=1,NPTS) 30 CONTINUE REWIND LU2 IF(LU4.NE.0)REWIND LU5 C C END OF THE LOOP FOR RECEIVERS C SMAXIM=0. DO 26 I=1,MMD IF(SMAXIM.GT.SMAX(I))GO TO 26 SMAXIM=SMAX(I) XXX=XX(I) 26 CONTINUE IF(LU4.NE.0)WRITE(LU4,113)XXX,SMAXIM WRITE(LOU,113)XXX,SMAXIM IF(LU4.EQ.0)GO TO 1 C IF(LU4.NE.0)READ(LU5)NPTS,TMIN,DT,NDST,DIST(1),RSTEP DO 31 I=1,MMD IF(LU4.NE.0)READ(LU5)(SEIS(J),J=1,NPTS) C C NORMALIZATION OF SEISMOGRAMS C IF1=0 IF2=0 SMAXI=SMAX(I) IF(SMAXI.LT..000001)GO TO 35 DO 34 J=1,NPTS SEIS(J)=999.1*SEIS(J)/SMAXI IF(ABS(SEIS(J)).LT.1.)GO TO 33 IF2=0 IF(IF1.EQ.0.AND.J.EQ.1)IF1=1 IF(IF1.EQ.0)IF1=J-1 GO TO 34 33 IF(IF1.EQ.0)GO TO 34 IF(IF2.EQ.0)IF2=J 34 CONTINUE TM=TMIN+FLOAT(IF1-1)*DT IF(IF2.EQ.0)IF2=NPTS NPS=IF2-IF1+1 35 IF(IF1.EQ.0)NPS=0 IF(IF1.EQ.0)TM=TMIN C WRITE(LU4,109)XX(I),SMAX(I),TM,NPS WRITE(LOU,107)XX(I),SMAX(I),TM,NPS ISS=20 IS=IF1-1 41 IS1=IF2-IS IF(IS1.LT.20)ISS=IS1 IF(IS.LT.0)IS=0 DO 40 K=1,ISS 40 JS(K)=SEIS(IS+K) WRITE(LU4,111)(JS(K),K=1,ISS) IF(MPRINT.EQ.2.AND.NPS.NE.0) 1WRITE(LOU,110)(JS(K),K=1,ISS) IF(IS1.LE.20)GO TO 31 IS=IS+20 GO TO 41 31 CONTINUE REWIND LU4 32 CONTINUE C C 100 FORMAT(26I3) 101 FORMAT(A) 102 FORMAT(8F10.5) 104 FORMAT(1X,'NUMBER OF READINGS FROM LU2 GREATER THAN 2000.'/, 11X,'COMPUTATION TERMINATED'/) 105 FORMAT(1X,'NUMBER OF RECEIVER POSITIONS GREATER THAN 100.',/, 11X,'COMPUTATION TERMINATED'/) 106 FORMAT(1X,'NUMBER OF POINTS IN THE SYNTHETIC SEISMOGRAM 1GREATER THAN 3001.'/1X,'COMPUTATION TERMINATED'/) 107 FORMAT(1X,'EPICENTRAL DISTANCE =',F10.5,E15.9,1X,F10.5,I5) 108 FORMAT(4I5,5F10.5) 109 FORMAT(F10.5,E15.8,F10.5,I5) 110 FORMAT(1X,20I4) 111 FORMAT(20I4) 112 FORMAT(1X,'FREQUENCY TOO HIGH OR GAMA TOO SMALL'/ 11X,'COMPUTATION TERMINATED'/) 113 FORMAT(1X,'EPICENTRAL DISTANCE =',F10.5,1X,'SMAXIM =',E15.9) 114 FORMAT(2I3,F10.5,E12.6,5F10.6) 116 FORMAT(2I3,F12.7,6E12.6,F10.6) 117 FORMAT(6E12.6) 118 FORMAT(3F10.5) C 99 CONTINUE STOP END C C======================================================================= C INCLUDE 'source.for' C source.for C C======================================================================= C