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 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. c c synthetic c c A synthetic input can be genrated for FOCMEC. This is done automartically c for a fps in S-file which has method SYNTET written in the F-file. Synthetic c polarities and amplitudes will then be generated for the same observations c as given in the S-file. These observations are normally written out to the file c focmec.dat as polarites and corrected amplitudes. The focmec.dat is then input c to the original focmec. The synthetic observations are now made by correcting the c focmec.dat file to the theortical polarities and amplitude ratios corresponding c to the SYNTET fault plane solutions. The theoretical values are therefore at the c source and not affected by Q and free surface. Since the theoretical values are c generated by focmec routines, focmec should exactly return the SYNTET solution. c Before makin the solution with focmec, it is possible to edit the focmec.dat file c to introduce erorrs or remove observations in order to test the sensitivity of c the solution errros in input. c vp/vs is hardwired to 1.74 c The same input is also used with HASH so HASH can also be tested in the same way. c c important files: c c focmec.inp: a nordic file with one or all events if composoite c solution, use to store all polarities and amplitudes c or plotting, also used to stores old fps and newly c selected solutions. the solution stored in s-file, c if run from eev, is taken from here. generated by foc_prepare c focmec.dat: Input file for the original focmec program, made by c routine foc_prepare from hyp.out and print.out. The c observations are now corrected for Q and free surface. c fps.out : Each time a solutiopn is selected (saved) the solution c is written to fps.out, solutions are accumulated. 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 aug 27 10 jh : skip questions, use option 1 if prompt input is o, plot different c solution in different colors, reformat graphics screen text c output c sep 21 jh : also different colors P and T c oct 1 10 jh : error in amp ratio read wrongly c oct 20 jh : fix foc-prepare so amplitudes can be used in composite solution c nov 4 jh : add argument p so only foc prepare is run c fix so no crash in finding minimum number of polarities when c amplitudes are present c dec 22 10 jh : gfortran pc: remove winplot include, modify read of print.out, c implicit none, unit check, computer type check, remove c call to tsend c jan 11 11 jh : fix so more than 9 amplitude ratio errors possible, fix c 80 chars written to to focmec.inp, was 79, pc gfortran problem c feb 10 11 jh : check overflow in q-corr in write out c feb 22 11 jh : size of plot from color.def c mar 6 11 jh : correct output format,was old c mar 17 11 jh : fix so fps can be plotted without readings, option o c apr 1 11 jh : change status of scratch file to unknow, made problem with permisisons c on a linux system c may 31 12 jh : add mt plot, look for mt plot c sep 19 13 jh : change error message when no phase with amplitude, c accpet Pg,Pn,Sg and Sn for finding travel time c (suggestion by Christian Baillard) c 2015.06.04 pv : small change due to compiler warning on GOTO c 2016 01 29 jh : implement automatic amplitudes c 2016 07 19 jh : implement synthetic test c 2016 07 27 jh : number of bad amp fit and amp error was not written to c output F-line, set default polarities 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 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. c INTEGER*2 NERRP,NERRA ! allowed errors for pol and amp ratio INTEGER 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*90 text ! gfortran pc character*80 data(max_data) character chr_solution *(80) ! Selected solution. integer deldeg ! search increment logical foc_plot_only real cutp, cuts ! hardwired for limiting calc. amp at nodal planes c -- tektronix options integer nars ! number of arguments character*80 arg(5) ! arguments C--- CHARACTER answer*1 character*80 editor real aa,bb,cc ! help var C --- LOGICAL merge real xx(4) ! store strike,dip and rake, number of bad polarities logical sun,pc,linux C -------------- c include 'version.inc' foc_plot_only=.false. 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 args... 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 or p for just plot 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.'o'.and.arg(1)(1:1).ne.'p') then write(6,*)' Invalid argument, use:' write(6,*)' a: Calculate angles' write(6,*)' o: Plot saved solution, no questions' write(6,*)' p: Just run foc preparation' write(6,*)' None: Fault plane solution after using HYP' stop endif endif if(arg(1)(1:1).eq.'o') foc_plot_only=.true. c--- limits for amplitudes at or near nodal planes cutp=0.05 cuts=0.15 call foc_prepare(npol,nratio,vpvs,foc_plot_only,editor,cutp,cuts) c c stop if just prepare c if(arg(1)(1:1).eq.'p') stop c c both hard copy file and screen plot c plotoption=1 c--- hard copy type is postscript hctype=1 c-- set window sixe wsize=80 call get_window_size if(size_focmec.gt.0) wsize=size_focmec ! from color.def c 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 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( unknown$, ! 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( unknown$, ! 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,'(a)') 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 c c write(6,'(a)') chr_solution ifoc=0 chr_solution(79:80)=' F' chr_solution(71:78)='FOCMEC ' c c put number of bad pol in correct position c chr_solution(61:62)=chr_solution(35:36) chr_solution(35:36)=' ' c c amplitude error fit c chr_solution(57:60)=chr_solution(63:66) chr_solution(63:66)=' ' c c c number of ratio errors c chr_solution(64:65)=chr_solution(55:56) chr_solution(55:56)=' ' c c delete what is not transferred c chr_solution(38:49)=' ' c write(6,'(a)') chr_solution do i=2,nhead c write(6,'(a)') data(i) if(data(i)(79:80).eq.' F') then ifoc=1 data(i)=chr_solution endif enddo c write(6,*) ifoc 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) c c also write to fps.out c read(chr_solution,'(3f10.1)') aa,bb,cc call add_fps(aa,bb,cc,'FOCMEC ','F') 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 c c bypass if all defaults c 150 continue if(arg(1)(1:1).eq.'o') then iselect=1 stat_plot=.false. else c150 write(6,*) 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 endif ! end of sel. loop 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' if(arg(1).eq.'o') then write(6,*) 'Enter to continue' read(5,'(a)') text stop ! stop in defualt mode endif 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 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. c write(6,*)'nerra',nerra 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. if(text.eq.' ') then max_amp_error=0.2 else read(text,*,err=136) max_amp_error endif 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, enter for default 2' read(5,'(a)') text if(text.eq.' ') then deldeg=2 else read(text,*) deldeg endif 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' c c make sure formatting is right, funny read in focmec_exe c if(nerra.lt.10) write(write6,'(i1,7x,a)') * nerra,'allowed amplitude errors' if(nerra.gt.9.and.nerra.lt.100) write(write6,'(i2,6x,a)') * nerra,'allowed amplitude errors' if(nerra.gt.99) write(write6,'(i3,5x,a)') * nerra,'allowed amplitude errors' c write(write6,'(i1,7x,a)') nerra,'allowed amplitude errors' write(write6,'(f4.2,4x,a)') cutp, & 'lower limit P cutoff' write(write6,'(f4.2,4x,a)') cuts, & '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------------------------------------------------------------------------- c find solutions with minimum polarity errors c------------------------------------------------------------------------- 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 9988 continue 191 continue rewind(write6,iostat=code) c c check if any solutions c if(minerrp.eq.99999) then write(6,*) write(6,*)' No solution, increase error limit for amplitudes' goto 100 endif 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 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) c 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 c for mt curves c integer nmt1,nmt2 ! number of points in two mt graphs real xmt(10000),ymt(10000) ! mt data sets, both in same vectors 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 c c plot p and t c c if(cbtp(jj).eq.'P') call xset_color(color_foc_p) c if(cbtp(jj).eq.'T') call xset_color(color_foc_t) ii=(icount+1)/2 if(ii.lt.5) then call xset_color(ii) else call xset_color(6) ! if more than 4 solution, use black endif c call xset_color(icount/2) 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) ! no effect now jh aug 2010 c c mt plot c c c read mt curves c nmt1=0 nmt2=0 open(11,file='beachball',status='old',err=445) read(11,*,end=445) nmt1 do i=1, nmt1 read(11,*,err=445)xmt(i),ymt(i) enddo read(11,*,end=445)nmt2 if(nmt2.gt.0) then do i=nmt1+1, nmt2+nmt1 read(11,*,err=445)xmt(i),ymt(i) enddo endif goto 446 445 continue nmt1=0 446 continue close(11,status='delete') c c plot mt c if(nmt1.gt.0) then c c scale up mt to r c do i=1,nmt1+nmt2 xmt(i)=xmt(i)*r+ix ymt(i)=ymt(i)*r+iy enddo c c plot mt curves c do i=1,nmt1 if (i.eq.1) then call xmovabs(xmt(i),ymt(i)) else call xdrwabs(xmt(i),ymt(i)) endif enddo do i=nmt2+1,nmt1+nmt2 if (i.eq.nmt1+1) then call xmovabs(xmt(i),ymt(i)) else call xdrwabs(xmt(i),ymt(i)) endif enddo endif c C ------- Plot mechanisms, that is fault lines c do 1000 j = 1 , icount ! icount is twice the number of mechanisms ii=(j+1)/2 if(ii.lt.5) then call xset_color(ii) else call xset_color(6) ! if more than 4 solution, use black endif 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 c 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,txt1 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 real strike,dip,rake 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 from s-file 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 c write the different solutions, but only the first ones (10) c c find first line c do i=1,20 read(write4,'(a)',end=3939)txt if(txt(5:10).eq.'Strike') goto 2222 enddo goto 3939 2222 continue txt=' STR DIP RAK Source ' call xchars(txt,21,0.0,740.0) do i=1,10 read(write4,'(3f10.1,40x,a8)',end=3939) * strike,dip,rake,txt1(1:8) i1=strike+0.5 i2=dip+0.5 i3=rake+0.5 txt=' ' write(txt(1:12),'(3i4)') i1,i2,i3 txt(14:21)=txt1(1:8) txt(22:80)=' ' c if(i.ge.7) call xchars(txt,79,0.0,740.0-(i-7)*20) c if(i.ge.10) call xchars(txt,79,0.0,740.0-(i-7)*20) j=i if(i.ge.5) j=6 ! black call xset_color(j) call xchars(txt,21,0.0,740.0-(i)*20) enddo call xset_color(6) ! back to normal 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) c write(6,*) 'number of polarities',icount c write(6,*) (data(i),i=1,nrecord) c c c find if run with synthetics, then take polarities from focmec.dat c to get synthetic polarities c do i=1,nhead if(data(i)(80:80).eq.'F'.and.data(i)(71:76).eq.'SYNTET') * then write(6,*) * 'Polarities plotted are from synthetic solution', * ' with polarities and angles given in focmec.dat' open(112,file='focmec.dat',status='old') read(112,'(a)') pol(i) ! just skip first line do j=1,icount read(112,'(20x,a1)') pol(j) enddo close(112) goto 3451 endif enddo 3451 continue do i=1,icount angle(i)=sangle(i) az(i)=saz(i) enddo 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 if 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 c call draw_plus(x(1),x(2),6.0) ! c jh call draw_circle(x(1),x(2),6.0) 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_dilat) c call draw_minus(x(1),x(2),6.) ! lot 21.02.2002 call draw_triangle(x(1),x(2),6.) ! jh oct 2010 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 c write(6,*)'pi' 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,foc_plot_only,editor,cutp,cuts) c c routine reads the hyp print.out and hyp.out files and make an c input file to focmec called focmec.dat. 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.foc_plot_only indicate only plot. c c the print.out file is used to read the model used for free surface coreection c and the travel times used for q-corrections c C J. HAVSKOV, october 1992 C C C LATEST UPDATE: c nov 93 latest bl routine c sep 94 multiple events for composite fault plane solutions c 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!! integer seiclen 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 character*80 editor 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,j 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 integer amp_manual ! number of manual amps integer amp_auto ! -----------automatic integer amp_auto_spec ! -------------------- spectral amps integer amp_type ! one of the above real strike,dip,rake ! for syntetic test real aasin logical foc_plot_only character*1 type,expl, answer real cutp,cuts ! for calcualting ratios 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*90 oneline(max_data+100) ! gfortran pc 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 amp fre surface amplifications real ar,av c angle real phr 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 write(6,*) 'indata',nrecord !cx 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,'(a)',iostat=code) oneline(1) ! gfortran pc c write(6,*) oneline(1) ! cx call sei code(fort$,code,read1,b_flag) ! Process the outcome. if( .NOT.b_flag ) goto 3424 c if( b_flag ) then 3423 continue ! Premature end of file. if(foc_plot_only) then write(write3,'(a80)') data(1)(2:79) write(write1,'(a80)',iostat=code) (data(i),i=1,nhead) goto 99 endif 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 does not correspond' write(6,*)' Make sure all events get located' write(6,*)' relocate 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. c end if 3424 continue c c read model, but only for first event c if(nevent.eq.1) then c ! 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,'(a)',end=3423) oneline(1) ! gfortran pc 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 !cx c write(*,*) i,vdepth(i),vp(i),vs(i) c enddo endif endif ! end reading model if(oneline(1)(2:4).eq.'stn') go to 3 go to 22 ! continue 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,'(a)',iostat=code) oneline(iprint) ! gfortran pc 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 write(6,*)' iprint',iprint 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 skip amplitude stuff if only plot of solution c if(foc_plot_only)goto 99 c c-------------------------------------------------------------------------------- c-------------------------------------------------------------------------------- c loop for reading and writing out amplitude ratios c--------------------------------------------------------------------------------- c---------------------------------------------------------------------------------- c c check what kind of amplitude data c amp_manual=0 amp_auto=0 amp_auto_spec=0 do i=nhead+1,nrecord-1 if(data(i)(11:14).eq.'AMSG'.or.data(i)(11:14).eq.'AMPN') * amp_manual=amp_manual+1 if(data(i)(11:14).eq.'ATPG'.or.data(i)(11:14).eq.'ATSG') * amp_auto=amp_auto+1 if(data(i)(11:14).eq.'ASPG'.or.data(i)(11:14).eq.'ASSG') * amp_auto_spec=amp_auto_spec+1 enddo 500 continue write(6,*) write(6,'(a,i5)') 'Number of polarities: ', npol write(6,'(a,3(a,i5,2x))')'Amplitude types: ', *' Manual: ',amp_manual, ' Automatic:',amp_auto, *' Spectral: ',amp_auto_spec c c skip question if no amplitudes c amp_type=1 if(amp_manual.eq.0.and.amp_auto.eq.0.and.amp_auto_spec.eq.0) *goto 502 if(npol.eq.0.and.amp_manual.eq.0.and.amp_auto.eq.0. *and.amp_auto_spec.eq.0) *then write(6,*)' No data ************** will stop' stop endif write(6,*) write(6,'(a,a)') *'Amplitude to use: Manual(1), Automatic(2), Spectral(3) ?' read(5,*) amp_type if(amp_type.ne.1.and.amp_type.ne.2.and.amp_type.ne.3) then write(6,*)'Wrong type' goto 500 endif 502 continue write(6,*) 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'.or. * data(i)(11:14).eq.'ATPG'.or.data(i)(11:14).eq.'ASPG'.or. * data(i)(11:14).eq.'ATSG'.or.data(i)(11:14).eq.'ASSG').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 c c check for chosen type of amplitude c if(amp_type.eq.1.and.data(i)(11:12).ne.'AM') goto 1200 if(amp_type.eq.2.and.data(i)(11:12).ne.'AT') goto 1200 if(amp_type.eq.3.and.data(i)(11:12).ne.'AS') goto 1200 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 phase name could be lower case, so correct c if(amp_phase(namp)(2:2).eq.'g') amp_phase(namp)(2:2)='G' if(amp_phase(namp)(2:2).eq.'n') amp_phase(namp)(2:2)='N' 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. remember, 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*aasin(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, use print file c to get travel times 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'.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'.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,'(a,1x,a,a,1x,a,1x,a)') 'For station ', * amp_stat(namp),' , ',amp_phase(namp), * 'Phase not found for travel time, read a P or S phase' namp=namp-1 goto 1200 ! skip this observation c 1112 continue c c-------------------- c correction for Q c-------------------- 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)) if(amp_per(namp).lt.1.0) then write(6,*) * ' *** Warning *** Period too low for distant event' endif 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', *' Fcor AZ DIST' do i=1,namp c c calculate free surface correction, just for write out, correction is done c later but could in fact have been done here c j=2 ! s phase if(amp_phase(i)(1:1).eq.'P') j=1 ! p phase call gdmot(57.3,j,1.0/vpvs,amp_angemg(i),ar,phr,av) if(j.eq.2) av=ar ! correction for free surface if(amp_phase(i)(1:1).eq.'S'.and.amp_comp(i).eq.'T') av=2.0 ! SH c c check q-corr c if(amp_qcor(i).gt.10000.) amp_qcor(i)=9999.0 write(6,'(1x,a5,1x,a1,1x,a2,1x,i9,2x,f5.2,1x, * f6.1,f7.1,2i7,f5.1,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)),av,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. gfortran: one blank was missing at end of line c write(write1,'(1x,a5,10x,a1,40x,i3,16x,i3,a1)') * 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 continue cxx if(foc_plot_only) return ! for plotting only, do not look for synthetic c c check if synthetic should be calculated, input from s-file c do i=1,nhead if(data(i)(80:80).eq.'F'.and.data(i)(71:76).eq.'SYNTET') then read(data(i)(1:30),'(3f10.1)') strike, dip,rake write(6,*) write(6,*)' This run with synthetic data' write(6,*)' Strike, dip and rake are: ',strike, dip,rake write(6,*) * ' Free surface and Q corrections are not affecting result' write(6,*) call make_fps_amp(strike,dip,rake,cutp,cuts) c c option to edit focmec.dat c write(*,*) ' Do you want to edit focmec.dat (y/n=enter) ? ' answer=' ' read(5,'(a1)') answer if (answer.eq.'y'.or.answer.eq.'Y' ) then text=editor(1:seiclen(editor))// ' focmec.dat' write(*,*) text(1:seiclen(text)) call systemc(text,seiclen(text)) endif c c check if number of observations have been changed by editing c open(112,file='focmec.dat',status='old') k=0 ! polarities l=0 ! amplitudes 4848 continue read(112,'(a)',end=4849) text if(text(21:21).eq.'D'.or.text(21:21).eq.'C') k=k+1 if(text(21:21).eq.'V'.or.text(21:21).eq.'S' * .or.text(21:21).eq.'H') l=l+1 goto 4848 4849 continue npol=k nratio=l close(112) return endif enddo 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 aasin 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 * aasin(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 FOCMEC.DEF, use defaults' c 9999 continue return end c function aasin(x) real x if(x.lt.-1.0) then write(6,*)' arcsin overflow', x x=-1.0 endif if(x.gt.1.0) then write(6,*)' arcsin overflow', x x=1.0 endif aasin=asin(x) return end subroutine make_fps_amp(strike,dip,rake,cutp,cuts) c c subroutine to generate synthetic input amplitude and polarity to c test fps programs focmec and hash in seisan c real xyz(9) ! help variables real lograt(3) ! log amplitude ratio for sv/p,sh/p, xx/xx real bot(3) ! p-amplitude with sign character*3 flag character*80 text(1000) ! temporary values of focmec.dat real strike,dip,rake RD = 45.0/ATAN(1.0) open(112,file='focmec.dat',status='old', err=10) goto 20 10 continue write(6,*)'error' write(6,*)'no focmec.dat file' stop 20 continue c c read to end c i=1 30 continue read(112,'(a)',end=40) text(i) i=i+1 goto 30 40 continue n=i-1 c write(6,*) n c write(6,*) strike,dip,rake c c write(6,'(a)')' STAT C AZIM AIN RATIO NOD' c c calculate the synthetics c do i=2,n c c get next data c read(text(i)(7:20),'(f6.2,f8.2)')az,ain c c calcualate using part of subroutine from focmec c azin=az ptoang=ain stoang=ain TREND = AZIN/RD PLUNGE = (90.0 - PTOANG)/RD COST = COS(TREND) SINT = SIN(TREND) COSP = COS(PLUNGE) SINP = SIN(PLUNGE) XYZ(1) = COST*COSP XYZ(2) = SINT*COSP XYZ(3) = SINP SPLUNG = (90.0 - STOANG)/RD SINP = SIN(SPLUNG) COSP = COS(SPLUNG) C C Next two vectors reversed in sign from normal convention because C of my (snoke) convention for SV and SH (down and left, facing the station) C XYZ(4) = -COST*SINP XYZ(5) = -SINT*SINP XYZ(6) = +COSP XYZ(7) = SINT XYZ(8) = -COST XYZ(9) = 0.0 vpvs3=5.26 dp1=dip/rd st1=strike/rd ra1=rake/rd c write(6,*) c write(6,*) strike,dip,rake,az,ain call LRATIO(1,DP1,ST1,RA1,XYZ,VPVS3,LOGRAT(1), 1 TOP,bot(1),CUTP,CUTS,FLAG) call LRATIO(2,DP1,ST1,RA1,XYZ,VPVS3,LOGRAT(2),TOP,bot(2), 1 CUTP,CUTS,FLAG) call LRATIO(3,DP1,ST1,RA1,XYZ,VPVS3,LOGRAT(3),TOP,bot(3), 1 CUTP,CUTS,FLAG) c c fish out what is needed c c c polarity c if(text(i)(21:21).eq.'C'.or.text(i)(21:21).eq.'D') then if(bot(2).lt.0.0) then text(i)(21:21)='D' else text(i)(21:21)='C' endif write (6,'(1x,a4,1x,a1,1x,2f7.1,7x,2x,a3)') * text(i)(1:4),text(i)(21:21),az,ain,flag endif c c amplitude c c SV(Z)/P(Z) if(text(i)(21:21).eq.'V') then write(text(i)(22:29),'(f8.4)') lograt(1) write (6,'(1x,a4,1x,a1,1x,2f7.1,1x,f7.4,1x,a3)') * text(i)(1:4),text(i)(21:21),az,ain,lograt(1),flag endif c SH(T)/P(Z) if(text(i)(21:21).eq.'H') then write(text(i)(22:29),'(f8.4)') lograt(2) write (6,'(1x,a4,1x,a1,1x,2f7.1,1x,f7.4,1x,a3)') * text(i)(1:4),text(i)(21:21),az,ain,lograt(2),flag endif c SV(Z)/SH(T) if(text(i)(21:21).eq.'S') then write(text(i)(22:29),'(f8.4)') lograt(3) write (6,'(1x,a4,1x,a1,1x,2f7.1,1x,f7.4,1x,a3)') * text(i)(1:4),text(i)(21:21),az,ain,lograt(3),flag endif enddo rewind(112) do i=1,n write(112,'(a80)') text(i) enddo close(112) return end