C C PROGRAM FOCMEC C C Calculates acceptable focal mechanisms based on polarities C and (S/P) amplitude ratios. C Inspired by the Georgia Tech program with statistics C from Kisslinger's. Conventions are those in C Herrmann as well as Aki & Richards C C Code is written in Fortran 77 for a PDP 11/34 operating C under RSX version 4.1 C --------------------------------------------------------------- C C Seisan program FOCMEC modified to run original FOCMEC program (focmec_exe) C L. Ottemoller, January 2002 C C Part of program converted for use under pc/ms dos, ms-fortran 5.0 C R. Arvidsson, 1992. c c amplitudes c c amplitudes P on Z and R, SV on Z and R and SH on T. c all events can be used, local and global c AMPN means using first arrival, global or local c AMPG forces use of PG in local case c Q correction made based on travel times, for global use calcualted c for first P or S- phase found in file, for local used first c arrival. If no S-time, calculate for both cases using Vp/Vs. c For now Q hardwired to 100*f for local and tstart=1 and 4 for P c and S for global. c the screen output list all amplitude observations and Q-corrections to be c multiplied with the amplitudes to get correct amplitudes. c the screen output of ratios give ratios and amplitudes (not corrected) c used to form ratios, c and the free surface corrections is to be multipled by the the ratio. changes c c oct 92 by jh: adopted to SEISAN c nov 8 92 by jh: Fix a stroke, focmec.out in lowe case c nov 27 ; fix initialize minbadp c feb 14 93 ; change dip strike rake to strike dip rake in all files, c circles for C an D, strike dip rake on plot c sep 14 93 : version 3.0, why not before ?????? c 17 : add focinp, change to foc_prepare c nov 93 : latest bl routine, plot station option c mar 22 93 : bug with ---, az c may 5 94 : new hypocenter c july 28 : fix common block to include minbadp in all c sep 94 : incxrease dimentions c nov 7 : do not run foc_prepare if any argument c dec 8 : dimensions c feb 7 95 : bug CJAB(BGS)Jan95 : Install error and file handling. CJAB(BGS)Feb95 : Tidy up some loose ends. c mar 28, 94 by jh : set window size correctly c apr 28 95 jh : add angle calculation as an option, c for composite c may 18 : bug c nov 23, 95 jh : make sure only P-phases are used c feb 13m 95 : only assume full weight polarities, however interanlly c full weight polarities still count for 2 and + anf - c for 1, so + and - can still be used. c nov 7 96 : adjust for large number of polarities when comp. f. sol c nov 24 : bug in --------- c mar 1 98 : enable plotting of several soluttion sform s-file, type 'OF' c mar 1 99 : ----------- version 7.0 check ---------------------------- C station to 5 chars c mar 23 99 bmt : include winplot.inc c apr 20 bmt : include winplot3.inc c sep 22 jh : also possible to plot faultplane sol. without a location c oct 28 jh : weighted phases now used, error with 5 char stations c nov 3 jh : when doing angles, do no tstop after one calculation c as again, mostly so result do not disappear on pc c aug 16 02 lot : changed SH(T) to SG(T) c aug 19 02 lot : fixed bug with increment c apr 29 03 lot : change in sterpol to ignore replot c jul 1 03 lot : delete routines that are no longer used c changed call to FMREPS c jul 15 jh : add momten to calls of pttpin and dsrin c jun 20 05 jh : plt to eps c dec 10 06 jh : fix so fps can be plotted without readings c Dec 18 08 jh : fix so all 3 amplitude ratios can be used, rewrite amplitude c part to include global case and put in attenuation, c fix a few bugs c Feb 6 09 jh : fix bug c feb 18 : group R as L C---------------------------------------------------------------- C C Arthur Snoke Virginia Tech 1984 C 23 October 1984: Added the possibility for using SH C amplitudes instead of only SV on the vertical. C Also, if the input line key is V (for SV) or H C (for SH), the program expects to find a polarity C for the S and a comment following the log(ratio). C Entering R for the key will work as before: assumes C the ratio is for SV on vertical and there is no C polarity or comment. The output now includens the C theoretical S polarity, but it does not use it in C the search. Another small change which only acts C to speed up the search is that for a B axis plunge C of 90 degrees there is only one trend tried. C The SH amplitude case has not been checked yet. C $ December 1984: Eliminated double solutions when plunge C of B is 0 degrees. (Only allowed trend to vary C between 0 and 180 degrees. This could affect C searches for small parts of the focal sphere.) C C MARCH 1990: CHANGED TO PC-DOS R.N.ARVIDSSON C C- include 'winplot1.inc' c include 'pc.inc' c PROGRAM FOCMEC C C SEISAN library/JAB inclusions... c -------------------------------- C include 'libsei.inc' ! Library definitions & data defns. include 'seidim.inc' ! array dimentions include 'seisan.inc' include 'seiplot.inc' ! plotting parameters c external & sei open, ! File open handler. & sei close, ! & closure. & sei clen, ! & length of string. & sei integer, ! Integer decoder. & sei code ! Error condition handler. integer sei clen, ! & function. & sei integer ! Ditto. integer code, ! Condition. & NSOL ! # solutions found. logical b_exit ! Force exit routine. c integer read1, ! Read unit1. & write1, ! Write unit1. & write2, ! Ditto 2. & write3, ! & 3. & write4, ! & 4. & write5, ! & 4. & write6 ! & 6 common /foc_unit$ / & read1, ! Read unit1. & write1, ! Write unit1. & write2, ! Ditto 2. & write3, ! & 3. & write4, ! & 4. & write5 ! & 5. c logical b_flag ! Flag existance or is a dummy!! C C ------- End of details ------- C C- PARAMETER (MAX = max_data) C---------! Max no. of polarities or polarities PARAMETER (MAX2 = 2*max) C----------! Twice MAX REAL*4 XYZ(6,MAX2) C----------! R & Theta station coordinates REAL*4 PSTATN(MAX) C----------! Stations with polarity data REAL*4 RSTATN(MAX) C----------! Stations with amplitude-ratio data C X = north, Y east, Z down INTEGER*2 POLRTY(MAX) C----------! 1 for C, -1 for D INTEGER*2 NPOL,NRATIO C----------! Number of Polarities read in INTEGER*2 KEYPOL(MAX) C----------! Keys polarity data relative to input REAL*4 LOGRAT(MAX) C----------! Log10 of (SV/P) or (SH/P) including C free surface correction, etc. INTEGER*2 NRAT C----------! Number of ratios read in INTEGER*2 KEYRAT(MAX) C-----------! Keys ratios data relative to input CHARACTER SVSH(2,max) C-----------! First element V or H for SV or SH, DIMENSION = MAX C Second element polarity: F or B for SV C L or R for SH (back to station) INTEGER*2 MAXSOL C------------! Exit after this many accept. aols. INTEGER*2 NERRP,NERRA ! allowed errors for pol and amp ratio REAL max_amp_error ! maximum allowed amplitude ratio error INTEGER MINERRP C-------------! Allowed number of polarity errors REAL*4 BADPOL(MAX) C------------! Stations with a polarity error INTEGER*2 NERRR C------------! Allowed number of ratio errors REAL*4 ERRRAT C------------! Maximum allowed Log10 ratio error REAL*4 THERAT(MAX) C------------! Theoretical Log10 of (SV/P) INTEGER*2 WTRAT(MAX) C------------! 1 if that ratio used, 0 otherwise REAL*4 VPVS3,VPVS C------------! Cube of P/S velocity ratio REAL*4 BTMIN,BTDEL,BTMAX C------------! Search range in B axis trend REAL*4 BPMIN,BPDEL,BPMAX C------------! Search range in B axis plunge REAL*4 AAMIN,AADEL,AAMAX C------------! Search range in A axis angle C measured from B trend C COMMON /FOCM/ PSTATN,RSTATN,XYZ,RD,POLRTY,NPOL,LOGRAT, 1 NRAT,MAXSOL,NERRP,NERRR,ERRRAT,THERAT,WTRAT,VPVS3, 2 BTMIN,BTDEL,BTMAX,BPMIN,BPDEL,BPMAX,AAMIN,AADEL,AAMAX, 3 BADPOL,KEYRAT,KEYPOL,minbadp COMMON /FOCCHAR/ SVSH logical stat_plot ! plot station codes if true real iwid character*1 xpol,cchar,type,exp character*80 text,data(max_data) character chr_solution *(80) ! Selected solution. integer deldeg ! search increment c -- tektronix options integer nars ! number of arguments character*80 arg(5) ! arguments C--- CHARACTER answer*1 character*80 editor C --- LOGICAL merge real xx(4) ! store strike,dip and rake, number of bad polarities logical sun,pc,linux C -------------- c include 'version.inc' include 'winplot2.inc' c c print version c out_version_date='July 23, 2001' if (version_new) out_version_date=version_date call print_ver call computer_type(sun,pc,linux) call get_seisan_def call get_editor(editor) c c Initialise... c ============= c code = e_ok$ ! Local condition. b_exit = .false. ! No forced exit. nsol = 0 ! No solutions found. write(*,*) write(*,*)'============ FOCMEC ============' read1 = 0 ! & unit. write1 = 0 ! Ditto. write2 = 0 ! Ditto. write3 = 0 ! Ditto. write4 = 0 ! Ditto. write5 = 0 ! Ditto. c chr_solution = ' ' ! Empty. c c Get command line arguments... c ----------------------------- c call get_arguments(nars,arg) c c prepare input from print.out from hyp unless there is an argument c which is either a for angle calculations or c for composite solution c if(nars.gt.0) then if(arg(1)(1:1).eq.'a'.or.arg(1)(1:1).eq.'A') then call angles stop endif if(arg(1)(1:1).ne.'c'.and.arg(1)(1:1).ne.'C') then write(6,*)' Invalid argument, use:' write(6,*)' a: Calculate angles' write(6,*)' c: Composite fault plane solution, from EEV' write(6,*)' None: Fault plane solution after using HYP' stop endif else call foc_prepare(npol,nratio,vpvs) endif c c both hard copy file and screen plot c plotoption=1 c--- hard copy type is postscript hctype=1 c-- set xwindow sixe wsize=80 c-- no graphics yet disp_open=0 c--- Merge polarities, t,p, and faultplane merge=.true. c--- Minimum bad polarities minbadp=1000 c--- size of circle rcircle=330.0 c--- no new solution inew=0 c c Open files... c ============= c focmec.inp is input in nordic format c call sei open( old$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.inp', ! File name to open. & read1, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). c call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.lst', ! File name to open. & write1, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.out', ! File name to open. & write2, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). call sei open( scratch$, ! Open (default stop on error). & ' ', ! No prompt. & 'scratch1.out', ! File name to open (n/a). & write3, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). c call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'scratch2.out', ! File name to open (n/a). & write4, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). c call sei open( scratch$, ! Open (default stop on error). & ' ', ! No prompt. & 'scratch3.out', ! File name to open (n/a). & write5, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). c c make sure focmec.out is empty c write(write2,*,iostat=code) call sei code(fort$,code,write2,b_flag) ! Process outcome. rewind(write2,iostat=code) call sei code(fort$,code,write2,b_flag) ! Process outcome. c c-------------------------------------------------------------------- c back here for decisions CJAB(BGS)Jan95 ....it seems logical to prompt user to save any CJAB(BGS)Jan95 solutions at this point, rather than when quiting CJAB(BGS)Jan95 the routine. The user can then plot the solution CJAB(BGS)Jan95 without having to come back in to the program c--------------------------------------------------------------------- C Treat any selected solutions... c =============================== c write in input file if so chosen and if a new solution has been selected c 100 if(inew.eq.1) then write(6,*)' Save solution (y/n)' read(5,'(a1)') answer if(answer.eq.'y'.or.answer.eq.'Y') then inew = 0 ! Re-set. rewind(read1,iostat=code) call sei code(fort$,code,read1,b_flag) ! Process outcome. call indata(read1,nstat,nphas,nhead, & nrecord,type,exp,data,id) c c check if any solution from before, if so overwrite, else add c ifoc=0 chr_solution(79:80)=' F' chr_solution(71:78)='FOCMEC ' do i=2,nhead if(data(i)(79:80).eq.' F') then ifoc=1 data(i)=chr_solution endif enddo if(ifoc.eq.0) then do i=nrecord,nhead,-1 data(i+1)=data(i) enddo data(nhead)=chr_solution nrecord=nrecord+1 endif rewind(read1,iostat=code) call sei code(fort$,code,read1,b_flag) write(read1,'(a80)',iostat=code)(data(i),i=1,nrecord) call sei code(fort$,code,read1,b_flag) endif c c If exit has been forced, then do so... c -------------------------------------- c if( b_exit ) goto 999 ! & exit. endif c c Main menu.... c ============= c 150 write(6,*) write(6,*) ' Stop (0)' write(6,*) ' Plot saved solution(s) (1)' write(6,*) ' Plot new solutions (2)' write(6,*) ' Plot selected solution (3)' write(6,*) ' Find new solutions (4)' write(6,*) ' -1, -2, -3 also plot station' c read(*,'(a)') text ! Get reply. iselect = sei integer( text, code ) ! & decode it. c c Process the outcome... c ====================== c Invalid alpha text... c ===================== c if( code .ne. e_ok$ ) then ! Unable to decode. write(*,*) write(*,*)'**** WARN: invalid option ...try again' goto 100 ! & menu. c c Numeric options... c ================== c else ! Otherwise. stat_plot = iselect .lt. 0 ! Plot stations?. if( iselect .ne. -4 .and. ! & re-set. & iselect .lt. 0 ) then ! Ensure positive. iselect = -iselect ! end if ! end if ! c c Finished... c =========== c If save solution flagged, then go back and try again... c ------------------------------------------------------- c if( iselect .eq. 0 ) then ! Exit. if( inew .eq. 1 ) then ! Try again. b_exit = .true. ! Flag exit. goto 100 ! c else ! Otherwise exit. goto 999 ! end if ! c c Plot saved solutions... c ======================= c else if( iselect .eq. 1 ) then ! Plot saved solutions. rewind(read1,iostat=code) ! Rewind all files. call sei code(fort$,code,read1,b_flag) ! Process outcome. rewind(write1,iostat=code) ! call sei code(fort$,code,write1,b_flag) ! Process outcome. rewind(write2,iostat=code) ! call sei code(fort$,code,write2,b_flag) ! Process outcome. c rewind(write4,iostat=code) ! call sei code(fort$,code,write4,b_flag) ! Process the outcome. c iold=0 call indata(read1,nstat,nphas,nhead,nrecord, & type,exp,data,id) c c select prime fault plane solution and other solutions 'OF' c write(write4,'(a)')' Strike Dip Rake' do i=1,nhead if(data(i)(79:80).eq.' F'.or.data(i)(79:80).eq.'OF') then if(iold.eq.0) then c do k=1,7 c do k=1,10 ! lot 18.01.2002 c lot Feb 4, 2002 write(write4,'(a)',iostat=code) data(i) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c enddo else write(write4,'(a)',iostat=code) data(i) endif c call sei code(fort$,code,write4,b_flag) ! Process the outcome. iold=1 endif enddo rewind(write4,iostat=code) c c check if any old solution, else go back to main menu c if( iold .ne. 1 ) then write(*,*) write(*,*)'**** WARN: no solution in S-file' goto 100 end if c c Plot new solutions... c --------------------- c ...there must be at least 7 lines in focmec.out file c else if( NSOL .gt. 0 .and. ! New solutions available. & iselect .eq. 2 ) then ! Plot new solutions. rewind(read1,iostat=code) ! call sei code(fort$,code,read1,b_flag) ! Process outcome. rewind(write1,iostat=code) ! call sei code(fort$,code,write1,b_flag) ! Process outcome. rewind(write2,iostat=code) ! call sei code(fort$,code,write2,b_flag) ! Process outcome. rewind(write4,iostat=code) ! call sei code(fort$,code,write4,b_flag) ! Process the outcome. c do i=1,9000 read(write2,'(a)',iostat=code) text call sei code(fort$,code,write2,b_flag) ! Process outcome. if( b_flag ) goto 95 ! End of file. write(write4,'(a)',iostat=code) text call sei code(fort$,code,write4,b_flag) ! Process the outcome. enddo c 95 rewind(write2,iostat=code) call sei code(fort$,code,write2,b_flag) ! Process outcome. c rewind(write4,iostat=code) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c c Invalid file information... c if( i .lt. 7 ) then ! Must be 7 lines. write(*,*) write(*,*)'**** WARN: no new solutions available' goto 100 end if c c No solutions available... c ------------------------- c else if( iselect .eq. 2 ) then ! write(*,*) write(*,*)'**** WARN: no new solutions available' goto 100 ! Back to menu. c c Plot selected solution... c ------------------------- c else if( NSOL .gt. 0 .and. ! Selection available. & iselect .eq. 3 ) then ! Plot selected solution. rewind(write4,iostat=code) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c rewind(write3,iostat=code) call sei code(fort$,code,write3,b_flag) ! Process outcome. c c Loop entries, must be at least 7... c c do i=1,7 c do i=1,10 ! lot 18.01.2002 write(write4,'(a)')' Strike Dip Rake' text = 'xxx' do while (text.ne.' Strike Dip Rake') read(write3,'(a)',iostat=code) text ! enddo read(write3,'(a)',iostat=code) text ! call sei code(fort$,code,write3,b_flag) ! Process outcome. c c ... end of file... c if( b_flag ) then ! End of file. write(*,*) ! write(*,*) & '**** WARN: no selected solution available' goto 100 ! Menu. end if ! c c Otherwise write out... c write(write4,'(a)',iostat=code) text call sei code(fort$,code,write4,b_flag) ! Process the outcome. c enddo c rewind(write3,iostat=code) call sei code(fort$,code,write3,b_flag) ! Process outcome. c rewind(write4,iostat=code) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c c No solutions available... c ------------------------- c else if( iselect .eq. 3 ) then ! write(*,*) write(*,*)'**** WARN: no selected solution available' goto 100 ! Back to menu. c c Get more solutions... c ===================== c else if( iselect .eq. 4 ) then ! Get more solutions rewind(read1,iostat=code) call sei code(fort$,code,read1,b_flag) ! Process outcome. rewind(write1,iostat=code) call sei code(fort$,code,write1,b_flag) ! Process outcome. rewind(write2,iostat=code) call sei code(fort$,code,write2,b_flag) ! Process outcome. c c input data for iterations c c CALL FOCINP c c c calculate focal mechanisms c c CALL SRCHFM( NSOL ) ! Get # solutions. c c c create input file to run focmec, lot 16.01.2002 c write(*,'(a,i3,a)') ' There are ',npol,' polarity readings ' 135 continue write(*,'(a)') ' Maximum number of allowed polarity errors or ' & // '-1 to show best solutions only ' c read(*,*,err=135) nerrp read(*,'(a)') text ! Get reply. nerrp = sei integer( text, code ) ! & decode it. minerrp=nerrp 136 continue if (nratio.gt.0) then write(*,'(a,i3,a)') ' There are ',nratio,' amp ratio readings' write(*,'(a)') ' Maximum number of allowed amplitude ' // & 'ratio errors ' c read(*,*,err=136) nerra read(*,'(a)') text ! Get reply. nerra = sei integer( text, code ) ! & decode it. if (nerra.gt.nratio) nerra=nratio write(*,'(a)') ' Maximum amplitude ratio error, ' & // ' return for default of .2 ' c read(*,*,err=136) max_amp_error read(*,'(a)') text ! Get reply. read(text,'(f5.2)',err=136) max_amp_error if (max_amp_error.eq.0.) max_amp_error=.2 endif if (nerrp.eq.-1) then ! search for solutions with nerrp, ! max number of solution is 99999, ! automatically find best solutions nerrp=int(sqrt(float(npol))+.5) minerrp=99999 endif c c since no longer half polarities counted, multiply by 2 c c nerrp=nerrp*2 ! taken out, lot 18.01.2002 IF (NERRP .GT. NPOL) NERRP = NPOL write(6,*)'Degree increment in search' read(5,*) deldeg c c write focmec command input c call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.run', ! File name to open (n/a). & write6, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). write(write6,'(a)') 'fm_new.out' write(write6,'(a)') '-----' if (linux) & write(write6,'(a)') ' ' c write(write6,'(a)') '/* Comment is previous line: Input file ' c & //'for focmec is next ' write(write6,'(a)') 'focmec.dat' write(write6,'(a)') ' correct file [y] ' write(write6,'(a)') ' relative weighting..[n] ' if (nerrp.le.9) then write(write6,'(i1,a)') nerrp, & ' allowed P polarity erors..[0] ' elseif (nerrp.le.99) then write(write6,'(i2,a)') nerrp, & ' allowed P polarity erors..[0] ' else write(write6,'(i3,a)') nerrp, & ' allowed P polarity erors..[0] ' endif if (nratio.gt.0) then write(write6,'(f4.2,4x,a)') vpvs,'vp/vs ratio' c write(write6,'(8x,a)') 'maximum allowed log10 of ratio' write(write6,'(f4.2,4x,a)') max_amp_error, & 'maximum allowed log10 of ratio' write(write6,'(i1,7x,a)') nerra,'allowed amplitude errors' write(write6,'(f4.2,4x,a)') .05, & 'lower limit P cutoff' write(write6,'(f4.2,4x,a)') .15, & 'lower limit S cutoff' endif if (minerrp.eq.99999) then ! in this case find best solutions write(write6,'(a)') & '99999 exit after this many acceptable solutions...[100] ' else write(write6,'(i4.4,a)') focmec_maxsol, & ' exit after this many acceptable solutions...[100] ' endif write(write6,'(a)') ' minimum B trend [0] ' write(write6,'(i3.3,a)') deldeg,' B increment [5] ' write(write6,'(a)') ' maximum B trend [355] ' write(write6,'(a)') ' min B plunge..[0] ' write(write6,'(i3.3,a)') deldeg,' increment [5] ' write(write6,'(a)') ' maximum..[90] ' write(write6,'(a)') ' minimum A angle..[0] ' write(write6,'(i3.3,a)') deldeg,' increment [5] ' write(write6,'(a)') ' maximum [85] ' call sei close( close$, write6, code ) ! Close plot file. c c option to edit focmec.run c write(*,*) ' Do you want to edit focmec.run (y/n) ? ' answer=' ' read(5,'(a1)') answer if (answer.eq.'y'.or.answer.eq.'Y' ) then text=editor(1:seiclen(editor))// ' focmec.run' write(*,*) text(1:seiclen(text)) call systemc(text,seiclen(text)) endif c c run focmec, lot 16.01.2002 c text='focmec_exe < focmec.run > focmec.log ' write(*,*) text(1:seiclen(text)) call sei close( close$, write1, code ) ! Close plot file. call sei close( close$, write2, code ) ! Close plot file. call systemc(text,seiclen(text)) call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.lst', ! File name to open. & write1, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.out', ! File name to open. & write2, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). c c copy output from fm_new.out to focmec.out c call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'fm_new.out', ! File name to open (n/a). & write6, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). if(minerrp.ne.99999) goto 191 c c find solutions with minimum polarity errors c i=0 text='x' c read headers do while (text(1:23).ne.' Dip Strike Rake') read(write6,'(a)') text enddo 190 continue read(write6,'(a)',end=191) text c i=i+1 c if (i.gt.9) then read(text,'(3f8.2,f11.1)') xx(1),xx(2),xx(3),xx(4) if (xx(4).le.float(minerrp)) then minerrp=int(xx(4)) endif c endif goto 190 191 continue rewind(write6,iostat=code) write(*,*) write(*,*)'...Minimum number of bad fits are ',minerrp write(*,*) i=0 c skip header do while (text(1:23).ne.' Dip Strike Rake') c i=i+1 read(write6,'(a)') text if (text(1:23).eq.' Dip Strike Rake') then write(*,'(a)') ' Strike Dip Rake Pol: P '// & 'SV SH Rat Err RMS RErr RErr (All)' write(write2,'(a)') & ' Strike Dip Rake Pol: P '// & 'SV SH Rat Err RMS RErr RErr (All)' else write(*,'(a)') text write(write2,'(a)') text endif enddo c if (i.eq.9) c text=' Strike Dip Rake ' c write(write2,'(a)') text c write(*,'(a)') text i=0 192 continue read(write6,'(a)',end=193) text read(text,'(3f8.2,f11.1)') xx(1),xx(2),xx(3),xx(4) if (xx(4).le.float(minerrp)) then write(text(1:36),'(3f10.4,i6)') xx(2),xx(1),xx(3),int(xx(4)) else text='skip' endif if (text.ne.'skip') then write(write2,'(a)') text write(*,'(a)') text endif goto 192 193 continue call sei close( close$, write6, code ) ! Close plot file. c c get number of new solutions c rewind(write2,iostat=code) nsol=-1 195 continue read(write2,'(a)',end=196) text if (text(1:20).eq.' Strike Dip') then nsol = 0 elseif (nsol.ne.-1) then nsol=nsol+1 endif goto 195 196 continue write(*,*)' There are ',nsol,' acceptable solutions. ' if (nsol.ge.100) then write(*,*) & ' NOTE: The maximum of 100 solutions has been reached' & // ' and the search is stoped ! ' endif rewind(write2,iostat=code) c write(*,*) c write(*,*)'...Minimum number of bad fits are ',minbadp/2 c c back to main menu c goto 100 c c Invalid option... c ================= c else write(*,*) write(*,*)'**** WARN: invalid option ...try again' goto 100 endif c c--------------------------------------------------------------------- c start of plotting section c--------------------------------------------------------------------- c 200 continue c c rewind all files before continuing c rewind(read1,iostat=code) call sei code(fort$,code,read1,b_flag) ! Process outcome. rewind(write1,iostat=code) call sei code(fort$,code,write1,b_flag) ! Process outcome. rewind(write2,iostat=code) call sei code(fort$,code,write2,b_flag) ! Process outcome. rewind(write3,iostat=code) call sei code(fort$,code,write3,b_flag) ! Process outcome. c rewind(write4,iostat=code) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c rewind(write5,iostat=code) call sei code(fort$,code,write5,b_flag) ! process the outcome. c c-- open plot file c call sei open( unknown$, ! Open (default stop on error). & ' ', ! No prompt. & 'focmec.eps', ! File name to open. & plotunit, ! On unit. & b_flag, ! Flag exists. & code ) ! Condition (n/a). c c open plotting c call open_display c c scale tec plot to postscript c write(plotunit,*,iostat=code)' 1.0 0.55 scale' call sei code(fort$,code,plotunit,b_flag) ! Process outcome. c CALL STERMEC(answer,iwid,merge,rcircle) CALL STERPOL(answer,IWID,merge,rcircle,stat_plot) c c Give prompt... c call xset_color(color_prompt) if( iselect .eq. 2 ) then text = 'Enter Q to quit or' call tchars(text,seiclen(text),670.0,700.0) text = 'press to continue or' call tchars(text,seiclen(text),670.0,680.0) text = 'move the cursor and enter P or T' call tchars(text,seiclen(text),670.0,660.0) text = 'to select your preferred solution.' call tchars(text,seiclen(text),670.0,640.0) else text = 'Enter Q to quit or' call tchars(text,seiclen(text),670.0,700.0) text = 'press to continue.' call tchars(text,seiclen(text),670.0,680.0) end if c c make sure plot is finished c call tsend c c get up cursor to select solution by selecting a p or t axis c 250 continue call xscursr(ichar,x,y) call tsend c c check if P, T or Q c cchar=char(ichar) if(cchar.eq.'p') cchar='P' if(cchar.eq.'t') cchar='T' if(cchar.eq.'q') cchar='Q' if(cchar.eq.'r') cchar='R' c if( cchar .eq. 'Q' ) then call clear_to_alpha call close_post call sei close( close$, plotunit, code ) ! Close plot file. goto 999 ! Finished. c c replot, go back to xscursor, lot April 2003 c elseif ( cchar .eq. 'R' .or. & cchar .eq. '#' ) then ! lot 01/07/2003 goto 250 elseif( cchar .ne. 'P' .and. ! Back to menu. & cchar .ne. 'T' ) then call clear_to_alpha call close_post call sei close( close$, plotunit, code ) ! Close plot file. if( iselect .eq. 1 ) then ! Do not prompt for save. goto 150 ! else ! Otherwise prompt if flagged! goto 100 ! end if ! c c check if a solution has been selected, if so find it c by finding the minimum distance to one of the tek coordinates c P or T axis, stored in unit "write5", a scratch file c else ! Solution selected. rewind(write5,iostat=code) call sei code(fort$,code,write5,b_flag) ! process the outcome. c dmin=10000.0 do i=1,1000 read(write5,'(a1,2f8.1)',iostat=code)xpol,xx1,yy1 call sei code(fort$,code,write5,b_flag) ! process the outcome. if( b_flag ) goto 900 c if(cchar.eq.xpol) then dist=sqrt((x-xx1)**2+(y-yy1)**2) if(dist.lt.dmin) then imin=i dmin=dist endif endif enddo 900 continue inew=1 c c now find corresponding fault plane solution, first calculate c fault plane solution number ipos c ipos=imin if(cchar.eq.'T') ipos=ipos+1 ipos=ipos/2 rewind(write2,iostat=code) call sei code(fort$,code,write2,b_flag) ! Process outcome. c c read forward to solution, there are 6 headers c c do i=1,6+ipos c do i=1,9+ipos ! lot, 17.01.2002 c read(write2,'(a)',iostat=code) chr_solution c call sei code(fort$,code,write2,b_flag) ! Process outcome. c enddo chr_solution = 'xxx' do while (chr_solution(1:29).ne. & ' Strike Dip Rake') read(write2,'(a)',iostat=code) chr_solution call sei code(fort$,code,write2,b_flag) ! Process outcome. enddo do i=1,ipos read(write2,'(a)',iostat=code) chr_solution call sei code(fort$,code,write2,b_flag) ! Process outcome. enddo c c write solution out in scratch file for future plotting c c do i=1,7 c do i=1,10 ! lot 17.01.2002 c write(write3,'(a)',iostat=code) chr_solution c call sei code(fort$,code,write3,b_flag) ! Process outcome. c enddo write(write3,'(a)',iostat=code) & ' Strike Dip Rake ' call sei code(fort$,code,write3,b_flag) ! Process outcome. write(write3,'(a)',iostat=code) chr_solution call sei code(fort$,code,write3,b_flag) ! Process outcome. c c Close PostScript and go back to alpha window... c call xmovabs(1.0,1.0) call clear_to_alpha call close_post call sei close( close$, plotunit, code ) ! Close the plot file. c c back to main menu c goto 100 end if ! c c Finished... c =========== c 999 write(*,*) write(*,*) ' Plot file is called focmec.eps' write(*,*) ' Input for GMT pspolar is called pspolar.inp ' c c Close all files before leaving... c ================================= c call sei close(delete$,write3, code ) ! Delete scratch. call sei close(delete$,write4, code ) ! Delete scratch. call sei close(delete$,write5, code ) ! Delete scratch. call sei close(close$+all$, 0, code ) ! Close everything. stop end C+ C SUBROUTINE DPSTRK(A,N,ANGS) C C C CALCULATES DIP, STRIKE AND RAKE (ANGS) FROM A AND N C- SUBROUTINE DPSTRK(A,N,ANGS) REAL N(3) DIMENSION A(3),ANGS(3) PI = 4.0*ATAN(1.0) IF (.NOT.(N(1) .EQ. 0.0 .AND.N(2) .EQ. 0.0)) THEN ANGS(2) = ATAN2(-N(1),N(2)) ELSE ANGS(2) = 0.0 END IF IF (ABS(SIN(ANGS(2))) .LT. 0.1) GO TO 100 ANGS(1) = ATAN2(-N(1)/SIN(ANGS(2)),-N(3)) GO TO 200 100 ANGS(1) = ATAN2(N(2)/COS(ANGS(2)),-N(3)) 200 A1 = A(1)*COS(ANGS(2)) + A(2)*SIN(ANGS(2)) IF (ABS(SIN(ANGS(1))) .LT. 0.1) GO TO 300 ANGS(3) = ATAN2(-A(3)/SIN(ANGS(1)),A1) GO TO 400 300 A2 = A(1)*SIN(ANGS(2)) - A(2)*COS(ANGS(2)) ANGS(3) = ATAN2(A2/COS(ANGS(2)),A1) 400 IF (ANGS(1) .GE. 0.0) GO TO 500 ANGS(1) = ANGS(1) + PI ANGS(3) = PI - ANGS(3) IF (ANGS(3) .GT. PI) ANGS(3) = ANGS(3) - 2*PI 500 IF(ANGS(1).LE.0.5*PI) GO TO 600 ANGS(1)=PI-ANGS(1) ANGS(2)=ANGS(2)+PI ANGS(3)=-ANGS(3) IF (ANGS(2) .GE. 2*PI) ANGS(2) = ANGS(2) - 2*PI 600 IF (ANGS(2) .LT. 0.0) ANGS(2) = ANGS(2) + 2.0*PI RETURN END c c******************SUBROUTINE INVERT********************* c subroutine invert(mat) c c inverts a 3x3 matrix c written by S.Ihnen, November 1980 c double precision mat(3,3),xmat(3,3) double precision a,b,c,det c c define elements of determinant c a=mat(1,1)*(mat(2,2)*mat(3,3)-mat(2,3)*mat(3,2)) b=mat(1,2)*(mat(2,1)*mat(3,3)-mat(3,1)*mat(2,3)) c=mat(1,3)*(mat(2,1)*mat(3,2)-mat(3,1)*mat(2,2)) det=a-b+c c c construct matrix of cofactors of original matrix c xmat(1,1)=mat(2,2)*mat(3,3)-mat(3,2)*mat(2,3) xmat(1,2)=-mat(2,1)*mat(3,3)+mat(3,1)*mat(2,3) xmat(1,3)=mat(2,1)*mat(3,2)-mat(2,2)*mat(3,1) xmat(2,1)=-mat(1,2)*mat(3,3)+mat(3,2)*mat(1,3) xmat(2,2)=-mat(3,1)*mat(1,3)+mat(3,3)*mat(1,1) xmat(2,3)=-mat(1,1)*mat(3,2)+mat(3,1)*mat(1,2) xmat(3,1)=mat(1,2)*mat(2,3)-mat(1,3)*mat(2,2) xmat(3,2)=-mat(1,1)*MAT(2,3)+mat(2,1)*mat(1,3) xmat(3,3)=mat(1,1)*mat(2,2)-mat(2,1)*mat(1,2) c c inverse of original matrix is just the transpose of c the adfoint matrix divided by the dterminant of the c original matrix c do 100 i=1,3 do 99 j=1,3 mat(i,j)=xmat(j,i)/det 99 continue 100 continue return end SUBROUTINE STERMEC(device,iwid,merge,rcircle) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Program to draw the sterographic projection, for focal mechan. C for a number of different solutions. C C R. Arvidsson. 1990-04-27 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c IMPLICIT INTEGER*2 (V,j-n) C --- include 'seiplot.inc' c c Libsei details... c ================= c include 'libsei.inc' ! Library definitions & data defns. external sei code ! Error condition handler. integer code ! Condition. logical b_flag ! Flag!! c integer read1, ! Read unit1. & write1, ! Write unit1. & write2, ! Ditto 2. & write3, ! & 3. & write4, ! & 4. & write5 ! & 5. common /foc_unit$ / & read1, ! Read unit1. & write1, ! Write unit1. & write2, ! Ditto 2. & write3, ! & 3. & write4, ! & 4. & write5 ! & 5. c c End of list... c ============== c real IX, IY, R, IX3(4), x(1000), 1 IDEVH, IN1, IN2, 2 XBTP(2), IXBTP, IYBTP, ixb(4), iwid C --- CHARACTER WHET(4)*1, 1 DUMMY*80, CBTP(3)*1, device*1 C --- c REAL ANBTP(6), ANGS(3), ANGS2(3), PTTP(4), DIP(200), STRIKE(200), c 1 RAKE(200), rcircle REAL ANBTP(6),ANGS(3),ANGS2(3),PTTP(4),DIP(1000),STRIKE(1000), 1 RAKE(1000),rcircle REAL*8 xo, yo, xp, yp, fi, delta, ox, alfa, az, ax, btp(2), 1 zo, alfa2, fi2 C --- REAL*4 MOMTEN(6) LOGICAL AN, PT, RIGHT, MT logical merge C --- c common /profil/ merge, rcircle C --- DATA WHET /'N','S','E','W'/ DATA CBTP /'B','T','P'/ DATA ZERO/0/,ONE/1/ c c Initialise... c ============= c code = e_ok$ ! Local condition. c C --- PI = 4.*atan(1.) C --- RIGHT=.TRUE. c call xset_color(color_frame) C -------Draw circle IX = 500 IY = 350 IX3(1)=IX IX3(2)=IY if (rcircle .lt. 0.0001 ) then R = 150 else r= rcircle endif IN1= 20 IN2 = IY + RCIRCLE*1.6 CALL FCIRCL(R, IX,IY) C --- Circle for B,P and T axes if ( merge ) then ixbtp = ix iybtp = iy else ixbtp = IX +2.25*RCIRCLE iybtp = iy endif cjens CALL FCIRCL(R,IXBTP,IYBTP) C --- Label the circle IX3(1)= IX IX3(2)= IY + R + 15 call xchars(whet(1),1,ix3(1),ix3(2)) c ISTAT = VGTEXT(IDEVH,IX3(1),IX3(2),WHET(1)) IX3(2)= IY + R IX3(3)= IX3(1) IX3(4)= IX3(2) - 5 call xmovabs(ix3(1),ix3(2)) call xdrwabs(ix3(3),ix3(4)) c ISTAT = VPLINE(IDEVH,2,IX3) IX3(2)= IY - R IX3(4)= IX3(2) + 5 call xmovabs(ix3(1),ix3(2)) call xdrwabs(ix3(3),ix3(4)) c ISTAT = VPLINE(IDEVH,2,IX3) C --- Label the circle for P, T and B axis IXB(1)= IXBTP IXB(2)= IYBTP + R + 15 call xchars(whet(1),1,ixb(1),ixb(2)) c ISTAT = VGTEXT(IDEVH,IXB(1),IXB(2),WHET(1)) IXB(2)= IXB(2) - 15 IXB(3)= IXB(1) IXB(4)= IXB(2) - 5 call xmovabs(ixb(1),ixb(2)) call xdrwabs(ixb(3),ixb(4)) c ISTAT = VPLINE(IDEVH,2,IXB) IXB(2)= IYBTP - R IXB(4)= IXB(2) + 5 call xmovabs(ixb(1),ixb(2)) call xdrwabs(ixb(3),ixb(4)) c ISTAT = VPLINE(IDEVH,2,IXB) C ------ Read input data C ------------------------ cjens open(unit=3,file='test') ICOUNT=1 C ---------------------------------------------- c do 20 k = 1, 6 c c read headers c DUMMY = ' xxx ' do while (DUMMY(1:29).ne. & ' Strike Dip Rake') read(write4,'(A)',iostat=code) DUMMY call sei code(fort$,code,write4,b_flag) ! Process the outcome. enddo c 20 continue C ---------------------------------------------- 10 READ(write4,112,iostat=code) strike(icount),DIP(icount), 1 RAKE(icount) c10 READ(write4,112,iostat=code) dip(icount),strike(icount), ! lot 16.01.2002 c & RAKE(icount) call sei code( warn$, code, write4, b_flag ) if( b_flag .or. code .ne. e_ok$ ) goto 99 ANGS(1)=DIP(ICOUNT) ANGS(2)=STRIKE(ICOUNT) ANGS(3)=RAKE(ICOUNT) C ------- Calculate auxiliary planes mt=.false. CALL FMREPS(ANBTP,ANGS,PTTP,ANGS2,AN,PT,RIGHT,MT,MOMTEN, & 0,0) DIP(ICOUNT+1)=ANGS2(1) STRIKE(ICOUNT+1)=ANGS2(2) RAKE(ICOUNT+1)=ANGS2(3) C ------- Plot B,T and B axes c ------- First set graphics text alignment, .i.e. center the symbol c istat=vstaln(one,one) do 30 jj = 2, 3 if ( jj .eq. 1 ) then C --------- B axis -------------------------------------------------------- C btp(1)=anbtp(5) C btp(2)=anbtp(6) elseif ( jj .eq. 2 ) then C --------- T axis btp(1)=pttp(3) btp(2)=pttp(4) elseif ( jj .eq. 3 ) then C --------- P axis btp(1)=pttp(1) btp(2)=pttp(2) endif C -------- plunge ----------------------------------------------------- btp(2)= btp(2)*pi/180 C ------- trend az = btp(1)*pi/180.0 C ------- calculate sterographic projection of a line ax= r*dtan(pi/4.0- btp(2)/2.0) xbtp(1) = ax*dsin(az) + ixbtp xbtp(2) = ax*dcos(az) + iybtp if(cbtp(jj).eq.'P') call xset_color(color_foc_p) if(cbtp(jj).eq.'T') call xset_color(color_foc_t) call xchars(cbtp(jj),1,xbtp(1),xbtp(2)) c if(cbtp(jj).eq.'P') then write(write5,'(a1,2f8.1)',iostat=code) & cbtp(jj),xbtp(1),xbtp(2) call sei code(fort$,code,write5,b_flag) ! process the outcome. c endif c istat = vgtext(idevh,xbtp(1),xbtp(2),cbtp(jj)) 30 continue c ------ reset graphics text alignment c istat=vstaln(zero,zero) C --- icount=icount+2 GOTO 10 C --------------------------------------------------------------------- 99 CONTINUE icount=icount-1 call xset_color(color_foc_plane) C ------- Plot mechanism do 1000 j = 1 , icount delta = dip(j)*pi/180. fi = strike(j)*pi/180. ii=1 do 100 i = 1, 181 alfa=(i-1)*pi/180. -pi/2. C ------- calculate sterographic projection of plane c fi2=pi/2- acos( cos(alfa)*cos(pi/2-fi) ) C ------- coordinate transformation to the x,y,z (E,N,Z) coordinates. C ------- The xp, yp are the coordinates in the plane of the fault. C ------- xp direction of dip and yp direction of strike. xp=dcos(alfa) yp=dsin(alfa) zo=xp*dsin(delta) xo=( yp*(dsin(fi)) + xp*(dcos(delta)*dcos(fi)) ) yo=( yp*(dcos(fi)) - xp*(dcos(delta)*dsin(fi)) ) fi2 = datan(zo/dsqrt(xo*xo + yo*yo)) if ( xo .gt. 0 .and. yo .gt. 0 ) then alfa2 = datan(xo/yo) elseif ( xo .gt. 0 .and. yo .lt. 0 ) then alfa2 = datan(dabs(yo)/dabs(xo)) + pi/2. elseif ( xo .lt. 0 .and. yo .lt. 0 ) then alfa2 = datan(dabs(xo)/dabs(yo)) + pi elseif ( xo .lt. 0 .and. yo .gt. 0 ) then alfa2 = datan(dabs(yo)/dabs(xo)) + 3.*pi/2. endif ox=dtan(pi/4.-fi2/2.) xo= real(ix) + real(r)*ox*dsin(alfa2) yo= real(iy) + real(r)*ox*dcos(alfa2) x(ii) = xo x(ii+1) = yo c ------------- Plot focal mechanism if (i .eq. 1) then call xmovabs(x(ii),x(ii+1)) else call xdrwabs(x(ii),x(ii+1)) endif ii=ii+2 100 continue c c make a stroke c call xout(10.0,10.0) C ------ Plot focal mechanism c ISTAT=VPLINE(IDEVH,180,X) 1000 continue 112 FORMAT(f10.4,f10.4,f10.4) c 112 FORMAT(f8.2,f8.2,f8.2) ! lot 16.01.2002 call xset_color(color_def) RETURN END C----------------------------------------------------------------- subroutine sterpol(device,iwid,merge,rcircle,stat_plot) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Plots polarities on the focal sphere, lower hemisphere projection C C R.N.Arvidsson, 1990-05-09 Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c IMPLICIT INTEGER*2 (V,J-N) C --- include 'seidim.inc' include 'seiplot.inc' c c Libsei details... c ================= c include 'libsei.inc' ! Library definitions & data defns. external sei code ! Error condition handler. integer code ! Condition. logical b_flag ! Flag!! c integer read1, ! Read unit1. & write1, ! Write unit1. & write2, ! Ditto 2. & write3, ! & 3. & write4, ! & 4. & write5 ! & 5. common /foc_unit$ / & read1, ! Read unit1. & write1, ! Write unit1. & write2, ! Ditto 2. & write3, ! & 3. & write4, ! & 4. & write5 ! & 5. c c End of list... c ============== c real IX, IY, R, IX3(4), x(2), 1 ISTAT, IDEVH, IBC, ISTEX, ISTEY, YE(8), IN1, IN2, 2 lx(4), lc(30), iwid, one C --- CHARACTER STA(max_data)*5, INPUT*40, TEXT*3, HEADER*60, WHET(4)*1, 1 POL(max_data)*1, device*1, type*1,exp*1 integer nhead,nrecord,nphas,nstat character*80 data(max_data),txt C --- REAL*8 ax, ANGLE(max_data), AZ(max_data), PI, SP(max_data) real sangle(max_data),saz(max_data) character*5 stat(max_data) logical stat_plot C --- real rcircle c --- LOGICAL DSR, merge real xsta1,xsta2,ysta1,ysta2,ysta3,ysta4 C --- COMMON /BLCK7/ ISTAT, IDEVH, IBC, ISTEY, ISTEX, YE c common /profil/ merge, rcircle c --- DATA WHET /'N','S','E','W'/ c c Initialise... c ============= c code = e_ok$ ! Local condition. c C --- PI = 4.*atan(1.) call xset_color(color_frame) c c initialize for station code plottong around mechanism c xsta1=75.0 xsta2=930.0 ysta1=670.0 ysta2=10.0 ysta3=650.0 ysta4=10.0 C -------Draw circle if ( merge ) then IX = 500.0 IY = 350.0 else ix = 100 + 4.5*RCIRCLE IY = 200 endif IX3(1)=IX IX3(2)=IY if ( rcircle .gt. 0.001 ) then r=rcircle else R = 150 endif C --- Label the circle IX3(1)= IX IX3(2)= IY + R + 15 call xchars(whet(1),1,ix3(1),ix3(2)) IX3(2)= IX3(2) - 15 IX3(3)= IX3(1) IX3(4)= IX3(2) - 5 call xmovabs(ix3(1),ix3(4)) call xdrwabs(ix3(3),ix3(4)) IX3(2)= IY - R IX3(4)= IX3(2) + 5 call xmovabs(ix3(1),ix3(4)) call xdrwabs(ix3(3),ix3(4)) c c plot polarity legend c c txt='COMPRESSION' c call draw_circle(650.,756.,6.) c call xchars(txt,11,670.,750.) c txt='EMERGENT COMPR.' c call draw_plus(800.,756.,6.) c call xchars(txt,15,820.,750.) c txt='DILATATION' c call draw_triangle(650.,736.,6.) c call xchars(txt,10,670.,730.) c txt='EMERGENT DILAT.' c call draw_minus(800.,736.,6.) c call xchars(txt,15,820.,730.) C ------------------------ C read input event data c call indata(read1,nstat,nphas,nhead,nrecord, & type,exp,data,id) c c plot header c call xset_color(color_title) call xchars(data(1),79,0.0,760.0) c c plot solution c rewind(write4,iostat=code) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c do i=1,20 read(write4,'(a)',end=3939) txt c if(i.ge.7) call xchars(txt,79,0.0,740.0-(i-7)*20) if(i.ge.10) call xchars(txt,79,0.0,740.0-(i-7)*20) enddo 3939 continue rewind(write4,iostat=code) call sei code(fort$,code,write4,b_flag) ! Process the outcome. c c get polarities etc c CALL GETPOL(NHEAD,NRECORD,DATA,ICOUNT,POL,SANGLE,SAZ,stat) do i=1,icount angle(i)=sangle(i) az(i)=saz(i) enddo call tsend C ------- Plot mechanism do 100 i = 1, icount if ( angle(i) .gt. 90.0 ) then angle(i) = 180.0 - angle(i) az(i) = 180.0 + az(i) endif angle(i) = pi/2.0 - angle(i)*pi/180.0 az(i) = az(i)*pi/180.0 C ------- calculate sterographic projection of plane ax= r*dtan(pi/4.0-angle(i)/2.0) x(1) = ax*dsin(az(i)) + ix x(2) = ax*dcos(az(i)) + iy c c plot station code is selected c if(stat_plot) then call xset_color(color_frame) c call xchars(stat(i),4,x(1)-30.0,x(2)-20.0) call xmovabs(x(1),x(2)) c c select where line go to for station code c if(x(1).lt.500.0.and.x(2).gt.350.0.and.ysta1.gt.0.0) * then call xdrwabs(xsta1,ysta1+7.0) call xchars(stat(i),5,xsta1-40.0,ysta1) call xchars(pol(i),1,xsta1-70.0,ysta1) ysta1=ysta1-20.0 endif if(x(1).lt.500.0.and.x(2).le.350.0.and.ysta1.lt.750.0) * then call xdrwabs(xsta1,ysta2+7.0) call xchars(stat(i),5,xsta1-40.0,ysta2) call xchars(pol(i),1,xsta1-70.0,ysta2) ysta2=ysta2+20.0 endif if(x(1).ge.500.0.and.x(2).gt.350.0.and.ysta3.gt.0.0) * then call xdrwabs(xsta2,ysta3+7.0) call xchars(stat(i),5,xsta2,ysta3) call xchars(pol(i),1,xsta2+65.0,ysta3) ysta3=ysta3-20.0 endif if(x(1).ge.500.0.and.x(2).le.350.0.and.ysta4.lt.750.0) * then call xdrwabs(xsta2,ysta4+7.0) call xchars(stat(i),5,xsta2,ysta4) call xchars(pol(i),1,xsta2+65.0,ysta4) ysta4=ysta4+20.0 endif endif if ( pol(i) .eq. '+' ) then call xset_color(color_foc_comp) ! lot 21.01.2002 call draw_plus(x(1),x(2),6.0) c call draw_circle(x(1),x(2),3.5) elseif ( pol(i) .eq. 'C' ) then call xset_color(color_foc_comp) call draw_circle(x(1),x(2),6.0) elseif ( pol(i) .eq. '-' ) then call xset_color(color_foc_comp) call draw_minus(x(1),x(2),6.) ! lot 21.02.2002 elseif ( pol(i) .eq. 'D' ) then call xset_color(color_foc_dilat) call draw_triangle(x(1),x(2),6.) ! lot 16.01.2002 elseif ( pol(i) .eq. 'V' ) then call xchars('V',1,x(1),x(2)) elseif ( pol(i) .eq. 'S' ) then call xchars('S',1,x(1),x(2)) elseif ( pol(i) .eq. 'H' ) then call xchars('H',1,x(1),x(2)) endif 100 continue 1000 continue call xset_color(color_def) C --- 111 FORMAT(A) 112 FORMAT(a4,f4.0,f8.0,4x,a1,f7.0) RETURN END c------------------------------------------------------------------ subroutine getpol(nhead,nrecord,data,npol,pol,ain,az,stat) c c get polarity data, jh, sep 92 c c input nhead,nrecord,data, see indata routine c output npol : number of polarity data c ain : angle of incidence c az : azimuths c stat : stations c implicit none integer nhead,nrecord character*80 data(*) character*5 stat(*) character*1 pol(*),sense integer npol real ain(*),az(*) c--reading same as integers integer iain,iaz c--counter integer i c c npol=0 do i=nhead+1,nrecord-1 c if((data(i)(17:17).eq.'D'.or.data(i)(17:17).eq.'C'.or. c * data(i)(17:17).eq.'-'.or.data(i)(17:17).eq.'+'). c * and.data(i)(11:11).eq.'P') then call get_sense(data(i),sense,0) if (sense.ne.' ') then if(data(i)(58:60).ne.' '.and.data(i)(77:79).ne.' ') * then npol=npol+1 c pol(npol)=data(i)(17:17) pol(npol)=sense ! lot 21.01.2002 read(data(i)(58:60),'(i3)') iain ain(npol)=iain read(data(i)(77:79),'(i3)') iaz az(npol)=iaz stat(npol)=data(i)(2:6) endif endif enddo c return end c--------------------------------------------------------- subroutine get_sense(text,sense,option) c c return sense as defined in original version of focmec c Lars Ottemoller, 21.01.2002 c implicit none character*80 text ! phase line character*1 sense ! output sense integer option ! 0: use E or I (used for plotting) ! 1: ignore E or I (used for computation) sense=' ' c compression if (text(17:17).eq.'C'.or.text(17:17).eq.'U') then if (text(11:11).eq.'P') then if (text(10:10).eq.'I'.or. & text(10:10).eq.' '.or. & option.eq.1) then sense='C' elseif (text(10:10).eq.'E') then sense='+' endif endif ! P c dilatation elseif (text(17:17).eq.'D') then if (text(11:11).eq.'P') then if (text(10:10).eq.'I'.or. & text(10:10).eq.' '.or. & option.eq.1) then sense='D' elseif (text(10:10).eq.'E') then sense='-' endif endif ! P elseif (text(17:17).eq.'V') then sense='V' elseif (text(17:17).eq.'S') then sense='S' elseif (text(17:17).eq.'H') then sense='H' endif return end c----------------------------------------------------------------- SUBROUTINE FCIRCL(R,X,Y) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INPUT R - RADIUS C X - CENTER OF CIRCLE C Y - CENTER OF CIRCLE C PURPOSE PLOTS A FILLED CIRCLE C C R.N. Arvidsson 1990-05-04 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c implicit integer*2 (v,j-n) real r, x, y, pi real ix(722), istat, idevh, ibc, istey, istex, ye(8) COMMON /BLCK7/ ISTAT, IDEVH, IBC, ISTEY, ISTEX, YE pi = 4.d0*asin(1.d0/(sqrt(2.d0))) ii=1 do 100 i = 1, 361 c ix(ii) = x + r*dcos(i*pi/180.) ix(ii) = x + r*cos(i*pi/180.) c ix(ii+1)= y + r*dsin(i*pi/180.) ix(ii+1)= y + r*sin(i*pi/180.) c ---------- Plot circe if (ii .eq. 1) then call xmovabs(ix(ii),ix(ii+1)) else call xdrwabs(ix(ii),ix(ii+1)) endif ii=ii+2 100 continue return end subroutine draw_circle(x,y,r) c c draws and open circle c call xmovabs(x,y) call xdrwabs(x,y+r) call xmovabs(x+r,y) k=12 z=6.28/k do i=1,k x1=r*cos(i*z)+x y1=r*sin(i*z)+y call xdrwabs(x1,y1) enddo return end subroutine closed_circle(x,y,r) ir=r call xmovabs(x,y) call xdrwabs(x,y) do i=2,ir call draw_circle(x,y,float(i)) enddo return end subroutine draw_triangle(x,y,r) c c draw triangle inside circle of radius r c lot 16.01.2002 c implicit none real x,y,r call xmovabs(x,y) call xdrwabs(x,y+r) call xdrwabs(x-r*sin(3.14/3.),y-r*cos(3.14/3.)) call xdrwabs(x+r*sin(3.14/3.),y-r*cos(3.14/3.)) call xdrwabs(x,y+r) return end subroutine draw_square(x,y,r) c c draw triangle inside circle of radius r c lot 16.01.2002 c implicit none real x,y,r,a a=r*cos(45./180.*3.14) call xmovabs(x-a,y+a) call xdrwabs(x-a,y-a) call xdrwabs(x+a,y-a) call xdrwabs(x+a,y+a) call xdrwabs(x-a,y+a) return end subroutine draw_plus(x,y,r) c c draw plus inside circle of radius r c lot 16.01.2002 c implicit none real x,y,r,a c the plus call xmovabs(x,y+r) call xdrwabs(x,y-r) call xmovabs(x-r,y) call xdrwabs(x+r,y) c square around it c call draw_square(x,y,r/cos(45./180.*3.14)) return end subroutine draw_minus(x,y,r) c c draw minus inside circle of radius r c lot 16.01.2002 c implicit none real x,y,r,a call xmovabs(x-r,y) call xdrwabs(x+r,y) c square around it c call draw_square(x,y,r/cos(45./180.*3.14)) return end subroutine foc_prepare(npol,nratio,vpvs) c c routine reads the hyp print.out and hyp.out files and make an c input file to focmec called focmec.inp having the same content c as the hyp.out except that angles of incidence are now c written in column under signal to noise ratio. If more than c one event is present, it is assumed that all should be used in a c composite fault plane solution and only header lines for the first c event is used. C C C J. HAVSKOV, october 1992 C C C LATEST UPDATE: c nov 93 lates bl routine c sep 94 multiple events for composite fault plane solutions c COMPUTERTYPE IMPLICIT NONE c c Libsei details... c ================= c include 'libsei.inc' ! Library definitions & data defns. include 'seidim.inc' external sei code, ! Error condition handler. & sei open, ! Open file handler. & sei close ! Close file handler. integer code ! Condition. logical b_flag ! Flag!! c integer read1, ! Read unit1 & read2, ! Ditto 2 & write1, ! Write unit1 & write2, ! Write unit2 & write3 ! Write unit3 c c End of list... c ============== c c--data array for one event CHARACTER*80 DATA(max_data) character*80 text C NUMBER OF DIFFERENT STATIONS IN DATA INTEGER NSTAT C NUMBER OF HEADERS AND RECORDS IN DATA INTEGER NHEAD,NRECORD C ANGLE OF INCIDENCE AND AZIMUTH REAL ANGINC,AZ integer iang C NUMBER OF DIFFERENT PHASE FOR EVENT INTEGER NPHASE INTEGER I,id,k,iprint,c,ind,l integer nevent ! number of events c NUMBER OF POLARITIES AND NUMBER OF AMPLITUDE RATIOS integer*2 npol,nratio c number of amplitude readings integer namp c amplitude real amp c station codes for amp readings character*5 amp_stat(max_data) c amplitude related parameters real amp_amp(max_data) ! amplitudes real amp_per(max_data) ! period real amp_dist(max_data) ! distance real amp_trtime(max_data) ! travel times real amp_az(max_data) ! azimuth real amp_anginc(max_data) ! angle of incidence real amp_angemg(max_data) ! angle of emergence real amp_qcor(max_data) ! q-correction character*1 amp_comp(max_data) ! component character*2 amp_phase(max_data) ! phase type character*1 ratio_type ! V, H or S c real qzero_p,qzero_s,qalpha_p,qalpha_s ! q-parameters real free_correction ! free surface correction real s_tstar,p_tstar ! tstar integer ibot,itop ! index for ratio character*1 type,expl c-- sense - one-character key for polarity, C=cpompression,D=Dilatation,+=emergent C,- = eme.dil. character*1 sense c-- one line in print file character*80 oneline(max_data+100) c velocity model integer max_layer parameter (max_layer=200) real vdepth(max_layer),vp(max_layer),vs(max_layer) integer imoho,nlayer real vpvs c distance and depth real dist,depth c p and s velocities from iasp91 real velp,vels c angle of emergance at station real aemrgp,aemrgs c phase name character*3 phasenm c free surface correction real correction c amp ratio real ratio c ratio name character*11 ratio_name c flag logical new real xx c c Initialise... c ============= c code = e_ok$ ! Local condition. nevent=0 nlayer=0 npol=0 nratio=0 namp=0 amp_stat(1)=' ' c c get q-paramters c call get_focmec_def(qzero_p,qzero_s, *qalpha_p,qalpha_s,p_tstar,s_tstar) c c open print file c call sei open( old$, ! Open (stop on error). & ' ', ! No prompt. & 'print.out', ! This file. & read1, ! On this unit. & b_flag, ! Exists?!! (n/a). & code ) ! Returned condition. c c open hyp.out file c call sei open( old$, ! Open (stop on error). & ' ', ! No prompt. & 'hyp.out', ! This file. & read2, ! On this unit. & b_flag, ! Exists?!! (n/a). & code ) ! Returned condition. c c open output file c call sei open( unknown$, ! Open (stop on error). & ' ', ! No prompt. & 'focmec.inp', ! This file. & write1, ! On this unit. & b_flag, ! Exists?!! (n/a). & code ) ! Returned condition. call sei open( unknown$, ! Open (stop on error). & ' ', ! No prompt. & 'pspolar.inp', ! This file. & write2, ! On this unit. & b_flag, ! Exists?!! (n/a). & code ) ! Returned condition. call sei open( unknown$, ! Open (stop on error). & ' ', ! No prompt. & 'focmec.dat', ! This file. & write3, ! On this unit. & b_flag, ! Exists?!! (n/a). & code ) ! Returned condition. c c here starts loop to read hyp.out and print.out c 1 continue call INDATA(read2,nstat,nphase,NHEAD, & NRECORD,TYPE,EXPL,DATA,id) c c check for end of file c if(nrecord.eq.0) goto 99 nevent=nevent+1 c c read forward in print file c 22 continue ! back here to read next line read(read1,'(a80)',iostat=code) oneline(1) c write(6,*) oneline(1) call sei code(fort$,code,read1,b_flag) ! Process the outcome. if( b_flag ) then 3423 continue ! Premature end of file. write(6,*) * ' No stations in file or no location,', * ' only old solutions can be ploted' write(6,*)' or' write(6,*)' Focmec operated outside EEV and hyp.out and' write(6,*)' print.out do not correspond, relocate with HYP' write(6,*)' and run focmec again' write(6,*)' Continue (y/n) ? ' read(5,'(a)')text if(text(1:1).eq.'y'.or.text(1:1).eq.'Y') then data(1)(72:72)='.' ! indicate only fp write(write1,'(a80)',iostat=code) (data(i),i=1,nhead) goto 3 endif call sei close(delete$,write1,code) ! Delete output file. call sei close(close$,read2, code) ! Close inputdata. call sei close(close$,read1, code) ! Ditto print.out. goto 9999 ! Return to caller. end if ! c c read vp/vs ratio, lot jan 2002 c if(oneline(1)(23:29).eq.'Vp/Vs =') then read(oneline(1)(31:35),'(f5.2)') vpvs endif c c read model, lot jan 2002 c if(oneline(1)(1:32).eq.' Depth, km Vp, km/s Vs, km/s'.and. & nlayer.eq.0) then do while (oneline(1).ne.' ') nlayer=nlayer+1 read(read1,'(a80)',end=3423) oneline(1) read(oneline(1)(1:32),'(f10.2,1x,f10.2,1x,f10.2)') & vdepth(nlayer),vp(nlayer),vs(nlayer) if(oneline(1)(36:36).eq.'N')imoho=nlayer enddo nlayer=nlayer-1 c do i=1,nlayer c write(*,*) i,vdepth(i),vp(i),vs(i) c enddo endif if(oneline(1)(2:4).eq.'stn') go to 3 go to 22 ! cntinue to read print.out until station is found c c station line found c 3 continue c c read print file for one event c iprint=1 4 read(read1,'(a80)',iostat=code) oneline(iprint) call sei code(fort$,code,read1,b_flag )! Process the outcome. if( b_flag ) goto 99 ! End of file. c c check for end of phases for one event c if(oneline(iprint)(1:13).eq.' ') go to 80 iprint=iprint+1 goto 4 c c event read c 80 iprint=iprint-1 c c read depth c read(data(1)(39:43),'(f5.1)') depth c c some info written out c if(nevent.eq.1) then ! write headers for first event write(write1,'(a80)',iostat=code) (data(i),i=1,nhead) write(write3,'(a80)',iostat=code) data(1)(2:79) ! lot 16.01.2002 endif call sei code(fort$,code,write1,b_flag) ! Process outcome. c c------------------------------------------------------------------------------ c loop for writing out polarities in focmec.dat, input for focmec_exe c------------------------------------------------------------------------------ c do i=nhead+1,nrecord-1 c c check that stations without station coordinates are not used, was possible c before dec 2008, jh c if(data(i)(58:79).ne.' ') then call get_sense(data(i),sense,1) ! fix polarity notation if (sense.ne.' ') then write(write1,'(a80)',iostat=code) data(i) ! focmec.inp call sei code(fort$,code,write1,b_flag) ! Process outcome. c-------------------------------------------------------- c polarity data for original focmec, in file focmec.dat c-------------------------------------------------------- c c get input for focmec_exe, lot 16.01.2002 c read(data(i)(11:12),'(a2)') phasenm read(data(i)(77:79),'(f3.0)') az read(data(i)(58:60),'(f3.0)') anginc read(data(i)(71:75),'(f5.0)') dist c c note, only 4 letter station codes can be used, does not matter here, only info c phase name written out for information c write(write3,'(a4,2f8.2,a1,18x,a)') * data(i)(2:5),az,anginc,sense,data(i)(10:15) npol=npol+1 ! count polarity readings c c c write to file pspolar.inp: station, epicentral distance, angle of incidence, c polarity, used for gmt plot c write(write2,'(a5,1x,a3,1x,a3,1x,a1)') & data(i)(2:6),data(i)(77:79), & data(i)(58:60),data(i)(17:17) endif endif enddo c c-------------------------------------------------------------------------------- c-------------------------------------------------------------------------------- c loop for reading and writing out amplitude ratios c--------------------------------------------------------------------------------- c---------------------------------------------------------------------------------- c---------- c check all lines for amplitude data for focmec c---------- namp=0 do i=nhead+1,nrecord-1 if ((data(i)(11:14).eq.'AMPG'.or.data(i)(11:14).eq.'AMPN'.or. * data(i)(11:14).eq.'AMSG'.or.data(i)(11:14).eq.'AMSN').and. * nevent.eq.1.and. * (data(i)(8:8).eq.'Z'.or.data(i)(8:8).eq.'R'. ! make sure correct component * or.data(i)(8:8).eq.'T')) then namp=namp+1 c c read station and component c amp_stat(namp)=data(i)(2:6) amp_comp(namp)=data(i)(8:8) c c amplitude type: , P or S for first letter, G or N for second letter. c G is direct phase (only local) and N is first arrival (local c or global c amp_phase(namp)=data(i)(13:14) c c read amplitude and period c read(data(i)(34:45),'(g7.1,f5.1)') amp_amp(namp), * amp_per(namp) c c check that not zero c if(amp_amp(namp).eq.0.0.or.amp_per(namp).eq.0.0) then write(6,*)amp_stat(namp), * ' zero amplitude or period, skipped' namp=namp-1 goto 1200 endif c c read corresponding azimuth and epicentral distance c read(data(i)(77:79),'(f3.0)') amp_az(namp) read(data(i)(71:75),'(f5.0)') amp_dist(namp) c c--------------------- c angle of incidence c--------------------- c c find angle of incidence at source and angle of emergence at station. c for first arrival, find angle of incidence in data array. c For local G phase calculate c c PROBLEM: get_ang does not take station elevation into account so angle claculated c might be 1-2 deg off if station has signicant elevation and distance c is short. To be fixed !!! c c c G phase for local model c if(amp_phase(namp)(2:2).eq.'G'.and.(data(1)(22:22).eq.'L'. * or.data(1)(22:22).eq.'R')) * then call get_ang_emrg(depth,amp_dist(namp),nlayer,vp,vs, * vdepth,imoho,'PG',amp_angemg(namp),amp_anginc(namp)) goto 1111 ! angle found, no more checking endif c c N-phase or first arriving phase, local and global. c first find angle of incidence c in data array, phase can be any refracted phase and not nescessarely c moho refracted. phase can also be PG if just given as P. If phase is c specifically given as PN, it is assumed it is the first arrival although c PN2, PN3 etc could have been first arrival. It is important here to c check for the user spcified phase, not program identified phase. remeber, c angle of incidence in data array is for program identified phase c assume same angle for P and S, so calculate for P. c do k=nhead+1,nrecord if(amp_stat(namp).eq.data(k)(2:5).and. ! station * (data(k)(11:12).eq.'PN'.or.data(k)(11:12).eq.'P ')) * then read(data(k)(58:60),'(i3)') l amp_anginc(namp)=l goto 1111 endif enddo c c if here, phase not found, skip this amplitude c write(6,*) amp_stat(namp), ' phase not found for ain' namp=namp-1 go to 1200 1111 continue c c------------------------------- c calculate angle of emergence for local model, assume first arriving P if c not PG c------------------------------- c if(data(1)(22:22).eq.'L'.or.data(1)(22:22).eq.'R') then if(amp_phase(namp)(2:2).eq.'G') then call get_ang_emrg(depth,amp_dist(namp),nlayer,vp,vs, * vdepth,imoho,'PG',amp_angemg(namp),amp_anginc(namp)) else call get_ang_emrg(depth,amp_dist(namp),nlayer,vp,vs, * vdepth,imoho,'P ',amp_angemg(namp),amp_anginc(namp)) endif else c c calculate angle of emergence for global case using ray parameter and the c iasp91 model c call iasp_vel(depth,velp,vels) ! get velocity at hypocentral depth amp_angemg(namp)=sin(amp_anginc(namp)/57.3)* * ((6371.0-depth)*5.8)/(velp*6371.0) ! assume top vp=5.8 amp_angemg(namp)=57.3*asin(amp_angemg(namp)) endif c c------------------------------- c get travel times c------------------------------- c c for Q-correction, just used first arrivng phase c since very little difference between Pg and Pn. In any case, c most observations use first arrivng phases c do k=1,iprint c c check if station and phase, use user identified phase c if(oneline(k)(2:6).eq.amp_stat(namp)) then ! same station if(amp_phase(namp)(1:1).eq.'P') then ! case of P if(oneline(k)(26:27).eq.'P '.or.oneline(k)(26:27). * eq.'PG'.or.oneline(k)(26:27).eq.'PN') then read(oneline(k)(58:64),'(f7.1)') * amp_trtime(namp) goto 1112 endif endif if(amp_phase(namp)(1:1).eq.'S') then !case of S if(oneline(k)(26:27).eq.'S '.or.oneline(k)(26:27). * eq.'SG'.or.oneline(k)(26:27).eq.'SN') then read(oneline(k)(58:64),'(f7.1)') * amp_trtime(namp) goto 1112 endif endif endif enddo c c if here, phase for travel time not found c c c if P travel time there, use that to calculate s-travel time c do k=1,namp if(amp_phase(namp)(1:1).eq.'S'.and.amp_stat(k). * eq.amp_stat(namp).and.amp_phase(k)(1:1).eq.'P') then amp_trtime(namp)=amp_trtime(k)*1.78 write(6,*)amp_stat(namp), * ' use P-travel time to calculate S-travel time' goto 1112 endif enddo c c no travel time found c write(6,'(1x,a,1x,a,1x,a)') amp_stat(namp),amp_phase(namp), * 'Phase not found for travel time, read phase' namp=namp-1 goto 1200 ! skip this observation c 1112 continue c c-------------------- c correction for Q c-------------------- c if(data(1)(22:22).eq.'L'.or.data(1)(22:22).eq.'R') then if(amp_phase(namp)(1:1).eq.'P') amp_qcor(namp)= * exp((3.14*amp_trtime(namp))/ * (amp_per(namp)*qzero_p*(1.0/amp_per(namp)**qalpha_p))) if(amp_phase(namp)(1:1).eq.'S') amp_qcor(namp)= * exp((3.14*amp_trtime(namp))/ * (amp_per(namp)*qzero_s*(1.0/amp_per(namp)**qalpha_s))) else if(amp_phase(namp)(1:1).eq.'P') * amp_qcor(namp)=exp(3.14*p_tstar/amp_per(namp)) if(amp_phase(namp)(1:1).eq.'S') * amp_qcor(namp)=exp(3.14*s_tstar/amp_per(namp)) endif endif ! endif of amplitude line 1200 continue ! jump to here if some information not found for amp enddo ! enndo for searching data array c c write q-values used c write(6,*) write(6,'(a,f6.1,a,f4.2,a,f6.1,a,f4.1,a,a,f4.2,a,f4.2)') *' Q: Local: Qp=',qzero_p,'**',qalpha_p, *' Qs=',qzero_s,'**',qalpha_s,' Global:', *' t*(P)=',p_tstar,' t*(S)=',s_tstar write(6,*) c c write out observations c write(6,'(a,a)') *' STAT C PH AMP PER TRTIME QCOR ANGINC ANGEMG', *' AZ DIST' do i=1,namp write(6,'(1x,a5,1x,a1,1x,a2,1x,i9,2x,f5.2,1x, * f6.1,1x,f6.1,2i7,1x,i4,1x,i5)') * amp_stat(i),amp_comp(i), * amp_phase(i),int(amp_amp(i)), * amp_per(i),amp_trtime(i),amp_qcor(i),int(amp_anginc(i)), * int(amp_angemg(i)),int(amp_az(i)),int(amp_dist(i)) enddo c------------------------------------------------------ c calculate amplitude ratios c------------------------------------------------------ c write(6,*) write(6,'(a)') * ' STAT Ratio type T Amp 1 Amp 2 Fcor LogRat' do i=1,namp do k=i+1,namp c c if same phase type and station, make ratio c if(amp_stat(i).eq.amp_stat(k). * and.amp_phase(i)(2:2).eq.amp_phase(k)(2:2)) then ratio_type=' ' ratio_name=' ' itop=0 c write(6,*) amp_stat(i),' ',amp_comp(i),' ',amp_comp(k), c * ' ',amp_phase(i),' ',amp_phase(k) c c case of dividing with P amplitude, but not on T c if((amp_phase(k)(1:1).eq.'P'.and.amp_comp(k).eq.'T') * .or.(amp_phase(i)(1:1).eq.'P'.and.amp_comp(i).eq.'T')) * then write(6,*) * amp_stat(i),' Do not use P on T, reading skipped' goto 777 endif c if(amp_phase(k)(1:1).eq.'P'.and. * amp_phase(i)(1:1).eq.'S') then itop=i ibot=k endif if(amp_phase(i)(1:1).eq.'P'.and. * amp_phase(k)(1:1).eq.'S') then itop=k ibot=i endif c c case of SV divided by SH c if(amp_phase(k)(1:1).eq.'S'.and. * (amp_comp(k).eq.'Z'.or.amp_comp(k).eq.'R').and. * amp_phase(i)(1:1).eq.'S'.and. * amp_comp(i).eq.'T') then itop=k ibot=i endif if(amp_phase(i)(1:1).eq.'S'.and. * (amp_comp(i).eq.'Z'.or.amp_comp(i).eq.'R').and. * amp_phase(k)(1:1).eq.'S'.and. * amp_comp(k).eq.'T') then itop=i ibot=k endif if(itop.gt.0) then ratio=(amp_amp(itop)*amp_qcor(itop))/ * (amp_amp(ibot)*amp_qcor(ibot)) else c c ratio S on Z to S R should not be used, so not an error c if(amp_phase(i)(1:1).eq.'S' * .and.amp_phase(k)(1:1).eq.'S' * .and.((amp_comp(i).eq.'Z'.and.amp_comp(k).eq.'R') * .or. (amp_comp(i).eq.'R'.and.amp_comp(k).eq.'Z'))) * goto 777 c c invalid names c write(6,*) ' Something wrong with names' write(6,*) amp_stat(i),' ',amp_comp(i),' ', * amp_comp(k),' ',amp_phase(i),' ',amp_phase(k) goto 777 endif c c ratio name c c c ratio with P(Z) at denominator c if(amp_phase(ibot)(1:1).eq.'P'.and.amp_comp(ibot). * eq.'Z') then if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'Z') then ratio_name = 'SV(Z)/P(Z)' ratio_type='V' endif if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'T') then ratio_name = 'SH(T)/P(Z)' ratio_type='H' endif if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'R') then ratio_name = 'SV(R)/P(Z)' ratio_type='V' endif endif c c ratio with P(R) at denominator c if(amp_phase(ibot)(1:1).eq.'P'.and.amp_comp(ibot). * eq.'R') then if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'Z') then ratio_name = 'SV(Z)/P(R)' ratio_type='V' endif if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'T') then ratio_name = 'SH(T)/P(R)' ratio_type='H' endif if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'R') then ratio_name = 'SV(R)/P(R)' ratio_type='V' endif endif c c ratio with SH at denominator c if(amp_phase(ibot)(1:1).eq.'S'.and.amp_comp(ibot). * eq.'T') then if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'R') then ratio_name = 'SV(R)/SH(T)' ratio_type='S' endif if(amp_phase(itop)(1:1).eq.'S'.and.amp_comp(itop). * eq.'Z') then ratio_name = 'SV(Z)/SH(T)' ratio_type='S' endif endif c c make free surface correction c free_correction=0.0 call freesurf_ratio(vpvs,ratio_name, * amp_angemg(i),amp_angemg(i),free_correction) if(free_correction.eq.0.0) then write(6,*) ' Unable to do free surface correction' write(6,*) ' Use 1.0' free_correction=1.0 endif c c correct amplitude ratio and take log c ratio=alog10(ratio*free_correction) c c write to screen c write(6,'(1x,a5,1x a11,1x,a1,1x,2i9,1x,f5.1,1x,f6.2)') * amp_stat(i),ratio_name, ratio_type,int(amp_amp(itop)), * int(amp_amp(ibot)),free_correction,ratio c c write to focmec input file, focmec.dat c nratio=nratio+1 write(write3,'(a4,2f8.2,a1,f8.4,1x,a1,1x,f6.2,1x,a)') * amp_stat(i)(1:4),amp_az(i),amp_anginc(i),ratio_type, * ratio,' ',amp_anginc(i),ratio_name c c write to focmec.inp c write(write1,'(1x,a5,10x,a1,40x,i3,16x,i3)') * amp_stat(i),ratio_type,int(amp_anginc(i)),int(amp_az(i)) endif 777 continue ! here if soemthing wrong with component name or skip enddo enddo c c back for next event c goto 1 c c End of file... c ============== c ...end of all events c 99 write(write1,*,iostat=code) & ' ' ! last line is blank call sei code(fort$,code,write1,b_flag) ! Process outcome. call sei close(close$,read1,code) ! Close print.out. call sei close(close$,read2,code) ! & data input. call sei close(close$,write1,code) ! & output. call sei close(close$,write2,code) ! & output. call sei close(close$,write3,code) ! & output. goto 9999 ! Return to caller. c c Return to Caller... c =================== c 9999 return end c subroutine angles c c interactively calculate different fault plane solution angles c using focmec subroutines c real PTTP(4) ! strike and dip of P and T real ANGS(3) ! dip, strike and rake of principal fault plane real ANGS2(3) ! dip, strike and rake of auilary ------------ real ANBTP(6) real momten(6) real pi,degrad integer i character*1 choice c PI = 4.0*ATAN(1.0) degrad=180.0/pi c 50 continue write(6,*) write(6,*) ' Relation between fault plane solution angles:' write(6,*) ' 1: Strike, dip and rake from P and T' write(6,*) ' 2: Other plane and P and T from first plane' write(6,*) ' Return to stop' read(5,'(a1)') choice if(choice.eq.'1') goto 1 if(choice.eq.'2') goto 2 if(choice.eq.' ') return c 1 continue write(6,*)' Enter strike and dip of P and strike and dip of T' read(5,*)pttp c aa=.true. c bb=.false. do i=1,4 pttp(i)=pttp(i)/degrad enddo call pttpin(pttp,angs,angs2,anbtp,momten,pi) do i=1,3 angs(i)=angs(i)*degrad angs2(i)=angs2(i)*degrad enddo write(6,'(a,2x,3f8.1)') *' Strike, dip and rake of principle fault plane' *,angs(2),angs(1),angs(3) write(6,'(a,2x,3f8.1)') *' Strike, dip and rake of auxilary fault plane' *,angs2(2),angs2(1),angs2(3) c call c *FMREPS(ANBTP,ANGS,PTTP,ANGS2,bb,aa,bb,-1,-1) c write(6,*)angs c return goto 50 c c 2 continue c c calculate from strike dip and rake c write(6,*)' Enter strike, dip and rake' read(5,*) angs(2),angs(1),angs(3) do i=1,3 angs(i)=angs(i)/degrad enddo call DSRIN (ANGS,ANBTP,ANGS2,PTTP,momten,PI) do i=1,3 angs2(i)=angs2(i)*degrad enddo do i=1,4 pttp(i)=pttp(i)*degrad enddo write(6,'(a,2x,3f8.1)') *' Strike, dip and rake of auxiliary fault plane' *,angs2(2),angs2(1),angs2(3) write(6,*)' Strike and dip of P and strike and dip of T' write(6,*) pttp goto 50 return end subroutine get_ang_emrg(evdepth,dist,mnlayer,mvp,mvs,mdepth,imoho, & phasenm,aemrg,ann) c c Lars Ottemoller, Jan 2002 c determine emergence angle at station for given velocity model, c distance and event depth c implicit none include 'hypparm.inc' integer seiclen real dist ! distance in km real evdepth ! hyp depth in km real mvp(*),mvs(*),mdepth(*) ! layer velocity model integer mnlayer ! number of layers integer imoho,iconrad ! index for moho real aemrg ! angle of emergance for p and s real rearth ! Earth radius character*(*) phasenm ! phase name, PN,PG,P,SN,SG,S real pi real radtodeg integer i character*8 phsid character*1 prmd,prm2 real tmin,delta,ann integer iflag,unit,minflag,iustat if (dist.eq.0.) then aemrg=0. return endif unit=99 pi=asin(1.)*2. radtodeg=45./atan(1.) rearth=6371. iconrad=0 c set test parameter defaults call settest(test) do i=1,mnlayer v(i)=0. if (phasenm(1:1).eq.'S') v(i)=mvs(i) if (phasenm(1:1).eq.'P') v(i)=mvp(i) parm(i)=v(i) d(i)=mdepth(i) if (v(i).eq.0.) return enddo do i=1,mnlayer-1 parm(mnlayer+i)=mdepth(i+1)-mdepth(i) enddo nl=mnlayer ! used in dtdx2 prmd=' ' c c if not a first arrival, specify arrival, G or N. first arrival could c e.g. PN3 or P c if (seiclen(phasenm).eq.2) prmd=phasenm(2:2) c write(6,*) 'phase', prmd,imoho xh(1)=dist/rearth xh(2)=0. xh(3)=evdepth x0(1)=0. x0(2)=0. x0(3)=0. call dtdx2(xh,x0,prmd,imoho,iconrad,unit,tmin, & dx,delta,ann,iflag,phsid) c c compute emergence angle from horizontal slowness c aemrg = radtodeg * asin(dx(1)/rearth*v(1)) c write(*,*) ' aemrg = ',aemrg,' ain= ',ann, c & ' tmin = ',tmin return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine freesurf_ratio(vpvs,ratio_name,aemrgp,aemrgs,ratio) implicit none c c compute free surface correction based on emergance angle, based c on program freesurf by A. Snoke c c Lars Ottemoller, 31.01.2002 c c input: vpvs - ratio vp/vs c ratio_name - ratio to be computed, see below c aemrgp,aemrgs - emergence angle for p and s c c output: ratio - amplitude correction factor c c hane c dec 08 jh change SG on T to SH, add missing factor for SH/V c real radtodeg ! factor for rad to deg real vpvs,vsvp ! ratios vp/vs and vs/vp real aemrgp,aemrgs ! emergance angle at station real apv,asv,apr,asr ! vertical and radial displacement real ash ! real phr ! for sv inc beyond critical, phv=phr+90. real ratio ! amplitude correction character*(*)ratio_name! name of ratio to be returned in ratio radtodeg=45./atan(1.) vsvp=1./vpvs ash = 2. c get p amplitudes call gdmot(radtodeg,1,vsvp,aemrgp,apr,phr,apv) c get s amplitudes call gdmot(radtodeg,2,vsvp,aemrgs,asr,phr,asv) asr = abs(asr) asv = abs(asv) c write(*,*) 'P amplitudes on Z and R:',apv, apr c write(*,*) 'S amplitudes on Z and R:',asv, asr if (ratio_name.eq.'SV(R)/P(Z)') then ratio = apv/asr elseif (ratio_name.eq.'SV(R)/SH(T)') then ratio = ash/asr elseif (ratio_name.eq.'SV(R)/P(R)') then ratio = apr/asr c elseif (ratio_name.eq.'SV(Z)/P(Z)') then ratio = apv/asv elseif (ratio_name.eq.'SV(Z)/P(R)') then ratio = apr/asv elseif (ratio_name.eq.'SV(Z)/SH(T)') then ratio = ash/asv c elseif (ratio_name.eq.'SH(T)/P(Z)') then ratio = apv/ash elseif (ratio_name.eq.'SH(T)/P(R)') then ratio = apr/ash endif return end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine iasp_vel(depth,velp,vels) c c finds iap91 velocity velp , vels, at depth depth c only hardwired down to 760 km c implicit none character*1 phase real depth real velp, vels real x ! normalized radius x=(6371.0-depth)/6371.0 if(depth.lt.20.0) then velp=5.8 vels=3.36 endif if(depth.ge.20.0.and.depth.lt.35.0) then velp=6.5 vels=3.75 endif c if(depth.ge.35.0.and.depth.lt.120.0) then velp=8.78541-x*0.74953 vels=6.706231-2.248585*x endif c if(depth.ge.120.0.and.depth.lt.210.0) then velp=25.41389-x*17.69722 vels=5.75020-1.27420*x endif c if(depth.ge.210.0.and.depth.lt.410.0) then velp=30.78765-x*23.25415 vels=15.24213-11.08552*x endif c if(depth.ge.410.0.and.depth.lt.660.0) then velp=29.38896-x*21.40656 vels=17.70732-13.50652*x endif c if(depth.ge.660.0.and.depth.lt.760.0) then velp=25.96984-x*16.93412 vels=20.76890-16.53147*x endif c if(depth.gt.760) then write(6,*) *' Too deep to calculate iasp velocities, will use vp=8 and vs=3.5' velp=8 vels=3.5 endif c return end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine get_focmec_def(qzero_p,qzero_s,qalpha_p,qalpha_s, *p_tstar,s_tstar) c C get focmec defaults, so far only for attenuation c c written by J. Havskov c c implicit none external sei get file, ! Search directories & open file. & sei close, ! Close file open. & sei code ! Error processor. include 'libsei.inc' ! Library definitions & data defns integer code ! Condition. logical b_eof ! End of file?. integer def_unit ! unit to read from c real qzero_p,qzero_s,qalpha_p,qalpha_s ! q-parameters real free_correction ! free surface correction real s_tstar,p_tstar ! tstar character*80 data(100) ! content of parameter file integer i,nline c c set default parameters in case there is no focmec.def c p_tstar=1.0 s_tstar=4.0 qzero_p=100.0 qalpha_p=1.0 qzero_s=100.0 qalpha_s=1.0 c c open and read default file with stations and components to use c -------------------------------------------------------------- c call sei get file( open$+ignore$, ! Find and open without messages. & def_unit, ! On file unit. & code, ! Condition (n/a). & 'DAT', ! Alternative directory to search. & 'FOCMEC.DEF' ) ! For this file. c c read file if there... c ------------ c if(code.ne.e_ok$) goto 666 i=1 10 read(def_unit,'(a)',iostat=code) data(i) ! Read from file. call sei code( fort$, ! Process fortran i/o condition. & code, ! Condition. & def_unit, ! On unit. & b_eof ) ! End of file?. c if( b_eof ) then ! End of file. nline=i ! Store records in file. call sei close( close$, def_unit, code ) ! Close (Default stop on error). c else ! Otherwise. i=i+1 ! Increment record number. goto 10 ! Read another record. end if ! c c read parameters c do i=1,nline c c local Q c if(data(i)(1:22).eq.'QZERO and QALPHA FOR P'.and. * data(i)(41:50).ne.' ') then read(data(i),'(40x,2f10.1)',err=99) qzero_p,qalpha_p endif c if(data(i)(1:22).eq.'QZERO and QALPHA FOR S'.and. * data(i)(41:50).ne.' ') then read(data(i),'(40x,2f10.1)',err=99) qzero_s,qalpha_s endif c c tstar if(data(i)(1:17).eq.'TSTAR for P AND S'.and. * data(i)(41:50).ne.' ') then read(data(i),'(40x,2f10.1)',err=99) p_tstar,s_tstar endif enddo goto 9999 c c no def file c 666 continue write(6,*) 'No FOCMEC.DEF file, use defaults' goto 9999 c c errror c 99 continue write(6,*)' Errror in MULPLT.DEF, use defaults' c 9999 continue return end c