C C PROGRAM TO CALCULATE PHASE VELOCITIES FOR A GIVEN SOURCE MECHANISM C AND A MODEL FOR WHICH EIGENFUNCTIONS ARE AVAILABLE FROM MINOS C SINGLE MODE. B.ROMANOWICZ OCTOBER 1986. C c IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*36 IKRI CHARACTER*36 IPH,ICONE CHARACTER*36 IDAT,IRIS,IVOIRE,IGLOO CHARACTER*36 INER CHARACTER*80 IPR,IIN,INREF INTEGER*2 DATFIL(18),IST(50,3) INTEGER*2 ITITLE(30),OUTFIL(18),IREST(18),IWX(18),IXI(18),ISTO(18) INTEGER*2 ILLV(18),IVG(18),JIAK(2),IVICT(4) DIMENSION T1(4),TE1(4),TE2(4),TW1(4),TW2(4),ERO(500),CTC(500) COMMON/OFT/CZ(500),CW(500),TJ(500),TV(500) COMMON/TDKP/TOUL(20),QI1(20),QI2(20),JOY,RAP(20) COMMON/COP/GLATE,GL1,GLATS,GL2 COMMON/PULSE/TJP(4,500),DWP(4,500),EP(4,500),UP(4,500) COMMON/DOCH/ DP(500),DON(1000),CU(500),GN1(5) COMMON/OUF/D(500),XM(6),TC(500),CC(500),P(500),CV(500,20) COMMON/MODD/TB(500),R(500),NRAD,NMOD,MOD(20),NTOT,EIP(500) *,NF(4),NLA(4) COMMON/BITS/PI,RN,VN,WN,GN,FL,FL1,FL2,FL3,SFL3,JCOM COMMON/DATA/ASS(5,500),FAIS(5,500),DIST(50),AZ(50) *,AM(50),DT(50) COMMON/EIGEN/AS(500),AP(500),AQ(500),CT(500),UU(500),Q(500) *,TT(500) COMMON/SPECTR/RER(5,500),AIMR(5,500),X(5,500),XX(5,500),XM0 COMMON/COOL/EIGL(500),IWW COMMON/VIEL/DWGC(4,500),VPGC(4,500),ER(4),TIMP,IWQ EQUIVALENCE(OUTFIL,IPH) EQUIVALENCE(IWX,IRIS) EQUIVALENCE(IXI,IVOIRE) EQUIVALENCE(ISTO,IGLOO) EQUIVALENCE(ILLV,ICONE) EQUIVALENCE(IVG,INER) EQUIVALENCE(DATFIL,IDAT) EQUIVALENCE(IREST,IKRI) DATA BIGG/6.6723D-11/,RHOBAR/5515.D0/,NMIN/5/ PI=3.14159265358979D0 RAD=180./PI TPI=2.*PI PRINT*,'CALCULATE PHASE VEL.(YES=1)' READ(5,*)IPHAS PRINT*,'ANTIPODAL DATA(YES=1)' READ(5,*)IPOD IF(IPOD.EQ.1)WRITE(1,115) C IF(IPHAS.EQ.0)GO TO 8 PRINT*,'PHASEVEL FILENAME' READ(5,100)OUTFIL OPEN(12,FILE=IPH) 8 CONTINUE 115 FORMAT(' ANTIPODAL DATA') print*,'MODEL FILE NAME' READ(5,'(a80)') IIN 100 FORMAT(18A2) OPEN(7,FILE=IIN) print*,'MODEL TITLE' READ(5,1000)ITITLE 1000 FORMAT(30A2) print*,'title',ititle print*,'DATA FILE NAME' READ(5,100) DATFIL OPEN(15,FILE=IDAT) print*,'ENTER JCOM :' READ(5,*) JCOM print*,' ENTER OPTION (IOPT=2 FOR MODES)' READ(5,*)IOPT print*,'iopt',iopt c PRINT*,'CORRECTION ELLIPTICITE O=1 OU N=0' c READ*,IWW PRINT*,'LAT. ET LONG. DE L"EPICENTRE ' READ*,NOE,GLATE,GL1 PRINT*,'INT. PERIODES A CONSERVER ' READ*,(TE1(I),TE2(I),I=1,4) READ*,(ER(I),I=1,4) 50 FORMAT(I5/(6E13.5)) C READ MODEL RADII READ(7,200)(ITITLE(I),I=1,20) WRITE(6,200)(ITITLE(I),I=1,20) 200 FORMAT(20A2) READ(7,112)IFANIS,IREF,IFDECK 112 FORMAT(I1,I4,I2) WRITE(6,111)IFANIS,IREF,IFDECK 111 FORMAT(' IFANIS',3I5) READ(7,113)NRAD,NIC,NOC 113 FORMAT(3I5) WRITE(6,111)NRAD,NIC,NOC READ(7,105)(R(I),I=1,NRAD) 105 FORMAT(F8.0) C WRITE(6,60)NRAD,(R(I),I=1,NRAD) 60 FORMAT(I5/(6E12.3)) C NORMALIZE RN=R(NRAD) GN=PI*BIGG*RHOBAR*RN VN2=GN*RN VN=SQRT(VN2) WN=VN/RN C C LECTURE DES VITESSES DE PHASE DE REFERENCE C PRINT*,'REFERENCE PHASEVEL FILENAME' READ(5,'(a80)')INREF OPEN(11,FILE=INREF) READ(11,120)IPER 120 FORMAT(I5) WRITE(6,120)IPER READ(11,56)IMOD WRITE(6,56)IMOD 56 FORMAT(36A2) READ(11,120)NC WRITE(6,120)NC READ(11,*)(TC(I),CC(I),I=1,NC) C WRITE(6,28)(TC(I),CC(I),I=1,NC) 28 FORMAT(F7.2,F8.4) IF(IPOD.EQ.1)GO TO 98 print*,'STRIKE DIP SLIP DEP TAU XM' READ(5,*)STRK,DIP,SLIP,DEP,TAU,XM0 WRITE(6,555)STRK,DIP,SLIP,dep,tau,xm 555 FORMAT(' STRK DIP SLIP DEP TAU MO'/(6F9.3)) C CALL MOMENT(DIP,SLIP,STRK,XM,XMT) CALL MOMENT(STRK,DIP,SLIP,XM) 98 CONTINUE C C LECTURE DES DONNEES SUR L'UNITE 5 C READ(15,1000)ITITLE JIAK(1)=ITITLE(1) JIAK(2)=ITITLE(2) WRITE(6,55)ITITLE WRITE(12,55)ITITLE WRITE(16,55)ITITLE 55 FORMAT(36A2) READ(15,300)NTOT C WRITE(1,300)NTOT 300 FORMAT(I5) READ(15,102)(TB(I),I=1,NTOT) c WRITE(6,102)(TB(I),I=1,NTOT) c write(6,102)(tc(i),i=1,nc) DO 99 L=1,NTOT PER=TB(L) CT(L)=YNTERP(TC,CC,PER,NC,0) c print*,'tb ct',tb(l),ct(l) 99 CONTINUE 102 FORMAT(6(1X,E12.5)) I=1 80 READ(15,*,END=900)NOPT,GN1(I),GN2 READ(15,102)(DON(KTA),KTA=1,NOPT) READ(15,12)IST(I,1),IST(I,2),IST(I,3),IY0,ID0,TS,DT(I), *NTOT,DIST(I),AZ(I),AM(I),STLAT,STLON 12 FORMAT(A2,2A1,1X,2I4,F10.3,F5.1,I6,2F10.3,F5.0,2F8.2) WRITE(6,12)IST(I,1),IST(I,2),IST(I,3),IY0,ID0,TS,DT(I),NTOT, *DIST(I),AZ(I),AM(I) READ(15,102)(ASS(I,K),FAIS(I,K),K=1,NTOT) C WRITE(1,102)(ASS(I,K),FAIS(I,K),K=1,NTOT) C T1(I)=DIST(I)*(GN1(I)-GN2)/GN2/GN1(I) C PRINT*,'T1',I,T1(I) C IF(T1(I).GT.TE2)T1(I)=TE2 I=I+1 GO TO 80 900 CONTINUE WRITE(12,5500)NOE,JIAK,IST(1,1),IST(1,2),STRK,DIP, *SLIP,DEP,TAU,XM0 WRITE(16,5500)NOE,JIAK,IST(1,1),IST(1,2),STRK,DIP, *SLIP,DEP,TAU,XM0 5500 FORMAT(I3,1X,2A2,1X,A2,A1,6(1X,F7.2)) WRITE(12,5510)GLATE,GL1,STLAT,STLON WRITE(16,5510)GLATE,GL1,STLAT,STLON 5510 FORMAT(4F8.2) GLATS=STLAT GL2=STLON NSTAT=I-1 WRITE(12,*)NSTAT WRITE(16,*)NSTAT 109 FORMAT(A2,2A1,2F8.3) 556 FORMAT('STRK',F9.3,' DIP',F9.3,' SLIP',F9.3) 557 FORMAT('DEPTH',F9.3,' S.TIME',F9.3,' S.MOMENT',F9.3) DO 84 I=1,NSTAT J=0 IH=INT(AM(I))+1 IDIOT=AM(I) IVICT(IH)=I TW1(I)=TE1(IH) TW2(I)=TE2(IH) IF(IDIOT.EQ.0)THEN DRIL=DIST(I) 561 FORMAT('DISTANCE',F9.3,' AZIMUT',F9.3) ENDIF 84 CONTINUE DO 247 I=1,4 TE1(I)=TW1(I) 247 TE2(I)=TW2(I) IF(IPOD.EQ.1)GO TO 63 C C ENTER EIGENFUNCTIONS C 1 CONTINUE CALL MODIN(JCOM,IOPT,DEP,DRIL,0.) 69 FORMAT(' TK SR PR QR',F9.3,3E12.5) 800 CONTINUE PRINT*,'ENTERING RADPAT' CALL RADPAT(XM,NSTAT,STRK) PRINT*,'EXIT RADPAT' DO 64 J=1,NTOT DO 64 I=1,NSTAT WRITE(1,65)I,J,RER(I,J),AIMR(I,J),TAU 65 FORMAT(' I J RER AIMR TAU',2I5,3E12.5) X(I,J)=RER(I,J)*RER(I,J)+AIMR(I,J)*AIMR(I,J) 1515 FORMAT(' I J Q DIST X ASS',2I5,E12.5,F10.4,2E12.5) X(I,J)=SQRT(X(I,J)) X(I,J)=X(I,J)*EXP(-Q(J)*DIST(I)) XX(I,J)=ATAN2(AIMR(I,J),RER(I,J)) COM=XX(I,J) XX(I,J)=XX(I,J)-(TAU*TPI/TB(J)) c print*,' dist x',tb(j),dist(i),x(i,j),xx(i,j) c modif 06/29/92: 666 continue if(xx(i,j).ge.0.)go to 677 xx(i,j)=xx(i,j)+tpi go to 666 677 continue c print*,' dist x',tb(j),dist(i),x(i,j),xx(i,j) C DFR=XX(I,J)*RAD c PRINT*,'X xX ',X(I,J),xx(i,J) 64 CONTINUE PRINT*,'OUT OF 64' PRINT*,'OUT OF 63' 63 CONTINUE C DO 66 J=1,NTOT C EIGL(J)=PI/TB(J)/UU(J) 66 CONTINUE 67 FORMAT(' XX',6E12.5) PRINT*,'ENTER T AND C(T) FOR 2 CONSEC.PERIODS' c READ(5,*)TT1,CT1,TT2,CT2 c modif hillegass c read(5,*)tt1 do l=1,ntot diff=abs(tt1-tb(l)) c print*,'tt1,tb',l,tt1,tb(l) if(diff.lt.1.e-4)go to 47 enddo 47 l1=l+1 tt2=tb(l1) ct1=ct(l) ct2=ct(l1) print*,'tt1 tt2 ct1 ct2',tt1,ct1,tt2,ct2,l,l1 IF(IPHAS.NE.1)THEN IF(IPHAS.EQ.0)GO TO 90 ENDIF c IF(IWW.EQ.1)THEN c CALL COPOLE(GTHE,PHIM,THER) c POL2=1.-3.*COS(GTHE)*COS(GTHE) C POL1=COS(GTHE) c READ(18,*)JOY c DO 48 IHI=1,JOY c48 READ(18,*)TOUL(IHI),QI1(IHI),QI2(IHI),RAP(IHI) c ENDIF C PRINT*,'THER',THER,'GTHE',GTHE,'PHIM',PHIM DO 70 I=1,NSTAT DIS=DIST(I) IDS=INT(AM(I)+1.) C JIQ=MOD(IDS,2) JIQ=IDS-IDS/2*2 IDS=IDS/2 IF(JIQ.EQ.0)JIQ=-1 C PRINT*,'THER',THER IF(IWW.EQ.1)DIS=6371.*(JIQ*THER+TPI*IDS) PRINT*,'DIS',DIS IF(IPOD.EQ.1)DIST(I)=6371.*PI*AM(I) FIREF=(TT2*CT2-TT1*CT1)*DIST(I)/TT1/TT2/CT1/CT2 NREF=-FIREF WRITE(6,779)NREF,FIREF 779 FORMAT(' NREF=',I5,'FIREF',F10.3) C FAC=DIST(I)/TT1 FAC=DIS/TT1 Q0=FAC/CT1*TPI C Q1=DIST(I)/TT2/CT2*TPI Q1=DIS/TT2/CT2*TPI QQ=ABS(Q1-Q0) WRITE(6,789)Q0,Q1,QQ EPS=1. IF(FAIS(I,1).LT.0.)EPS=-1. IP=1 IP0=0 P(1)=FAIS(I,1) ITEST=0 IF((AM(I)-1.).LT.1.E-4)GO TO 101 DO 72 J=2,NTOT EPS1=1. DIF=(FAIS(I,J)-FAIS(I,J-1))*EPS IF(DIF.GT.0.)EPS1=0. P(J)=FAIS(I,J)+EPS*(NREF*TPI*FLOAT(IP)+EPS1*TPI+IP0*TPI) PP=ABS(P(J)-P(J-1)) PRINT*,'PP QQ,EPS1',PP,QQ,EPS1 IF(ABS(EPS1).LT.0.1)GO TO 53 IF(PP.GT.QQ*1.2)THEN P(J)=P(J)-TPI*EPS EPS1=0. IF(ABS(FAIS(I,J)).LT.0.1)IP0=IP0-1 PRINT*,'IP0',IP0 ENDIF GO TO 54 53 CONTINUE IF(PP.LT.QQ*0.8)THEN EPS1=1. P(J)=P(J)+TPI*EPS ENDIF IF(PP.GT.QQ*1.3)THEN P(J)=P(J)-TPI*EPS IP0=IP0-1 ENDIF ITEST=1 C PRINT*,'ITEST,EPS1,IP0',ITEST,EPS1,IP0 54 CONTINUE 544 FORMAT(' TB P FAIS IP EPS1 IP0',3E12.5,I5,F7.2,I5) IP=IP+1 IF(EPS1.GT.0.)IP0=IP0+1 72 CONTINUE GO TO 103 101 CONTINUE print*,'gone to 101' DO 107 J=2,NTOT EPS1=1. DIF=(FAIS(I,J)-FAIS(I,J-1))*EPS IF(DIF.GT.0.)EPS1=0. P(J)=FAIS(I,J)+EPS*(NREF*TPI*FLOAT(IP)+EPS1*TPI+IP0*TPI) IF(J.EQ.2)GO TO 104 Z1=ABS(P(J-1)-P(J-2)) Z2=ABS(P(J)-P(J-1)) zzs=z2-z1 c PRINT*,'Z1,Z2',tb(j),Z1,Z2,zzs IF(Z2.GE.Z1*0.8)GO TO 106 EPS1=1. P(J)=P(J)+TPI*EPS 106 CONTINUE c print*,'fais,p',tb(j),fais(i,j),p(j) 104 CONTINUE IP=IP+1 IF(EPS1.GT.0.)IP0=IP0+1 107 CONTINUE 103 CONTINUE DO 71 J=1,NTOT 71 FAIS(I,J)=P(J) PRINT*,'APRES LISSAGE' 777 FORMAT(' FAIS',6E12.5) 789 FORMAT(' Q0 Q1',3E12.5) C Q0=Q0-TPI*AINT(Q0/TPI) ip=0 eps=1. DO 77 J=1,NTOT DP(J)=XX(I,J)-FAIS(I,J) ccc modif 06/29/92: dp(j)=dp(j)-eps*tpi*ip if(j.eq.1)go to 775 eps=1. dff=dp(j)-dp(j-1) if(dff.lt.0.)eps=-1. print*,tb(j),xx(i,j),fais(i,j),ip,eps,dff,dp(j-1),dp(j) if(abs(dff).gt.abs(qq))then ip=ip+1 dp(j)=dp(j)-eps*tpi endif 775 continue c write(6,778)tb(j),xx(i,j),fais(i,j),dp(j),dff ccc end modif 778 FORMAT(' TB XX FAIS DP',F9.3,4E12.5) 77 CONTINUE DFI=YNTERP(TB,DP,TT1,NTOT,0) dfir=ynterp(tb,dp,tt2,ntot,0) WRITE(6,780)Q0,DFI,dfir 780 FORMAT(' Q0,DFI',3E12.5) NM=IDINT((Q0-DFI)/TPI) WRITE(6,118)NM 118 FORMAT(' NM=',I5) DO 73 J=1,NTOT c modif c DFI=FAIS(I,J)-XX(I,J) c DFI=-DFI dfi=dp(j) C WRITE(1,771)I,J,FAIS(I,J),XX(I,J),DFI,DIST(I) 771 FORMAT(2I5,3E12.3,F8.2) DISTA=DIST(I) FAC=DISTA/TB(J) DO 74 K=1,10 NK=NM-K+4 CV(J,K)=FAC/(DFI/TPI+NK) IF(CV(J,K).GT.50.OR.CV(J,K).LT.0.)CV(J,K)=0. 74 CONTINUE 73 CONTINUE WRITE(6,13)IST(I,1),IST(I,2),DIST(I),AZ(I) 1003 FORMAT(' PERIOD VPOBS VPMOD SCATT DOMGA VGOBS VGMOD SCATT *NUM') 13 FORMAT(1X,2A2,1X,'DISTANCE',F11.2,2X,'AZIMUTH',F8.2,F8.2,I5) c J=0 c DO 76 L=1,NT0 c IF(TK(L).GT.TB(1))GO TO 76 c DO 86 KK=1,NTOT c DIFF=ABS(TK(L)-TB(KK)) c IF(DIFF.LT.1.E-4)GO TO 87 86 CONTINUE c GO TO 76 87 CONTINUE c J=J+1 C PRINT*,'J L,K TB TK',J,L,K,TB(K),TK(L) c DO 88 K=1,10 88 CV(J,K)=CV(KK,K) c CTC(J)=CT(KK) c TK(J)=TK(L) 76 CONTINUE c NT0=J WRITE(12,13)IST(I,1),IST(I,2),DIST(I),AZ(I),AM(I),NTOT DO 761 J=1,NTOT WRITE(12,15)TB(J),(CV(J,K),K=1,10),CT(J) 761 CONTINUE 15 FORMAT(1X,'T=',F6.2,10F6.3,1X,'CT',F7.3) 70 CONTINUE 90 CONTINUE 5560 FORMAT(F8.2,3(1X,F6.2),1X,I3) 5570 FORMAT(10F7.4) WRITE(6,998)IPHAS 998 FORMAT(' IPHAS',I5) 999 CONTINUE CLOSE(7) CLOSE(13) CLOSE(12) CLOSE(8) CLOSE(15) stop end