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. modified 06/30/92 C CHARACTER*36 IPH CHARACTER*36 IDAT CHARACTER*80 IIN,INREF INTEGER*2 DATFIL(18),IST(50,3) INTEGER*2 ITITLE(30),OUTFIL(18) INTEGER*2 JIAK(2) DIMENSION TE1(4),TE2(4),er(4) COMMON/TDKP/TOUL(20),QI1(20),QI2(20),JOY,RAP(20) COMMON/COP/GLATE,GL1,GLATS,GL2 COMMON/DOCH/ DP(600),DON(1000),CU(600),GN1(5) COMMON/OUF/D(600),XM(6),TC(600),CC(600),P(600),CV(600,20) COMMON/MODD/TB(600),R(600),NRAD,NMOD,MOD(20),NTOT,EIP(600) *,NF(4),NLA(4) COMMON/BITS/PI,RN,VN,WN,GN,FL,FL1,FL2,FL3,SFL3,JCOM COMMON/DATA/ASS(5,600),FAIS(5,600),DIST(50),AZ(50) *,AM(50),DT(50) COMMON/EIGEN/AS(600),AP(600),AQ(600),CT(600),UU(600),Q(600) *,TT(600) COMMON/SPECTR/RER(5,600),AIMR(5,600),X(5,600),XX(5,600),XM0 EQUIVALENCE(OUTFIL,IPH) EQUIVALENCE(DATFIL,IDAT) DATA BIGG/6.6723D-11/,RHOBAR/5515.D0/,NMIN/5/ PI=3.14159265358979D0 RAD=180./PI TPI=2.*PI c c read filenames, epicentral data,period ranges etc from unit 5 c PRINT*,'CALCULATE PHASE VEL.(YES=1)' READ(5,*)IPHAS PRINT*,'ANTIPODAL DATA(YES=1)' READ(5,*)IPOD IF(IPOD.EQ.1)WRITE(6,115) PRINT*,'PHASEVEL FILENAME' READ(5,100)OUTFIL OPEN(12,FILE=IPH) print*,'MODEL FILE NAME' READ(5,'(a80)') IIN OPEN(7,FILE=IIN) print*,'MODEL TITLE' READ(5,1000)ITITLE 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*,'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) c C READ MODEL RADII from unit 7 c READ(7,200)(ITITLE(I),I=1,20) WRITE(6,200)(ITITLE(I),I=1,20) READ(7,112)IFANIS,IREF,IFDECK READ(7,113)NRAD,NIC,NOC READ(7,105)(R(I),I=1,NRAD) C NORMALIZE RN=R(NRAD) GN=PI*BIGG*RHOBAR*RN VN2=GN*RN VN=SQRT(VN2) WN=VN/RN C C read reference phase velocities from unit 11 C PRINT*,'REFERENCE PHASEVEL FILENAME' READ(5,'(a80)')INREF OPEN(11,FILE=INREF) READ(11,120)IPER WRITE(6,120)IPER READ(11,56)IMOD WRITE(6,56)IMOD READ(11,120)NC WRITE(6,120)NC READ(11,*)(TC(I),CC(I),I=1,NC) IF(IPOD.EQ.1)GO TO 98 c c read source parameters and calculate moment c print*,'STRIKE DIP SLIP DEP TAU XM' READ(5,*)STRK,DIP,SLIP,DEP,TAU,XM0 WRITE(6,555)STRK,DIP,SLIP,dep,tau,xm CALL MOMENT(STRK,DIP,SLIP,XM) 98 CONTINUE C C read data on unit 1more run.vit5 : first read title and periods of spectra C READ(15,1000)ITITLE JIAK(1)=ITITLE(1) JIAK(2)=ITITLE(2) WRITE(6,55)ITITLE WRITE(12,55)ITITLE READ(15,300)NTOT READ(15,102)(TB(I),I=1,NTOT) c c read data train by train c 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 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) I=I+1 GO TO 80 900 CONTINUE NSTAT=I-1 WRITE(12,5500)NOE,JIAK,IST(1,1),IST(1,2),STRK,DIP, *SLIP,DEP,TAU,XM0 WRITE(12,5510)GLATE,GL1,STLAT,STLON WRITE(12,*)NSTAT IF(IPOD.EQ.1)GO TO 63 C C ENTER EIGENFUNCTIONS and calculate source term C 1 CONTINUE CALL MODIN(JCOM,IOPT,DEP,DRIL,0.) CALL RADPAT(XM,NSTAT,STRK) DO 64 J=1,NTOT DO 64 I=1,NSTAT X(I,J)=RER(I,J)*RER(I,J)+AIMR(I,J)*AIMR(I,J) 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)) XX(I,J)=XX(I,J)-(TAU*TPI/TB(J)) c modif 06/30/92: fais(i,j)=xx(i,j)-fais(i,j) fais(i,j)=fais(i,j)-aint(fais(i,j)/tpi)*tpi 666 continue if(fais(i,j).ge.0.)go to 677 fais(i,j)=fais(i,j)+tpi go to 666 677 continue 64 CONTINUE 63 CONTINUE c c interpolate reference model to periods of spectra c DO 99 L=1,NTOT PER=TB(L) CT(L)=YNTERP(TC,CC,PER,NC,0) 99 CONTINUE PRINT*,'ENTER reference period' read(5,*)tt1 print*,'tt1',tt1 do l=1,ntot diff=abs(tt1-tb(l)) if(diff.lt.1.e-4)go to 47 enddo 47 l2=l+1 l1=l tt2=tb(l2) ct1=ct(l1) ct2=ct(l2) print*,'tt1 tt2 ct1 ct2',tt1,ct1,tt2,ct2,l1,l2 c c calculate phase velocities c DO 70 I=1,NSTAT DIS=DIST(I) IDS=INT(AM(I)+1.) 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 FAC=DIS/TT1 Q0=FAC/CT1*TPI Q1=DIS/TT2/CT2*TPI QQ=ABS(Q1-Q0) c c start unwrapping phase c c modif dfi=fais(i,l1) nm=int((q0-dfi)/tpi) print*,'nm=',nm 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 c IF((AM(I)-1.).LT.1.E-4)GO TO 101 if(dist(i).lt.15000.)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+nm*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 54 CONTINUE 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 fref= tpi*dist(i)/tb(j)/ct(j) 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+nm*tpi) IF(J.EQ.2)GO TO 104 Z1=ABS(P(J-1)-P(J-2)) Z2=ABS(P(J)-P(J-1)) izzs=(fref-fais(i,j))/tpi c PRINT*,'Z1,Z2',tb(j),Z1,Z2,zzs c IF(Z2.GE.Z1*0.8)GO TO 106 c EPS1=1. c P(J)=P(J)+TPI*EPS 106 CONTINUE print*,'fais,p',tb(j),fais(i,j),p(j),fref,izzs,z1,z2 104 CONTINUE IP=IP+1 IF(EPS1.GT.0.)IP0=IP0+1 107 CONTINUE 103 CONTINUE c c end unwrapping phase c DO 71 J=1,NTOT 71 dp(J)=P(J) cccccc c DO 77 J=1,NTOT c DP(J)=XX(I,J)-FAIS(I,J) c77 CONTINUE cccccc c c DFI=YNTERP(TB,DP,TT1,NTOT,0) dfir=ynterp(tb,dp,tt2,ntot,0) WRITE(6,780)Q0,DFI,dfir c NM=IDINT((Q0-DFI)/TPI) c WRITE(6,118)NM c c calculate 10 consecutive phase velocity curves c DO 73 J=1,NTOT DISTA=DIST(I) FAC=DISTA/TB(J) dfi=dp(j) DO 74 K=1,10 c NK=NM-K+4 nk=-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 c c end writing in output file c WRITE(6,13)IST(I,1),IST(I,2),DIST(I),AZ(I) 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 70 CONTINUE 90 CONTINUE WRITE(6,998)IPHAS 998 FORMAT(' IPHAS',I5) 999 CONTINUE CLOSE(7) CLOSE(13) CLOSE(12) CLOSE(8) CLOSE(15) 12 FORMAT(A2,2A1,1X,2I4,F10.3,F5.1,I6,2F10.3,F5.0,2F8.2) 13 FORMAT(1X,2A2,1X,'DISTANCE',F11.2,2X,'AZIMUTH',F8.2,F8.2,I5) 15 FORMAT(1X,'T=',F6.2,10F6.3,1X,'CT',F7.3) 55 FORMAT(36A2) 56 FORMAT(36A2) 100 FORMAT(18A2) 102 FORMAT(6(1X,E12.5)) 105 FORMAT(F8.0) 112 FORMAT(I1,I4,I2) 113 FORMAT(3I5) 115 FORMAT(' ANTIPODAL DATA') 118 FORMAT(' NM=',I5) 120 FORMAT(I5) 200 FORMAT(20A2) 300 FORMAT(I5) 555 FORMAT(' STRK DIP SLIP DEP TAU MO'/(6F9.3)) 778 FORMAT(' TB XX FAIS DP',F9.3,4E12.5) 779 FORMAT(' NREF=',I5,'FIREF',F10.3) 780 FORMAT(' Q0,DFI',3E12.5) 789 FORMAT(' Q0 Q1',3E12.5) 1000 FORMAT(30A2) 5500 FORMAT(I3,1X,2A2,1X,A2,A1,6(1X,F7.2)) 5510 FORMAT(4F8.2) stop end