program autosig c c input files c c 1 - one waveform file c 2 - s-file and 1 waveform file c 3 - s-file, multiple waveform files c c c input data c c 1 - no location, no phase readings -> find phase arrival c 2 - location, no phase readings -> compute phase arrivals, or force 1 c 3 - phase readings -> use phase arrivals, or force 1/2 c c c input_mode c c 1 - input file is Nordic file, name given through environmental variable c 2 - S-file given through environmental variable, used by seisnet c 3 - input is Nordic file c 4 - continuous data c c c functions c c 1 - pick phase based on sta/lta c 2 - check if spectral level of signal is above noise at its maximum c 3 - determine type (L,R,D,N) based on spectral levels at different c frequencies c 4 - determine spectral flat level and corner frequency based on c spectrum of P or S, if distance known to correct for geom spr and Q, c using genetic algorithm c ( 5 - read amplitude for Ml and Mb ) c ( 6 - location ) c c c by Lars Ottemoeller, Dec 1999 c c changes: c 08.02.2000 lo: use file_ready c 09.02.2000 lo: if no detrigger use last sample for coda c 10.02.2000 lo: put event id 'R' c 23.05.2000 lo: use new wav reading routines c 23.05.2000 lo: compute spectrum, new routine auto_moment c 14.02.2001 lo: change in sta_lta_nrec, at BGS c 15.02.2001 lo: add call to get_seisan_def, at BGS c 21.02.2001 lo: fix call to spectrum c taken out mulplt.inc c 05.03.2001 lo: add routine to determine type of event from spectrum c 08.03.2001 lo: fixes in sta_lta_nrec and find_coda c 14.03.2001 lo: add genetic algortihm to determine om and f0 from spectrum c ----------------------------------------- c 20.03.2001 lo: name of program changed to autosig c 21.03.2001 lo: changed variable odata to raw_signal c 23.03.2001 lo: add correction for near surface attenuation using kappa c 30.03.2001 lo: add get_amp_per, auto amplitude routine c 20.06.2001 lo: added function to use const Q for long periods c 03.07.2001 lo: removed call to dtime c 23.09.2002 lo: add all phase output c 11.11.2002 lo: fix bug with no output for waveform file input c 11.11.2002 lo: add arguments c 13.11.2002 lo: bug fix on input_mode c 16.06.2003 lo: add reading of continuous data, not yet finished c 24.05.2005 lo: put in recent changes c 07.2008 lo: some changes, compute spec average c 22.02.2009 lo: search for delta kappa as well c change to finding low frequency c 12.03.2010 jh: add poles and zeros to subroutine get_amp_per, c currently not used, also added to get_om_f0, c freq_event_type c 22 03 2013 jh: add wav_mem_init, remeove msdif, proc_time and other unnused pars in c call to get_om_f0, comment out run of expect script c 22.05.2013 lo: added plotunit, but no plotting yet c 23 12 2015 jh: add test to call update_spec, only test(116) is c used and set locally since normally read form c STATION0.HYP c c on Linux routines auto_moment and invert_spectrum cause problems c implicit none C Seisan library inserts and routines... C ====================================== C include 'libsei.inc' ! Seisan definitions include 'seidim.inc' ! ------------------ include 'seisan.inc' ! ------------------ include 'waveform.inc' ! ------------------ include 'autosig.inc' ! autosig include file include 'rea.inc' ! rea include file external computer_type ! Get platform type. external sei open, ! Open file routine. & sei close, ! Close file routine. & sei get file ! Find files in diff dirs integer sei clen ! Length of string integer r_index real test(200) ! reset test form hyp character*1 dchar integer write01,write02,err_unit,! File units, & read01,read02,code ! Local condition. logical b_flag ! Flag end of file c -------------------------------------------------------------- character*80 infile,outfile,autosig_out, ! file names & outf,xfile,autosig_err integer*4 idata(max_sample) ! integer data real raw_signal(max_sample),fil_signal(max_sample) ! original and filtered data integer it ! trace counter character*80 fheader(max_trace) ! file headers character*1040 theader ! trace header integer c,i,x,z,j,y,ind ! counters integer coda(max_trace) ! coda duration double precision drate ! rate as double prec. character*80 line,text1,text2 ! some text line character*80 sfile ! name of sfile character*80 question ! text character*1 cbyte(max_trace) ! 2 or 4 byte character*5 net_code ! network code character*29 mainhead_text ! mainheader text double precision phase_time(max_phase) ! phase times double precision msec,msec2 character*80 aphases(max_phase) ! phase array real sta,lta,ratio,mincoda,dtrle,fill,filh ! variables from parameter file ! for individual stations integer in ! 1 if filenr.lis integer filec ! file counter real cof(8) ! filter coefficients real gain ! filter gain integer npasses ! number of passes, 1 forward, ! 2 both ways real fsn1,fsn2 ! spectral band with good s/n character*80 data(max_data) ! event arrays integer nstat ! no of stations for event c integer nphase ! no of phases for one event c integer nrecord ! no of records for event c integer nhead ! no of headers for event character*1 exp,evid ! event ids-------------- integer id ! event id line number character*2 fr_type ! event type based on spectrum integer dist_cnt,local_cnt,reg_cnt! counter for type of event real om,f0,omi,f0i ! flat level and corner frequency real ka ! computed delta kappa integer nspec ! number of spectra c real spfmin,spfmax ! frequnecy range in which signal above noise logical clipped ! true if signal clipped logical lfreeze ! true if to freeze lta at trigger integer nfiles ! number of files real duration ! signal duration integer phase_index ! index to phase in rea structure integer spec_index ! index to spec phase integer jday ! julian day logical hyp_flag ! true if hypocenter known real slat,slon,sheight,sdist ! station coords and distance real degtorad ! constant to convert degree to radians real distdeg,az0,xx(4),yy double precision abs_time character*1 phase character*1 which_phase real geo_corr integer trace_unit,event_unit, ! file units & mom_unit,f0_unit,mag_unit,res_unit,plotunit real res_sum real proc_count,time real amp,per real mag_sum,mom_sum,mag_save(100) character*80 data_save(max_data) integer nrecord_save logical event_done integer mag_n real res_m_a(2) double precision time_amp real pi real max_dist,min_dist integer mag_ml_index parameter (pi=3.141592654) character*12 time1 ! time strings character*14 time2 integer narg ! number of arguments character*80 arg(10) ! arguments logical onlyonce ! true if input file given on prompt or ! from env variable real cmin,rsec integer year,month,day,doy,hour,min,sec character*14 end_time logical flag character*80 file_name character*3 agency c c print version c include 'version.inc' out_version_date='July 23, 2001' if (version_new) out_version_date=version_date call print_ver test(116)=0.0 ! average mags, set 1.0 for median c c get arguments c call get_arguments(narg,arg) call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'autosig.trace', ! Filename. & trace_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. write(trace_unit,*) ' AUTOSIG TRACE FILE ' call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'autosig.event', ! Filename. & event_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'moment.comp', ! Filename. & mom_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'f0.comp', ! Filename. & f0_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'mag.comp', ! Filename. & mag_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'res.comp', ! Filename. & res_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'autosig.eps', ! Filename. & plotunit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. c c get parameters from parameter file c call get_autosig_par(trace_unit) call get_seisan_def call dir_char(dchar) call get_agency(agency) c c init variables c infile = ' ' outfile='autosig.wav.001' autosig_out ='autosig.out' autosig_err ='autosig.err' filec = 0 npasses = 1 clipped=.FALSE. input_mode = 1 ! default waveform data degtorad=0.0174533 proc_count=0. res_sum=0. max_dist=-1. min_dist=99999. in=0 onlyonce=.false. cwav=.false. c c overwrite settings if set as arguments c do i=1,narg if (arg(i)(1:5).eq.'-spec') then if (arg(i+1)(1:2).eq.'on') then lautospec=.true. else lautospec=.false. endif elseif (arg(i)(1:6).eq.'-clear') then if (arg(i+1)(1:2).eq.'on') then lclear=.true. else lclear=.false. endif elseif (arg(i)(1:6).eq.'-phase') then if (arg(i+1)(1:2).eq.'on') then lautophase=.true. else lautophase=.false. endif elseif (arg(i)(1:7).eq.'-infile') then infile=arg(i+1)(1:seiclen(arg(i+1))) elseif (arg(i)(1:5).eq.'-help') then write(*,*) 'autosig [-infile -spec on/off '// & '-phase on/off -clear on/off' endif enddo c c open error output file c call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & autosig_err, ! Filename. & err_unit, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. line=' ' c c check if started from other program c call get_message('autosig',line) 40 continue c c in case autosig is started from other program c if (line(1:7).eq.'autosig'.and.line(9:9).ne.' ') then call get_env_event(sfile) c c read sfile c c call file_ready(sfile) call sei open( old$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & sfile, ! Filename. & write02, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. call rea_event_in(write02,.true.,data,code) if (code.ne.0) then write(*,*) ' error reading file ' call sei close( close$, write02, code) stop endif call sei close( close$, write02, code) read(line(9:79),'(a73)') infile input_mode = 2 onlyonce=.true. lwavout= .FALSE. elseif (infile.ne.' '.and.in.ne.1) then onlyonce=.true. else c c if autosig not started from other program read input file name c question=' Filename or number, filenr.lis for all, or cont' call filename(question,infile) endif if (infile(1:3).eq.'EOF') STOP if(infile(1:10).eq.'filenr.lis') then c call file_ready('filenr.lis') call sei open( old$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & 'filenr.lis', ! Filename. & read02, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. in=1 elseif (infile(1:4).eq.'cont') then if (n_cont_base.eq.0) then write(*,*) ' no continuous database defined in SEISAN.DEF ' stop endif c c write out the database names c write(*,*)'Databases are:' do i=1,n_cont_base write(*,*)' ',cont_base(i) enddo write(*,*) ' Start time (yyyymmddhhmmss) ' read(5,'(a)') cwav_start_time write(*,*) ' End time (yyyymmddhhmmss) ' read(5,'(a)') end_time write(*,*) ' Time interval in minutes ' read(5,*) cmin c c set interval length c cont_interval=cmin*60. c c set end time for continuous data c read(cwav_start_time,'(i4,4i2,f2.0)') year,month,day, & hour,min,rsec call timsec(year,month,day,hour,min,rsec,msec) msec=msec+cont_interval call sectim(msec,year,doy,month,day,hour,min,rsec) write(cwav_end_time,'(i4.0,5i2.0)') & year,month,day,hour,min,int(rsec) lwavout=.false. lautoamp=.false. lautospec=.false. lall_phases=.true. onlyonce=.true. input_mode=4 cwav=.true. call cwav_time_limits(0) endif c c open ascii output c if (input_mode.ne.2) then call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & autosig_out, ! Filename. & write02, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. elseif (input_mode.eq.2) then call sei open( old$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & sfile, ! Filename. & write02, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. endif 50 continue c c if filenr.lis c if(in.eq.1) then read(read02,'(7x,a)',end=40) infile if(infile(1:4).eq.' ') stop filec=filec+1 write(outfile(13:15),'(i3)') filec do i=13,15 if (outfile(i:i).eq.' ') outfile(i:i)='0' enddo endif c c open input file c if (input_mode.ne.4) then call wav_init call wav_mem_init wav_filename(1)=' ' wav_filename(1)=infile c c now read first header for one event c call get_full_wav_name(infile, wav_filename(1)) write(*,*) ' INPUT FILE : '// & wav_filename(1)(1:seiclen(wav_filename(1))) if (lwavout) then write(*,*) ' OUTPUT FILE: '//outfile(1:seiclen(outfile)) endif call read_wav_header(1) else wav_error_message=' ' endif 60 continue if (cwav) then rea_nwav=0 n_cont_trace=0 endif if (input_mode.eq.3) wav_error_message='MODE 3' write(*,*) ' wav_error_message= '//wav_error_message if(wav_error_message.ne.' ') then c c file could be nordic c if (input_mode.ne.3) then c c open when checking 1st time c call sei open( old$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & infile, ! Filename. & read01, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. endif nspec=0 mag_n=0 mag_sum=0. mom_sum=0. mag_ml_index=0 call rea_event_in(read01,.true.,data,code) write(*,*) data(1) c c save data vector c do i=1,rea_nrecord data_save(i)=data(i) enddo nrecord_save=rea_nrecord if ((rea_nrecord.eq.0.or.code.ne.0).and.input_mode.eq.3) then call sei close( close$, read01, code) goto 9999 endif onlyonce=.true. if (rea_nrecord.eq.0.and.in.ne.1) then write(6,'(1x,a)') wav_error_message goto 9999 else input_mode=3 write(*,*) ' Nordic input file ... ' c c read hyp parameters c hyp_flag=.false. if (hyp_lat(1).ne.0..or. & hyp_lon(1).ne.0..or. & hyp_depth(1).ne.0.) hyp_flag=.true. call timsec(hyp_year(1),hyp_month(1),hyp_day(1), & hyp_hour(1),hyp_min(1),hyp_sec(1),abs_time) c write(*,*) ' lat/lon/dep ',hyp_lat(1),hyp_lon(1),hyp_depth(1) j=0 c c find index to ml value c do i=1,6 if (hyp_mag_type(i,1).eq.'L') then mag_ml_index=i endif enddo c c get waveform filenames c call wav_init nfiles=0 wav_nfiles=0 do i=1,rea_nwav xfile=' ' xfile=rea_wav(i)(2:79) xfile=xfile(1:seiclen(xfile)) call get_full_wav_name(xfile,wav_filename(i)) enddo c c read headers c do i=1,rea_nwav call read_wav_header(i) if(wav_error_message.ne.' ') then write(6,'(1x,a)') wav_error_message endif enddo endif if (input_mode.eq.1) then if (in.eq.1) goto 50 goto 9999 endif c c file is waveform file, set content for output file c elseif (input_mode.eq.1) then input_mode=1 hyp_flag=.false. c hyp_fix_org(1)=' ' c hyp_depth_flag(1)=' ' call rea_hyp_clear(1) hyp_lat(1)=0. hyp_lon(1)=0. hyp_depth(1)=0. hyp_year(1)=wav_year(1) hyp_month(1)=wav_month(1) hyp_day(1)=wav_day(1) hyp_hour(1)=wav_hour(1) hyp_min(1)=wav_min(1) hyp_sec(1)=wav_sec(1) rea_nphase=0 rea_nhyp=1 rea_nwav=1 rea_wav(1)=' '// & ' ' write(rea_wav(1)(2:),'(a)') infile(1:seiclen(infile)) rea_wav(1)(80:80)='6' call systime(time1,time2) rea_id_line=' ACTION:ASI '//time2//' OP: '// & 'STATUS: ID:20021001025954 I' cACTION:REE 02-10-01 10:02 OP:gdf STATUS: ID:20021001025954 write(rea_id_line(61:74),'(i4.4,5i2.2)') & wav_year(1),wav_month(1),wav_day(1),wav_hour(1),wav_min(1), & int(wav_sec(1)) elseif (input_mode.eq.4) then call cwav_time_limits(0) call cwav_read_bases wav_nchan=n_cont_trace hyp_flag=.false. call rea_hyp_clear(1) hyp_lat(1)=0. hyp_lon(1)=0. hyp_depth(1)=0. read(cwav_start_time,'(i4,4i2,f2.0)') hyp_year(1),hyp_month(1), & hyp_day(1),hyp_hour(1),hyp_min(1),hyp_sec(1) rea_nphase=0 rea_nhyp=1 c rea_nwav=1 c rea_wav(1)=' '// c & ' ' c write(rea_wav(1)(2:),'(a)') infile(1:seiclen(infile)) c rea_wav(1)(80:80)='6' call systime(time1,time2) rea_id_line=' ACTION:ASI '//time2//' OP: '// & 'STATUS: ID:20021001025954 I' write(rea_id_line(61:74),'(i4.4,5i2.2)') & wav_year(1),wav_month(1),wav_day(1),wav_hour(1),wav_min(1), & int(wav_sec(1)) endif c c clear out phases if set c if (lclear) then rea_nphase=0 endif c c open output files c if (wav_nchan.ge.100) then lwavout= .FALSE. write(*,*) ' too many traces to create output file ' endif c waveform output if (lwavout) then chr_f_form$ = 'unformatted' c call file_ready(outfile) call sei open( unknown$, ! Open file & stop if erro. & ' ', ! Prompt (n/a). & outfile, ! Filename. & write01, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. endif c c number of output traces is twice the input traces, since trigger trace is added c wav_out_nchan = wav_nchan * 2 write(*,*) ' number of traces: ',wav_nchan c c read trace header(s) to make new file header c dist_cnt=0 local_cnt=0 reg_cnt=0 if (input_mode.eq.1) nspec = 0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c start loop over all traces c do i=1,wav_nchan c c original trace is c c c c sta/lta ratio signal is index x c if (lwavout) then ! double traces if cretae output file c=i*2 -1 x=c+1 else c=i x=i endif it=c c c set output variables c wav_out_stat(c)=wav_stat(i) wav_out_stat(x)=wav_stat(i) wav_out_comp(c)=wav_comp(i) wav_out_comp(x)=wav_comp(i) wav_out_year(c)=wav_year(i) wav_out_year(x)=wav_year(i) wav_out_month(c)=wav_month(i) wav_out_month(x)=wav_month(i) wav_out_day(c)=wav_day(i) wav_out_day(x)=wav_day(i) wav_out_hour(c) = wav_hour(i) wav_out_hour(x) = wav_hour(i) wav_out_min(c) = wav_min(i) wav_out_min(x) = wav_min(i) wav_out_sec(c)= wav_sec(i) wav_out_sec(x)= wav_sec(i) wav_out_nsamp(c)= wav_nsamp(i) wav_out_nsamp(x)= wav_nsamp(i) wav_out_rate(c)=wav_rate(i) wav_out_rate(x)=wav_rate(i) cbyte(x)='4' cbyte(c)='4' c c set waveform names if continuous data c if (cwav) then do j=1,cwav_nseg(i) ind=0 do z=seiclen(cwav_filename(i,j)),1,-1 if (cwav_filename(i,j)(z:z).eq.dchar.and. & ind.eq.0) ind=z enddo file_name=cwav_filename(i,j) & (ind+1:seiclen(cwav_filename(i,j))) flag=.true. do c=1,rea_nwav if (rea_wav(c)(2:seiclen(file_name)+1).eq.file_name) & flag=.false. enddo if (flag) then rea_nwav=rea_nwav+1 rea_wav(rea_nwav)= & ' '// & ' 6' write(rea_wav(rea_nwav)(2:seiclen(file_name)+1),'(a)') & file_name(1:seiclen(file_name)) endif enddo endif enddo c c make seisan headers c if (lwavout) then call sheads(wav_out_year,wav_out_month,wav_out_day, * wav_out_hour,wav_out_min,wav_out_sec,wav_out_nchan,c, * net_code,mainhead_text,wav_out_stat,wav_out_comp, * wav_out_nsamp,wav_out_rate,cbyte, * outf,fheader,theader) endif c c write header c if (lwavout) then x=int(((2*wav_nchan)-28)/3) if (x.lt.0) x=0 do z=1,12+x write(write01) fheader(z) enddo endif c c delete files c call systemc('rm -f *.obs',11) event_done=.false. do i=1,wav_nchan c c read response into common block c wav_current_chan(1)=i if (lwavout) then ! double traces if cretae output file c=i*2 -1 x=c+1 else c=i x=i endif write(*,*) 'trace: ',i,' ' // wav_stat(i) // wav_comp(i) it=c c c get parameters for station c c call get_station(wav_out_stat(c),wav_out_comp(c),sta,lta,ratio, & mincoda,dtrle,fill,filh,code) if (code.eq.1) then c c set flag that some station in parameter file c event_done=.true. c c check if P or S c phase=which_phase(phase_mode) c c use P for Q waves if set and P waves used c if (phase.eq.'P') then par_q0=par_p_q0 par_qalpha=par_p_qalpha else par_q0=par_s_q0 par_qalpha=par_s_qalpha endif c c check for phase time c phase_index=0 if (phase_mode.eq.3..or.phase_mode.eq.4..or. & phase_mode.eq.5.) then do j=1,rea_nphase if (rea_stat(j).eq.wav_stat(i).and. & rea_phase(j)(1:4).ne.'SPEC') then if (rea_phase(j)(1:1).eq.phase(1:1).or. c & rea_phase(j)(1:1).eq.phase(1:1).or. & phase_mode.eq.5.) then phase_index=j phase=rea_phase(j)(1:1) call timsec(hyp_year(1),hyp_month(1),hyp_day(1), & rea_hour(j),rea_min(j),rea_sec(j),rea_abs_time(j)) endif endif enddo c c check that phase there if required c elseif ((phase_mode.eq.1..or.phase_mode.eq.2.).and. & phase_infile.eq.1.) then do j=1,rea_nphase if (rea_stat(j).eq.wav_stat(i).and. & rea_phase(j)(1:4).ne.'SPEC') then phase_index=j endif enddo if (phase_index.eq.0.) code=0 endif c c check if spec already done and keep index c spec_index=0 do j=1,rea_nphase if (rea_stat(j).eq.wav_stat(i).and. & rea_phase(j)(1:4).eq.'SPEC') then if (spec_index.eq.0) spec_index=j endif enddo c c read trace, signal stored in signal1 c call wav_read_channel(i) call read_resp if (wav_nsamp(i).gt.max_sample) code=0 endif if (code.eq.1) then do z=1,wav_nsamp(i) idata(z)=int(signal1(z)) c c save raw signal c raw_signal(z)=signal1(z) enddo endif c c make trace header c if (lwavout) then call sheads(wav_out_year,wav_out_month, * wav_out_day,wav_out_hour, * wav_out_min,wav_out_sec,wav_out_nchan,c, * net_code,mainhead_text,wav_out_stat,wav_out_comp, * wav_out_nsamp,wav_out_rate,cbyte, * outf,fheader,theader) write(write01) theader write(write01) (idata(z),z=1,wav_out_nsamp(c)) call sheads(wav_out_year,wav_out_month, * wav_out_day,wav_out_hour, * wav_out_min,wav_out_sec,wav_out_nchan,x, * net_code,mainhead_text,wav_out_stat,wav_out_comp, * wav_out_nsamp,wav_out_rate,cbyte, * outf,fheader,theader) write(write01) theader endif if (code.eq.1) then ! station defined in parameter file c c compute distance if hypocenter given c if (hyp_flag) then call stat_loc(wav_stat(i),data(1)(21:21),slat,slon,sheight) if (slat.eq.0..and.slon.eq.0..and.sheight.eq.0.) then write(*,*) ' location unknown for station: ' // & wav_stat(i) sdist=0. else if (hyp_lat(1).ne.-999.) then xx(1)=slat*degtorad xx(2)=slon*degtorad xx(3)=hyp_lat(1)*degtorad xx(4)=hyp_lon(1)*degtorad call delaz(xx(1),xx(2), & sdist,distdeg,az0,xx(3),xx(4)) write(*,*) ' dist for station ' // wav_stat(i) & // ' in km: ', sdist else sdist=0.1 endif endif endif c c calculations depending on which processes are activated in parameter file c drate=wav_out_rate(c) ! store rate as double precision c c remove linear trend or mean c if (lrtr) then write(*,*) ' removing linear trend' call rtr(signal1,wav_out_nsamp(c),drate) else if (lrmean) then call rmean(signal1,wav_out_nsamp(c)) endif endif c c check if signal is clipped c call check_clipped(signal1,wav_out_nsamp(c),clipped) if (clipped) write(*,*) ' *** signal clipped *** ' c c filter if not lower filter is -1. c if (par_fill(c).ne.-1.) then write(*,*) 'filtering: ',fill,filh c-- calc. filt. cof. call bndpas(fill,filh,1000.0/wav_out_rate(c), * cof,gain) c-- and filter call filter(signal1,wav_out_nsamp(c),cof,gain,npasses) endif do z=1,wav_out_nsamp(c) fil_signal(z)=signal1(z) enddo c c normalize signal to max of 1000 c call normalize_signal(signal1,wav_out_nsamp(c),1000.) c c compute characteristic function, y**2 + k * y'**2 if (lchfu) then call char_fu(signal1,wav_out_nsamp(c),kchfu) endif c c compute sta/lta, ratio signal in signal1 c if (.not.lrec) then lfreeze=.true. call sta_lta_nrec(signal1,wav_out_nsamp(c), * drate,sta,lta,lstasq,ratio,lfreeze) else call sta_lta_rec(signal1,wav_out_nsamp(c), * drate,sta,lta,lstasq) endif c c find phases, ratio given in signal1, spectral analysis to check if signal done c on raw_signal c duration=0. if (lautophase.or. & phase_mode.eq.0..or. & (phase_mode.ge.3..and.phase_index.eq.0.and.lautophase)) then lautoamp=.false. call find_phases(i,signal1,raw_signal,wav_out_nsamp(c), & drate,sta,ratio,dtrle,wav_abs_time(i), & aphases,phase_time,dist_f_sel,fr_type, & spike_test) c c loop over all phases c y=1 if (lall_phases) y=max_phase do z=1,y c c find coda on filtered signal c call find_coda(fil_signal,wav_out_nsamp(c),drate, & wav_abs_time(i),sta, & dtrle,phase_time(z),coda(z),lautocodasq) c c write phases to array, signal has to either above noise spectrum or longer c than minimum coda duration c c if (coda(1).gt.mincoda.or.fr_type(1:1).ne.'N') then if (coda(z).ge.mincoda.and.phase_time(z).ne.0.) then if (aphases(z)(22:22).ne.' ') then c c new phase c lautoamp=.true. rea_nphase=rea_nphase+1 phase_index=rea_nphase call rea_phase_clear(phase_index) rea_auto(phase_index)='autosig' rea_stat(phase_index)=wav_out_stat(c) rea_comp(phase_index)=wav_out_comp(c) rea_co(phase_index)(1:2)=wav_out_comp(c)(1:1)// & wav_out_comp(c)(4:4) rea_phase(phase_index)='P' phase(1:1)='P' rea_onset(phase_index)='E' rea_weight_out(phase_index)=' ' rea_weight_in(phase_index)=' ' rea_polarity(phase_index)=' ' call sectim(phase_time(z),rea_year(phase_index), & jday,rea_month(phase_index),rea_day(phase_index), & rea_hour(phase_index),rea_min(phase_index), & rea_sec(phase_index)) if (coda(z).gt.0..and.coda(z).lt.10000) then rea_coda(phase_index)=coda(z) endif endif endif enddo c c set duration either for P or S, as given in parameter file c if (sdist.ne.0..and.spec_dur_choice.eq.1.) then call gvel_to_time(sdist,phase_time(1), & msec,duration) else duration=spectrum_length_p endif elseif (phase_mode.ge.3..and. & phase_index.ne.0) then c c if phase to be taken from file c phase_time(1) = rea_abs_time(phase_index) if (spec_dur_choice.eq.1.) then call gvel_to_time(sdist,rea_abs_time(phase_index), & msec,duration) else if (which_phase(phase_mode).eq.'P') then phase='P' duration=spectrum_length_p elseif (which_phase(phase_mode).eq.'S') then phase='S' duration=spectrum_length_s endif endif elseif ((phase_mode.eq.1..or.phase_mode.eq.2.) c & .and.phase_time(1).ne.0.) then ! changed lot 20/03/2003 & ) then c c set group velocity window c rea_nphase=rea_nphase+1 phase_index=rea_nphase phase=which_phase(phase_mode) call rea_phase_clear(phase_index) call gvel_to_time(sdist,abs_time, & phase_time(1),duration) rea_auto(phase_index)='autosig' rea_stat(phase_index)=wav_out_stat(c) rea_comp(phase_index)=wav_out_comp(c) rea_co(phase_index)(1:2)=wav_out_comp(c)(1:1)// & wav_out_comp(c)(4:4) rea_phase(phase_index)=' ' rea_phase(phase_index)(1:1)=which_phase(phase_mode) rea_onset(phase_index)='E' rea_weight_out(phase_index)=' ' rea_weight_in(phase_index)=' ' c xxx rea_coda(phase_index)=int(duration) rea_dist(phase_index)=sdist rea_polarity(phase_index)=' ' call sectim(phase_time(1),rea_year(phase_index), & jday,rea_month(phase_index),rea_day(phase_index), & rea_hour(phase_index),rea_min(phase_index), & rea_sec(phase_index)) endif c write(*,*) ' phase_mode, phase_time ',phase_mode,phase_time(1) c c check event type from frequency content c c write(*,*) ' xxx duration: ',duration if (.not.clipped.and.duration.ne.0.) then call freq_event_type(raw_signal,i,dist_f_sel, & phase_time(1), & spectrum_pre_length,duration,fr_type) if (fr_type(1:1).eq.'D') then dist_cnt= dist_cnt+1 elseif (fr_type(1:1).eq.'R'.or.fr_type(1:1).eq.'N') then reg_cnt = reg_cnt + 1 else local_cnt= local_cnt+1 endif endif if (.not.clipped.and.duration.ne.0.) then c c determine auto amplitude c if (lautoamp.and.phase(1:1).eq.'S') then call get_amp_per(i,'WA',raw_signal,phase_time(1), & duration,amp,per,time_amp) c write(*,*) ' xxx ',amp,per,time_amp endif c c add amp and per to common block, add to old phase or make new phase c if (amp.ne.0..and.per.ne.0.) then if ((phase_index.eq.0.or. & phase_mode.eq.1..or.phase_mode.eq.2.).and. & time_amp.ne.0.) then ! if new phase rea_nphase=rea_nphase+1 phase_index=rea_nphase call rea_phase_clear(phase_index) rea_stat(phase_index)=wav_stat(i) rea_comp(phase_index)=wav_comp(i) rea_onset(phase_index)='E' call sectim(time_amp,rea_year(phase_index), & jday,rea_month(phase_index),rea_day(phase_index), & rea_hour(phase_index),rea_min(phase_index), & rea_sec(phase_index)) endif rea_auto(phase_index)='autosig' rea_amp(phase_index)=amp rea_per(phase_index)=per endif c write(*,*) ' xxx autospec ',lautospec c c determine spectral parameters, flat level and corner frequency c if (lautospec) then c c set values from s-file to initial values, used to compare residual c if (spec_index.gt.0) then omi=rea_omega0(spec_index) f0i=rea_cornerf(spec_index) else omi=0. f0i=0. endif write(*,*) ' trying to match spectrum ... ' call get_om_f0(raw_signal,i,spfmin,wav_rate(i)/2.5, & phase_time(1),spectrum_pre_length, & duration,om,f0,ka,omi,f0i, & abs_time,par_q0,par_qalpha,qcorner,par_kappa, & dkappa, & ga_pop_size,ga_n_generations, & ngrid_f,ngrid_om,grid_nloop, & res_m_a,search_method,norm_fit,ffixl,ffixh,fsn1,fsn2) endif res_sum=res_sum+res_m_a(2) proc_count=proc_count+1. write(res_unit,*) proc_count,res_m_a(2) c c write results to common block c if (om.ne.0..and.f0.ne.0.) then rea_nphase=rea_nphase+1 phase_index=rea_nphase rea_stat(phase_index)=wav_stat(i) rea_comp(phase_index)=wav_comp(i) rea_auto(phase_index)='autosig' rea_phase(phase_index)= ' ' rea_phase(phase_index)(1:4)='SPEC' rea_phase(phase_index)(5:5)=phase c rea_phase(phase_index)(5:5)=rea_phase(phase_index)(1:1) c rea_coda(phase_index)=int(duration) rea_omega0(phase_index)=om rea_cornerf(phase_index)=f0 rea_vp(phase_index) = par_vp rea_vs(phase_index) = par_vs rea_q0(phase_index)=par_q0 rea_qalpha(phase_index)=par_qalpha rea_kappa(phase_index)=par_kappa+ka ! changed lot 22/02/2009 rea_density(phase_index)=par_density rea_swin(phase_index)=duration c c source dimension c rea_radius(phase_index)=.37*par_vs/f0 ! changed to 0.36 apr 2013 jh call spec_dist(phase,sdist,hyp_depth(1),geo_corr) rea_geo_dist(phase_index)=1./geo_corr if (sdist.gt.max_dist) max_dist=sdist if (sdist.lt.min_dist) min_dist=sdist c c log moment c rea_moment(phase_index)=4.*pi* & rea_density(phase_index)*10.**rea_omega0(phase_index) & * .833 ! kk & * rea_geo_dist(phase_index) & * (1000.**5) ! to convert units of density and distance & / (1.e9) ! to convert to Nm if (rea_phase(phase_index)(5:5).eq.'P') then rea_moment(phase_index)=rea_moment(phase_index)* & rea_vp(phase_index)**3 else rea_moment(phase_index)=rea_moment(phase_index)* & rea_vs(phase_index)**3 endif rea_moment(phase_index)=log10(rea_moment(phase_index)) rea_mw(phase_index)=2./3.*rea_moment(phase_index) & -6.06 c c xxx remove if loop later c if (spec_index.ne.0) then if (rea_moment(spec_index).ne.0.) then mag_n=mag_n+1 mag_sum=mag_sum+rea_mw(phase_index) mag_save(mag_n)=rea_mw(phase_index) mom_sum=mom_sum+rea_moment(phase_index) c write(*,*) mag_n,mag_sum,rea_mw(phase_index) endif endif c c stress drop c rea_sdrop(phase_index)=(.44*10.**rea_moment(phase_index)) & / (1.e14*rea_radius(phase_index)**3) c c write out comparison c if (rea_moment(spec_index).ne.0.) then write(mom_unit,*) rea_moment(spec_index), & rea_moment(phase_index) write(f0_unit,*) rea_cornerf(spec_index), & rea_cornerf(phase_index) endif endif endif else write(err_unit,*) ' station/comp ',wav_out_stat(c),' ', & wav_out_comp(c), & ' not defined in autosig.par' endif c c write out data c if (lwavout) then do z=1,wav_out_nsamp(c) idata(z)=int(signal1(z)) enddo write(write01) (idata(z),z=1,wav_out_nsamp(c)) endif enddo ! end of loop over all traces cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (mag_n.gt.0) then c c get variance c do i=1,mag_n yy=yy+(mag_save(i)-mag_sum/float(mag_n))**2 enddo yy=(1./float(mag_n-1)*yy)**.5 write(mag_unit,*) rea_av_mw,mag_sum/float(mag_n),yy write(65,*) mom_sum/float(mag_n) if (mag_ml_index.ne.0) then write(66,*) hyp_mag(mag_ml_index,1),mag_sum/float(mag_n) endif endif if (dist_cnt.gt.local_cnt.and.dist_cnt.gt.reg_cnt) then fr_type='D ' write(*,*) ' *************** Distant *********** ' elseif (local_cnt.gt.dist_cnt.and.local_cnt.gt.reg_cnt) then fr_type='L ' write(*,*) ' *************** Local ************* ' else fr_type='R ' write(*,*) ' *************** Regional ********** ' endif if (input_mode.ne.2) then c data(1)(22:22)=fr_type(1:1) hyp_dist_id(1)=fr_type(1:1) endif c c compute averages and standard deviation and store in data c call rea_event_out(0,.true.,data,code) ! update array call update_spec(data,rea_nrecord,rea_nhead,agency,test,.false.) c c write out nordic file c do i=1,rea_nrecord write(write02,'(a)') data(i) enddo c c write out event if at least 1 channel processed c if (event_done) then do i=1,nrecord_save write(event_unit,'(a80)') data_save(i) enddo endif c c call gmt script c c change if (input_mode.ne.2.and.input_mode.ne.4) then call systemc('spec.exp',8) endif if (lwavout) then call sei close( close$, write01, code) endif write(*,*) ' min distance = ',min_dist write(*,*) ' max distance = ',max_dist if (input_mode.eq.3) then goto 60 endif if (input_mode.eq.4) then c c calculate new start time c read(cwav_start_time,'(i4,4i2,f2.0)') year,month,day,hour, & min,rsec call timsec(year,month,day,hour,min,rsec,msec) msec=msec+(cont_interval*.9) call sectim(msec,year,jday,month,day,hour,min,rsec) write(cwav_start_time,'(i4.4,4i2.2,i2.2)') & year,month,day,hour,min,int(sec) write(*,*) ' new start time: ',cwav_start_time read(end_time,'(i4,4i2,f2.0)') year,month,day,hour, & min,sec call timsec(year,month,day,hour,min,rsec,msec2) if (msec.ge.msec2) goto 9999 goto 60 endif if (in.eq.1) then goto 50 endif 9999 continue write(*,*) ' residual ', & res_sum/proc_count if (.not.onlyonce) then input_mode=1 goto 40 endif c c close files c if (in.eq.1) call sei close( close$, read02, code) call sei close( close$, err_unit, code) call sei close( close$, trace_unit, code) call sei close( close$, event_unit, code) call sei close( close$, mag_unit, code) call sei close( close$, f0_unit, code) call sei close( close$, mom_unit, code) call sei close( close$, res_unit, code) call sei close( close$, write02, code) call close_post call sei close( close$, plotunit, code) stop end