c - program to get waveform data from all waveform files in one s-file, c or a list ow waveform file in a file of filenr.lis format c the data can be selected in verious ways like channels, times etc c - the information is written out to a file in the local directory c - in seisweb mode, data is written to a file specified c and at end a seisweb message written to screen c - output format is SEISAN, GSE and SAC c - program can also only give file info c - optinally points can be skipped so resampled files can also be made, c it is up to the user to use correct antialias filters c - filter option is a 4 poles butterworth running one way, ther si currently c no check if the output after filtering is larger than one so funny c results might occur in the write out c c arguments: c -seisweb: use with seisweb, seisweb screen message C c -chan_out_file: file name of file with channel description c if this option is given, program terminates c after writing out the file c -wav_out_file: complete filename including path to wav output, c if not set, output to extract.out c If set to SEISAN, a seisan standard output file name c will be used, net code from SEISAN.DEF c -sfile : s-filename with complete path c -wav_files: list of waveform files, filenr.lis format c overrides input from -sfile c -format : output format, SEISAN, GSECM6, GSEINT, MSEED, SEED c default is SEISAN c -maxpoints : skip points to get approximately maxpoints c -filter : bandpass filter limits c -ground : compute displacement, velocity or acceleration (1,2,3) c -chansel: file with selection parameters, if -chansel not used, c all channels selected with whole window. c Format: First line gives number of channels, free form c Each following line give channel #, start and c duration of selected channels. if start is c zero use first sample, if interval is zero, c use whole interval, free format c -start : (1) Start time relative earliest channnel, c (2) Abs start time string yy...s.sss for all channels c IF ABS time used, string MUST contain a '.' c (3) Abs start time string yyyymmddhhmmss (integer), c used to define cont start time c -duration: length of all channels to select c start and duration not used if chansel file given c -duration will not have an effect if -start c is not used. c c -ichan : selection of one channel c -npole n: number of poles used for filter c -stat_out: write out station location file c -resp_out: write out list of response files c -cwav : read from cont data base c -cseed : interval from large seed volume c -seed_location: seed location code NOT USED c -wav_in_file: input waveform file c -interactive : enable interactive mode c -rsam: comput 1 min rsam data c -cbase: name of file with selected continuous databases, c all is default c c examples: c c wavetool -sfile sfile c all channels selected from all files c wavetool -sfile sfile -start 100 c start 100 secs after earliest channel, take all channels and all c remaining data c wavetool -sfile sfile -start 19950101202211.1 c same as above except now start time is in absolute time c wavetool -wav_in_file lsa.iris.big.seed -start 19990509000000 -duration 500 -cseed c extract from a seed volume c c c changes: c c april 19, 00 lo fix gse output c april 20, 00 lo add some arg options and fix gse output c april 25 jh many changes c may 3 lo routine signal_to_visual c may 8 00 lo some fixes with the chansel input file c may 10 jh set network code to MERGE if no code defined c may 15 lo some smaller changes, bug net code(jh) c may 25 bmt add ichan option (select 1 channel) c oct 20 jh merge_wav changed to 5 chars c feb 01 01 lo ground motion output c mar 05 01 lo use i3 in channel list c mar 06 01 lo keep time of max and min in signal_to_visual c may 01 01 lo started to extract from waveform files directly c july 13 01 lo extract now also wave conversion, any input, any output c nov 11 01 jh extract from cont base c dec 03 01 bm clean chan_out_file when nothing is found c bm added new option "-command_file" for reading prompt arguments c from a file, one line per argument c mar 14 02 jh comment out write(17 c april 30 02 lo seisan output always 4 byte c may 8 02 bm multiply by 10 when output is GSE format (only seisweb) c aug 22 02 lo added option -wav_in_file to give waveform file nam c apr 9 03 jh put in gain in seisan output if filter or ground correction to c make sure values are more than 1 count, make it possibele to filter c also when no response removal, check if filtering is possible c may 14 03 jh add text to header when corrected output written c jun 30 03 lo add rsam output c sep 10 03 jh check value of sample in GSE format, if response, was sometimes too big c apr 04 lo: some format fixing c dec 2 jh: large seed volume c mar 10 05 jh: add component to output file name if only one channel c jun 20 05 jh: seed location code input c feb 17 06 jh: change references to extract. NOTE: seisweb uses old c extract !!!!!!! c mar 28 06 jh: bug with overwrite seisan file c bugs in defining waveform file name etc pointed out c by wayne crawford c mar 30 06 lo: option to select continuous database, used from mulplt c jan 17 07 jh: fix bug when reading seed files with time gap, if converting c to seisan file, main header will be wrong, but channel c header ok, to be fixed c mar 14 07 lo: added seed output using gse2seed c mar 21 07 lo: added option for dataless seed c may 30 07 lo: use file to read cbase c oct 24 07 jh: put in writing of location and networ in seisan and mseed c remove old INPUT OF LOCATION CODE, should now come from c input file c dec 19 07 jh fix probelm of wrong interval using selection out of a c seed file, cseed in common block c aug 01 08 lo set wav_interactive to true if filenr.lis c aug 27 08 jh put in more comments c oct 23 08 jh fix output for seisweb from chan_out_text c april 17 09 jh remove option s-file form option cseed, should not have c been there and made cseed option not work implicit none include 'seidim.inc' include 'libsei.inc' include 'waveform.inc' include 'seisan.inc' include 'write_mseed.inc.f' integer seiclen integer nstat,nphase,nhead,nrecord,id ! for reading s-file character*1 type,exp ! ----------------- character*80 sfilname ! s-file name character*80 chan_out_file ! file name for channel info character*80 data(max_data) ! one s-file character*80 chansel ! file with selection integer npresent ! number of wav files present character*80 text ! general text string character*80 wav_out_file ! full name of output file character*80 wav_files ! input files with wav f. names character*80 file_out ! output file name character*80 seed_file_out ! output file name character*5 net_code ! network code, for output file name integer narg ! number of arguments character*80 arg(40) ! arguments logical seisweb ! true if call from seisweb logical exist ! file existance logical station_out_flag ! true if write out station file logical resp_out_flag ! true if write out list of c logical cseed now in common block ! large seed file ! resp files c character*2 seed_location_code ! seed location code character*80 stat_list(max_trace) ! station list integer stat_cnt ! station counter character*80 format_out ! output format integer max_points ! maximum number of points out integer nskip(max_trace) ! number of points to skip integer syear,smonth,sday,shour,smin ! common start real ssec ! common start integer code ! return error code real flow,fhigh ! filter band pass f integer npole ! number of poles real cof(8) ! filter coefficients real gain ! filter gain logical ground ! true if ground motion integer ground_type ! dis=1, vel=2 or acc=3 complex y_com(max_sample) ! complex spectrum integer file_cnt ! file counter real slat,slon,sheight ! station location character*1 resp_flag ! Y if response there character*1 statfile_char ! character in station file character*120 chan_out_text(1000) ! output text array real dc,ndc ! used for remove_dc character*80 question ! question for filenames logical sfile_flag ! true if sfile logical rsam_flag ! true for rsam output character*5 cbase_select(max_cont) ! selected databases integer n_cbase_select ! number of the above character*80 cbase_list ! file with cont databases c c following for seisan format output files c integer year(max_trace),month(max_trace),day(max_trace), * hour(max_trace),min(max_trace),nsamp(max_trace) real rate(max_trace),sec(max_trace) character*80 mainhead(max_trace) ! seisan main header character*1040 chahead ! seisan channel header character*29 h_text ! seisan header integer signal4b(max_sample) ! 4 byte seisan data integer*2 signal2b(max_sample) ! 2 byte seisan data character*1 cbyte(max_trace) ! indicator of 2 or 4 bytes real ymax ! maximum output value real factor ! gain factor if used c c following parameters for selected output c double precision start,duration ! start and duration one channel character*18 pstart ! absolute start from prompt,ymdhms double precision abs_start_time ! pstart in seconds real pduration ! duratiopn from prompt character*5 stat(max_trace) character*4 comp(max_trace) character*2 location(max_trace) ! seed location character*2 network(max_trace) ! seed network integer i,k,j,l,ichan,in character*80 file real evla,evlo,evel ! event loctaion integer nerr ! error code logical open_file,op_file ! for opening mseed file integer mseed_block ! block counter for wring miniseed files character*160 system_call logical sun,linux,pc 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 call computer_type(sun,pc,linux) chan_out_file=' ' sfilname=' ' wav_files=' ' wav_out_file=' ' seisweb=.false. c seed_location_code=' ' pstart='0' flow=0.0 fhigh=0.0 npole=4 pduration=0.0 n_cbase_select=0 chansel='_' format_out='SEISAN' ! default format to write is seisan evla=0. evlo=0. evel=0. ichan=-1 ground=.false. ground_type=0 in=0 statfile_char=' ' stat_cnt=0 sfile_flag=.false. station_out_flag=.false. resp_out_flag=.false. max_points=0 wav_interactive=.false. rsam_flag=.false. cseed=.false. cwav=.false. mseed_block=1 ! start with first block cbase_list=' ' c c get seisan defaults for network code merge_wav c call get_seisan_def c c get arguments c call get_arguments(narg,arg) 10 continue if(narg.lt.2) then c c get s-file name or waveform file name c question=' Filename of s-file or waveform file, ' // & 'number or filenr.lis ' call filename(question,sfilname) if (seiclen(sfilname).le.0.or.sfilname.eq.'EOF') stop c c filenr.lis, interactive input c if (sfilname(1:10).eq.'filenr.lis') then open(2,file='filenr.lis',status='old', err=999) in=1 wav_interactive=.true. ! lot 1/8/08 file_cnt=0 11 continue read(2,'(7x,a)',end=12) sfilname if(sfilname(1:1).eq.' ') go to 12 goto 13 12 continue close(1) goto 10 endif 13 continue file_cnt=file_cnt+1 c c check if wav file exist c if(sfilname.ne.' ') then inquire(file=sfilname,exist=exist) if(.not.exist) then write(6,*) ' No such file' goto 10 endif endif c c interactive input of parameters c text=' ' if (in.eq.0) then write(6,*) * 'Maximum number of points in output trace, return for all' read(5,'(a)') text endif if(text.eq.' ') then max_points=0 else read(text,'(i8)') max_points endif text=' ' if (file_cnt.eq.1) then write(6,*) ' Ground motion output (dis = 1, ' & // 'vel = 2, acc = 3, none = return)' read(5,'(a)') text if(text.eq.' ') then ground_type=0 else read(text,'(i1)') ground_type if (ground_type.eq.1.or. & ground_type.eq.2.or. & ground_type.eq.3) then ground = .true. endif c c read filter c c write(*,*) ' filter low and high ' c read(5,*) flow,fhigh endif c c read filter c write(*,*) ' Filter low and high, return for no filter ' read(5,'(a)') text if(text.eq.' ') then flow=0.0 fhigh=0.0 else call sei get values(2,text, code ) flow=array$(1) fhigh=array$(2) endif endif else c c read argument from a file c if( arg(1)(1:13) .eq. '-command_file' ) then call get_argument_from_file(arg,narg) endif c c read arguments to fill in choices c do i=1,narg if(arg(i)(1:13).eq.'-wav_out_file') then wav_out_file=arg(i+1) endif if(arg(i)(1:6).eq.'-sfile') then sfilname=arg(i+1) sfile_flag=.true. endif if(arg(i)(1:6).eq.'-cbase') then c n_cbase_select=n_cbase_select+1 c cbase_select(n_cbase_select)=arg(i+1) cbase_list=arg(i+1) endif if(arg(i)(1:7).eq.'-filter') then call sei get values(1,arg(i+1), code ) flow=array$(1) call sei get values(1,arg(i+2), code ) fhigh=array$(1) endif if(arg(i)(1:5).eq.'-rsam') then rsam_flag=.true. endif if(arg(i)(1:6).eq.'-npole') then call sei get values(1,arg(i+1), code ) npole=int(array$(1)) endif if(arg(i)(1:7).eq.'-ground') then call sei get values(1,arg(i+1), code ) ground_type = int(array$(1)) if (ground_type.eq.1.or. & ground_type.eq.2.or. & ground_type.eq.3) then ground = .true. endif endif if(arg(i)(1:8).eq.'-seisweb') then seisweb=.true. endif if(arg(i)(1:9).eq.'-stat_out') then station_out_flag=.true. endif if(arg(i)(1:9).eq.'-resp_out') then resp_out_flag=.true. endif if(arg(i)(1:10).eq.'-wav_files') then wav_files=arg(i+1) endif if(arg(i)(1:7).eq.'-format') then format_out=arg(i+1) endif if(arg(i)(1:8).eq.'-chansel') then chansel=arg(i+1) endif if(arg(i)(1:14).eq.'-chan_out_file') then chan_out_file=arg(i+1) open(1,file=chan_out_file,status='unknown') endfile 1 rewind 1 close(1) endif if(arg(i)(1:6).eq.'-start') then pstart=arg(i+1) endif if(arg(i)(1:9).eq.'-duration') then call sei get values(1,arg(i+1), code ) ! Extract 1 value pduration=array$(1) endif if(arg(i)(1:10).eq.'-maxpoints') then call sei get values(1,arg(i+1), code ) ! Extract 1 value max_points=array$(1) endif if(arg(i)(1:6).eq.'-ichan') then call sei get values(1,arg(i+1), code ) ! Extract 1 value ichan=array$(1) endif if(arg(i)(1:5).eq.'-cwav') then cwav=.true. ! cont base endif if(arg(i)(1:6).eq.'-cseed') then cseed=.true. ! large seed endif c if(arg(i)(1:6).eq.'-seed_location') then c seed_location_code=arg(i+1)(1:2) ! seed location code c endif if(arg(i)(1:12).eq.'-interactive') then wav_interactive=.true. endif if(arg(i)(1:12).eq.'-wav_in_file') then sfilname=arg(i+1) endif enddo endif c c-------------------------------------------------------- c--------------------------------------------------------- c get all info if cont data base or a large seed file c--------------------------------------------------------- c--------------------------------------------------------- c if(cwav.or.cseed) then if(cseed) write(6,*)' Large SEED' c c check if cont databases c if(cwav) then if (n_cont_base.eq.0) then write(6,*) ' No continuous database defined in SEISAN.DEF ' stop endif c c read file with cont databases c if (cbase_list.ne.' ') then open(1,file=cbase_list,status='old') 18 continue n_cbase_select=n_cbase_select+1 c write(*,*) n_cbase_select read(1,'(a)',end=19) cbase_select(n_cbase_select) c write(*,*) ' database selected: '// c & cbase_select(n_cbase_select) goto 18 19 continue close(1) n_cbase_select=n_cbase_select-1 endif c c check if databases selected c if (n_cbase_select.gt.0) then do i=1,n_cbase_select cont_base(i)=cbase_select(i) enddo n_cont_base=n_cbase_select endif c c write out the databases c write(*,*)'Databases are:' do i=1,n_cont_base write(*,*)' ',cont_base(i) enddo endif c c get start time and time interval c if(pduration.eq.0.0.or.pstart.eq.' ') then write(6,*) * 'Abs start time and interval must be given for'// * ' cwav or cseed option' stop endif cwav_start_time=pstart ! also used for seed pstart='0' ! indicate that whole window used, no further ! selection in time cont_interval=pduration c c calculate end time and extended start time c call cwav_time_limits(0) c c read the header information for all files in all bases in time c interval, assume info available in common block, not large seed file c if(cwav) call cwav_read_bases c c read data for section of seed file c if(cseed) then wav_filename(1)=' ' wav_filename(1)=sfilname wav_nfiles=1 ! only one file can be used write(*,*) ' reading waveform headers ... ' write(6,*) wav_filename(1) call wav_init call read_wav_header(1) write(6,*) 'headers read' call wav_seed_read_header_interval ! find position in file of window write(6,*) ' large seed read' endif c c for continuous/seed data, if number of samples not set, do it here c do k=1,wav_nchan if (wav_nsamp(k).le.0) then write(6,*)'set nsamp' wav_nsamp(k)=int(pduration*wav_rate(k)) wav_duration(k)=pduration endif enddo c c terminate if no waveform files c if(wav_nchan.eq.0) goto 999 goto 700 ! jump all normal reading, straight to output endif c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c ccccccccccccccc end of section for large seed or cont data cccccccccccc c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c----------------------------------------------------------------------- c check if waveform file by reading the header c this is done because sfilname can be a waveform file name c and a sfile name, should be changed ! c at this point, the waveform file name is transferred from c sfilname to wav_file_name c----------------------------------------------------------------------- if (.not.sfile_flag.and.sfilname.ne.' ') then wav_error_message=' ' call wav_init wav_nfiles=1 wav_filename(1)=' ' wav_filename(1)=sfilname call get_full_wav_name(sfilname, wav_filename(1)) c write(*,*) ' reading waveform headers ... ' call read_wav_header(1) endif c c--------------------------------------------------------- c case of waveform files from sfile c--------------------------------------------------------- c c check if sfile, if not waveform c if(wav_error_message.ne.' ') then if(sfilname.ne.' ') then ! assume s-file open(1,file=sfilname,status='old', err=999) c c read s-file c call indata * (1,nstat,nphase,nhead,nrecord,type,exp,data,id) close(1) c c read location and depth from first line c read(data(1)(24:43),'(f7.3,f8.3,f5.2)') evla,evlo,evel read(data(1)(16:16),'(a1)') statfile_char c c get waveform file names c call auto_tr(data,nhead,wav_nfiles,wav_filename) endif endif c c---------------------------------------------------------------- c case of input from a list of waveform files c--------------------------------------------------------------- c c check for waveform files in a filenr.lis format input file c wav_files is from argument list c if(wav_files.ne.' ') then write(6,*) wav_files open(1,file=wav_files,status='old', err=999) j=0 20 continue j=j+1 read(1,'(7x,a)',end=21) wav_filename(j) if(wav_filename(j).ne.' ') go to 20 21 continue wav_nfiles=j-1 endif c c---------------------------- c if(.not.seisweb) then write(6,*)' Number of wav-files',wav_nfiles endif c c find how many waveform files are present, replace the origial c names with the present files full path c npresent=0 ! no files initially do i=1,wav_nfiles call get_full_wav_name(wav_filename(i),text) if(.not.seisweb) then c write(*,*) i,wav_filename(i),text endif if(text.ne.' ') then npresent=npresent+1 wav_filename(npresent)=text endif enddo wav_nfiles=npresent if(.not.seisweb) then write(6,*)' Number of wav-files present', wav_nfiles endif c c terminate if no waveform files c if(wav_nfiles.eq.0) goto 999 c wav_nchan=0 c c loop to read all headers of all files c do i=1,wav_nfiles c write(6,'(1x,a)') wav_filename(i) call read_wav_header(i) if(wav_error_message.ne.' ') write(6,'(1x,a)') * wav_error_message write(6,'(1x,a40,2x,a)') wav_filename(i)(1:40), * wav_file_format(i) enddo c if(.not.seisweb) then write(6,*)' Total number of channels available:',wav_nchan endif c c terminate if no channels c if(wav_nchan.eq.0) goto 999 c--------------------------------------------------------------- c section for output, if cont data base enter here directly c since all headers read above c--------------------------------------------------------------- c 700 continue ! from cont reading above c c write input out, optionally also in a file c if(.not.seisweb) then write(6,*)' Total duration:',wav_total_time open(3,file='respfile_list.out',status='unknown') do i=1,wav_nchan c c check if response available for channel and get station location c c j was not set in 2. "if(station_out..)", pointed out by c wayne crawford, mar 06, so in following line j=1 so j is set c when station_out is false c j=1 c if (station_out_flag) then call stat_loc(wav_stat(i),statfile_char,slat,slon,sheight) c c sort into station location array c j=1 do while(j.le.stat_cnt.and. & wav_stat(i).ne.stat_list(j)(1:5)) j=j+1 enddo endif if (station_out_flag.and. & wav_stat(i).ne.stat_list(j)(1:5)) then stat_cnt=stat_cnt+1 stat_list(stat_cnt)=' ' write(stat_list(stat_cnt),'(a5,1x,f8.3,1x,f8.3)') & wav_stat(i),slat,slon endif resp_flag='N' if (ground.or.resp_out_flag) then wav_current_chan(1) = i wav_resp_file = ' ' c xxx call read_resp if (wav_resp_filename.ne.' ') then resp_flag='Y' write(3,'(a)') & wav_resp_filename(1:seiclen(wav_resp_filename)) c write(3,'(a)') 'cp /net/seismo/seismo/CAL/'// c & wav_resp_filename(1:seiclen(wav_resp_filename)) //' .' else c write(555,'(a)') 'cp /net/seismo/seismo/CAL/'// c & wav_resp_file(1:seiclen(wav_resp_file)) //' .' endif wav_resp_file = ' ' endif c write(6,'(1x,i2,1x,a,1x,a,1x,i5,4i3,f7.3,i6,3f9.3)') c 05 march 2001, lo c 08 June 2001 lo, added location and resp c chan_out_text(i)=' ' write(chan_out_text(i), * '(i3,1x,a,1x,a,1x,i5,4i3,f7.3,i6,3f9.3, * 2(1x,f8.3),1x,f5.0,1x,a1)') * i, wav_stat(i),wav_comp(i), * wav_year(i),wav_month(i),wav_day(i),wav_hour(i),wav_min(i), * wav_sec(i),wav_nsamp(i),wav_rate(i),wav_delay(i), * wav_duration(i), * slat,slon,sheight,resp_flag write(*,'(a)') chan_out_text(i)(1:seiclen(chan_out_text(i))) enddo close(3) ! close response list file if (station_out_flag) then open(3,file='station_list.out',status='unknown') open(4,file='station_list.err',status='unknown') do j=1,stat_cnt read(stat_list(j),'(6x,f8.3,1x,f8.3)') slat,slon if (slat.ne.0..and.slon.ne.0.) then write(3,'(a)') stat_list(j)(1:seiclen(stat_list(j))) else write(4,'(a)') stat_list(j)(1:seiclen(stat_list(j))) endif enddo close(3) close(4) endif endif if(seisweb) then do i=1,wav_nchan chan_out_text(i)=' ' write(chan_out_text(i), * '(i3,1x,a,1x,a,1x,i5,4i3,f7.3,i6,3f9.3, * 2(1x,f8.3),1x,f5.0,1x,a1)') * i, wav_stat(i),wav_comp(i), * wav_year(i),wav_month(i),wav_day(i),wav_hour(i),wav_min(i), * wav_sec(i),wav_nsamp(i),wav_rate(i),wav_delay(i), * wav_duration(i), * slat,slon,sheight,resp_flag enddo endif c c read channel selection file c if(chan_out_file.ne.' ') then open(1,file=chan_out_file,status='unknown') do i=1,wav_nchan c write(1,'(1x,i2,1x,a,1x,a,1x,i5,4i3,f7.3,i6,3f9.3)') c 05 march 2001, lo write(1,'(a)') chan_out_text(i)(1:seiclen(chan_out_text(i))) c write(1,'(i3,1x,a,1x,a,1x,i5,4i3,f7.3,i6,3f9.3)') c write(1,'(i3,1x,a,1x,a,1x,i5,4i3,f7.3,i6,3f9.3, c * 2(1x,f8.3),1x,f5.1,1x,a1)') c * i, wav_stat(i),wav_comp(i), c * wav_year(i),wav_month(i),wav_day(i),wav_hour(i),wav_min(i), c * wav_sec(i),wav_nsamp(i),wav_rate(i),wav_delay(i), c * wav_duration(i), c * slat,slon,sheight,resp_flag enddo close(1) if (wav_out_file.eq.' ' ) then if(seisweb) then write(6,*)'extract channel info terminated: ok ', wav_nchan endif stop endif endif c c------------------------------------------------------------------------- c select and possibly process data for output file, can be done c from keyboard if no input on prompt line c following is what is selected using keyboard input c------------------------------------------------------------------------- c c select channels to read, by default all data are selected c call wav_index_total do i=1,wav_nchan wav_out_chan(i)=i wav_out_start(i)=wav_abs_time(i)-wav_abs_time(wav_first) wav_out_duration(i)=wav_duration(i) c write(6,*) 'org wav duration', wav_duration(i) enddo wav_out_nchan=wav_nchan c c----------------------------------------------------------------------- c keyboard input c----------------------------------------------------------------------- c if(narg.lt.2) then if(in.ne.1) then write(6,*)' Select all data, y=return,n' read(5,'(a)') text if(text.ne.' '.and.text.ne.'y'.and.text.ne.'Y') then write(6,*)' Number of channels to read' read(5,*) wav_out_nchan write(6,*)' Channel number, start time and window' do i=1,wav_out_nchan read(5,*,err=555) wav_out_chan(i),start,duration c c find out if a small or a large number, which indicate if abs start c if(start.gt.1.0e10) then write(pstart,'(f18.3)',err=555) start read(pstart,'(i4,4i2,f6.3))',err=555) * syear,smonth,sday,shour,smin,ssec goto 556 555 continue write(6,*)' Error in input' stop 556 continue call timsec(syear,smonth,sday,shour,smin,ssec, * abs_start_time) start=abs_start_time-wav_abs_time(wav_first) ! cal. rel start endif wav_out_start(wav_out_chan(i))=start ! put in correct channel if(duration.ne.0.0) then wav_out_duration(wav_out_chan(i))=duration endif enddo endif endif if ((in.eq.1.and.file_cnt.eq.1).or.(in.eq.0)) then if (sun) then write(*,*) *' Output formats '// *'(SEISAN, GSE (def=CM6), GSEINT, SAC, MSEED, SEED, SEEDDL) ' else write(*,*) *' Output formats '// *'(SEISAN, GSE (def=CM6), GSEINT, SAC, MSEED) ' endif write(*,*) ' Default is SEISAN=return' read(*,'(a)') format_out call casefold(format_out) if(format_out.eq.' ') format_out='SEISAN' if(format_out.eq.'SEED'.and..not.sun) then write(*,*) ' SEED output format only works on sun ' stop endif c if(format_out.eq.'MSEED'.and.seed_location_code.eq.' ') then c write(6,*)'SEED location code, 2 characters' c read(5,'(a2)') seed_location_code c endif endif else c c------------------------------------------------------------------------------ c selection is from prompt, if no selection file, all data selected initially c there could also be a start and duration, see below c------------------------------------------------------------------------------ c c------------------------- c option selection file c------------------------- c if(chansel.ne.'_') then ! from a file open(1,file=chansel,status='old',err=901) goto 902 901 write(6,*)' Channel selection input file missing' goto 999 902 continue c c first get number of channels c read(1,*,err=903,end=903) wav_out_nchan c c read info if abs time used c do i=1,wav_out_nchan read(1,*,err=903) wav_out_chan(i),start,duration c c find out if a small or a large number, which indicate if abs start c start is used to select the start of the window in abs time c if(start.gt.1.0e10) then ! start is abs start write(pstart,'(f18.3)') start read(pstart,'(i4,4i2,f6.3))',err=903) * syear,smonth,sday,shour,smin,ssec call timsec(syear,smonth,sday,shour,smin,ssec, * abs_start_time) start=abs_start_time-wav_abs_time(wav_first) ! cal. rel start wav_out_start(wav_out_chan(i))=start else if (start.eq.0.) then ! use first sample wav_out_start(wav_out_chan(i))= * -wav_abs_time(wav_first)+wav_abs_time(wav_out_chan(i)) ! relative first channel else wav_out_start(wav_out_chan(i))=start endif endif if(duration.ne.0.0) then wav_out_duration(wav_out_chan(i))=duration endif enddo close(1) goto 904 903 continue write(6,'(a,a)')' Error in channel select file' goto 999 904 continue else c c---------------------------------------------------------------- c option same start and window all channels as given on prompt c----------------------------------------------------------------- c if(pstart.ne.'0') then read(pstart,'(f20.3)')start if(start.gt.1.0e10) then ! abs start times write(pstart,'(f18.3)') start read(pstart,'(i4,4i2,f6.3))',err=950) * syear,smonth,sday,shour,smin,ssec call timsec(syear,smonth,sday,shour,smin,ssec, * abs_start_time) if(abs_start_time.lt.wav_abs_time(wav_first)) then write(6,*)' Start is before first data point, ', * ' will use first point' abs_start_time=wav_abs_time(wav_first) endif start=abs_start_time-wav_abs_time(wav_first) ! cal. rel start else call sei get values(1,pstart, code ) ! Extract 1 value start=array$(1) endif goto 951 950 continue write(6,*)' Error in time string' goto 999 951 continue do i=1,wav_out_nchan wav_out_start(i)=start if(pduration.ne.0.0) wav_out_duration(i)=pduration enddo endif c c---------------------------------------------------------------- c option for select 1 channel c----------------------------------------------------------------- c if(ichan.ne.-1) then wav_out_nchan=1 wav_out_chan(1)=ichan endif endif if(pstart.ne.' '.and.pduration.ne.0.0) then endif endif c c------------------------------------------------------------------------- c calculate new start times and intervals, check if data is available c the call to wav_get_interval also calculates new number of samples c in variable wav_out_nsamp. wav_out_start is input and gives start time relative c to first sample in file, wav_out duration gives the lenght of the window. c both set above c------------------------------------------------------------------------- c c write(6,*) 'wav-out-start,duration',wav_out_start(1), c *wav_out_duration(1) call wav_get_interval c c get rid of empty channels c k=0 do i=1,wav_out_nchan if(wav_out_status(wav_out_chan(i)).gt.0) then k=k+1 wav_out_chan(k)=wav_out_chan(i) endif enddo wav_out_nchan=k c c check if any data left c cc write(6,*)' Number of channels after time window selection', cc *wav_out_nchaN if(wav_out_nchan.eq.0) then write(6,*)' No data left after window selection' goto 999 endif c if(.not.seisweb) then c do i=1,wav_out_nchan c k=wav_out_chan(i) c nskip(k)=1 ! not less than every point c write(6,'(i5,1x,a5,a4,i4)')k,wav_stat(k),wav_comp(k), c * wav_out_status(k) c enddo c endif c c find if skipping, if so correct number of points and sample rate c xxx c if(max_points.gt.0) then c do i=1,wav_out_nchan c k=wav_out_chan(i) c nskip(k)=int(wav_out_nsamp(k)/max_points) c if(nskip(k).lt.1) nskip(k)=1 c enddo c endif if(max_points.gt.0) then do i=1,wav_out_nchan call signal_to_visual(wav_out_chan(i),max_points,1) enddo endif c c set headers for rsam data c if (rsam_flag) then do i=1,wav_out_nchan call signal_to_rsam(wav_out_chan(i),1) enddo endif c c------------------------------------------------------------------------- c section to write out, different formats possible c------------------------------------------------------------------------- c c Make main header if needed and open output file, decide name of c output file c c if(wav_out_file.ne.' '.and.wav_out_file.ne.'SEISAN') then !lo 2004 file_out=wav_out_file c elseif (format_out.eq.'SEISAN') then ! lo 11-18-2002 c file_out=format_out elseif (wav_out_file.eq.'SEISAN') then file_out='SEISAN' elseif (format_out.ne.' '.and.narg.lt.2) then file_out=format_out else file_out='wavetool.out' endif c c fix sample rate and number of samples according to skipping c xxx c do i=1,wav_out_nchan c k=wav_out_chan(i) c wav_out_rate(k)=wav_out_rate(k)/nskip(k) c wav_out_nsamp(k)=wav_out_nsamp(k)/nskip(k) c enddo net_code='MERGE' c c take net code from SEISAN.DEF c if(merge_wav.ne.' ') then net_code=merge_wav endif c c if all channels have the same station name, use that for net code c text(1:5)=wav_stat(wav_out_chan(1)) do i=1,wav_out_nchan k=wav_out_chan(i) if(text(1:5).ne.wav_stat(k)) goto 888 enddo net_code=text(1:5) ! use station code for net code 888 continue c c c c make main header, the out variabels are used since the file might have been c cut so no longer the read start time c c do i=1,wav_out_nchan k=wav_out_chan(i) year(i)=wav_out_year(k) month(i)=wav_out_month(k) day(i)=wav_out_day(k) hour(i)=wav_out_hour(k) min(i)=wav_out_min(k) sec(i)=wav_out_sec(k) stat(i)=wav_stat(k) comp(i)=wav_comp(k) location(i)=wav_location(k) network(i)=wav_network(k) c xxx rate(i)=wav_out_rate(k)/nskip(k) c nsamp(i)=wav_out_nsamp(k)/nskip(k) rate(i)=wav_out_rate(k) nsamp(i)=wav_out_nsamp(k) c cbyte(i)=wav_cbyte(k) cbyte(i)='4' ! lot 30 April 2002 enddo c c make seisan main header, also used to make file name c call sheads(year,month,day,hour,min,sec,wav_out_nchan,1, * net_code,h_text,stat,comp, * nsamp,rate,cbyte, * file,mainhead,chahead) c c if only one chan, add component so no overwrite c if(wav_out_nchan.eq.1) then text=file file=text(1:seiclen(text))//comp(1) do i=1,seiclen(file) if(file(i:i).eq.' ') file(i:i)='_' enddo endif c---------------------- c SEISAN format c---------------------- if(format_out.eq.'SEISAN') then h_text='Extracted waveform file ' ! title c write(6,'(1x,a)') file if(file_out.eq.'SEISAN') file_out=file c c next 3 lines seems a repeat of above ??? comm out mar 28 ,06 jh c if (wav_nchan.eq.1) then c file_out=file_out(1:seiclen(file_out)) //'_'// c1 & wav_comp(1) c endif c c check if file exist, if so add SEI to name, jh mar 06 c inquire(file=file,exist=exist) if(exist) then file_out=file_out(1:seiclen(file_out))//'_SEI' endif do i=seiclen(file_out)-3,seiclen(file_out) if (file_out(i:i).eq.' ') file_out(i:i)='_' enddo c c if a correction is made, write in main header c if(ground.or.flow.ne.0.0.or.fhigh.ne.0.0) then mainhead(1)(2:13)='Filtered raw' if(ground_type.eq.1) mainhead(1)(2:13)='Displacement' if(ground_type.eq.2) mainhead(1)(2:13)='Velocity ' if(ground_type.eq.3) mainhead(1)(2:13)='Acceleration' if(flow.lt.10.0) write(mainhead(1)(15:19),'(f5.3)') flow if(flow.ge.10.0) write(mainhead(1)(15:19),'(f5.1)') flow if(fhigh.lt.10.0.and.fhigh.ge.0.0) * write(mainhead(1)(21:25),'(f5.3)') fhigh if(fhigh.ge.10.0) write(mainhead(1)(21:25),'(f5.1)') fhigh if(flow.ne.0.0.or.fhigh.ne.0.0) then mainhead(1)(27:28)='hz' mainhead(1)(20:20)='-' endif endif do i=1,8 write(6,'(1x,a)') mainhead(i)(1:78) enddo open(96,file=file_out,status='unknown',form='unformatted') do k=1,12 write(96)mainhead(k) enddo if(wav_out_nchan.gt.30) then k=(wav_out_nchan-31)/3+1 do i=13,k+12 write(96)mainhead(i) enddo endif endif c c---------------------------------------------------------- c GSE format has no main header c---------------------------------------------------------- c if(format_out(1:3).eq.'GSE'.or. & format_out(1:4).eq.'SEED') then c delete output file if (file_out(1:4).eq.'SEED') then seed_file_out=' ' seed_file_out=file(1:seiclen(file))//'_SEED' endif if (file_out(1:3).eq.'GSE'.or. & file_out(1:4).eq.'SEED') then file_out=file(1:seiclen(file))//'_GSE' endif open(96,file=file_out,status='unknown') close(96) open(96,file=file_out,status='unknown') endif c c---------------------------------------------------------- c SAC nothing done about header c---------------------------------------------------------- c c c----------------------------------------------------------- c MSEED, just open file c----------------------------------------------------------- c if(format_out(1:5).eq.'MSEED') then if (file_out(1:5).eq.'MSEED') then file_out=file(1:seiclen(file))//'_MSEED' endif op_file=open_file(96,file_out,12) ! use block size 4096=2**12 if(.not.op_file) then write(6,*) 'Error opening miniseed output file' stop endif endif c c---------------------------------------------------------- c Channel processing and write out loop c----------------------------------------------------------- c mseed_block=1 do i=1,wav_out_nchan k=wav_out_chan(i) c read whole channel c call wav_read_channel(k) c c correct number of samples, reading process might have found a gap c commnet out dec 19 2007 jh: this was wrong since the the number c of samples in whole interval was was used, not the number in c selected interval, seemed to only affect seed since this number was c reset after reading c c nsamp(k)=wav_out_nsamp(k) cjh c write(6,*) c * 'nsamp,out,wav',nsamp(k),wav_out_nsamp(k),wav_nsamp(k) c write(6,*)' Samples after correction',nsamp(k) c c filter but only if ground motion not done in which case c filtering is done in spectral domain c if((flow.ne.0.0.or.fhigh.ne.0.0).and..not.ground) then c c check sample rate before filtering c if(fhigh.ge.wav_rate(k)/2.0) then write(6,'(a,f8.2)')' Sample rate is: ',wav_rate(k) write(6,'(a,f8.2)') ' Filter is: ',fhigh write(6,'(a)')' Cannot filter **********' else c c remove DC c ndc=0 call remove_dc(signal1,ndc,dc,wav_nsamp(k)) c write(6,*) ' filter ',flow,'-',fhigh call bndpas(flow,fhigh,1000.0/wav_rate(k), * cof,gain) call filter(signal1,wav_nsamp(k),cof,gain,1) ! 4 pole one way endif endif c c compute ground motion c if (ground) then wav_current_chan(1) = k wav_resp_file = ' ' call read_resp if (wav_resp_status.ne.'9') then call remove_resp(signal1,y_com,wav_nsamp(k), & wav_rate(k),ground_type, & flow,fhigh,8,8) endif endif c c compute rsam data c if (rsam_flag) then call remove_dc(signal1,ndc,dc,wav_nsamp(k)) call signal_to_rsam(k,2) endif c c reduce signal to what is seen c if(max_points.gt.0) then call signal_to_visual(k,max_points,2) endif c c cut out data and put data in beginning of array c j=0 do l=wav_out_first_sample(k),wav_out_nsamp(k)+ * wav_out_first_sample(k)-1 j=j+1 signal1(j)=signal1(l) enddo c c skip points c xxx c j=0 c if(nskip(k).gt.1) then c do l=1,wav_out_nsamp(k),nskip(k) c j=j+1 c signal1(j)=signal1(l) c enddo c c fix sample rate and number of samples according to skipping c c wav_out_rate(k)=wav_out_rate(k)/nskip(k) c wav_out_nsamp(k)=wav_out_nsamp(k)/nskip(k) c endif c c write out, ddifferent formats c c c----------------- c SEISAN format c----------------- c if(format_out.eq.'SEISAN') then call sheads(year,month,day,hour,min,sec,wav_out_nchan,i, * 'NEWAV',h_text,stat,comp, * nsamp,rate,cbyte, * file,mainhead,chahead) c c put in location and network c chahead(8:8)=location(i)(1:1) chahead(13:13)=location(i)(2:2) chahead(17:17)=network(i)(1:1) chahead(20:20)=network(i)(2:2) write(6,'(1x,a)') chahead(1:79) c c if filtered or correction for ground response, values might be very small so c numbers are multiplied by a gain factor and the gain factor is put in header c if(ground.or.flow.ne.0.0.or.fhigh.ne.0.0) then c c find largest value to be sure data can be scaled up or down c ymax=0.0 do l=1,nsamp(i) if(abs(signal1(l)).gt.ymax) ymax=abs(signal1(l)) enddo c c find factor c factor=ymax/1.0e9 do l=1,nsamp(i) signal1(l)=signal1(l)/factor enddo chahead(76:76)='G' write(chahead(148:159),'(g12.7)') factor write(chahead(148:159),'(g12.7)') factor ! save factor in header endif write(96) chahead c cjh apr 2003 if(cbyte(i).eq.' '.or.cbyte(k).eq.'2') then cjh do l=1,nsamp(i) cjh signal2b(l)=signal1(l) cjh enddo cjh write(96)(signal2b(l),l=1,nsamp(i)) cjh else do l=1,nsamp(i) signal4b(l)=signal1(l) enddo write(96)(signal4b(l),l=1,nsamp(i)) cjh endif endif c c--------------------------------------- c GSE format c--------------------------------------- c if(format_out(1:3).eq.'GSE'.or. & format_out(1:4).eq.'SEED') then c c make sure number is not too big c do l=1,nsamp(i) if(signal1(l).gt. 100000000.0) signal1(l)= 100000000.0 if(signal1(l).lt. -100000000.0) signal1(l)=-100000000.0 enddo if (seisweb) then do l=1,nsamp(i) signal1(l)=signal1(l)*10. enddo endif if(format_out(1:5).eq.'SEED ') then call write_seisan_to_gse(96,k,.true.,'GSE2SEED') elseif(format_out(1:6).eq.'SEEDDL') then call write_seisan_to_gse(96,k,.true.,'GSE2SEEDDL') else call write_seisan_to_gse(96,k,.true.,format_out) endif c c--------------------------------------- c SAC format c--------------------------------------- c elseif(format_out(1:3).eq.'SAC') then call write_seisan_to_sacbin(k,.false., * .true.,evla,evlo,evel,nerr) c c----------------------------------------- c MINISEED format c------------------------------------------ c elseif(format_out(1:5).eq.'MSEED') then mseed_year=year(i) mseed_month=month(i) mseed_day=day(i) mseed_hour=hour(i) mseed_minute=min(i) mseed_second=sec(i) mseed_station=stat(i) mseed_channel(1:2)=comp(i)(1:2) mseed_channel(3:3)=comp(i)(4:4) mseed_total_samples=nsamp(i) mseed_sample_rate=rate(i) c mseed_location=seed_location_code mseed_location=location(i) mseed_network=network(i) mseed_encoding=steim1 ! steim1 do l=1,nsamp(i) signal4b(l)=signal1(l) enddo call seed_write_chn(96,signal4b,mseed_block) endif enddo c c close output waveform file c if (format_out(1:6).eq.'SEISAN'.or.format_out(1:3).eq.'GSE'. * or.format_out(1:5).eq.'MSEED'.or.format_out(1:4).eq.'SEED') *close(96) if (format_out(1:4).eq.'SEED') then system_call = ' ' system_call = 'gse2seed -i '//file_out(1:seiclen(file_out)) & //' -o '//seed_file_out(1:seiclen(seed_file_out)) write(*,*) system_call(1:seiclen(system_call)) call systemc(system_call,seiclen(system_call)) endif c----------------------------------------------------- goto 500 c 999 continue write(6,*)' Wavetool terminated: error' if(seisweb) then write(6,*)'Wavetool terminated: error' endif goto 600 c 500 continue if(seisweb) then write(6,*)' Waveform extraction terminated: ok' endif text=' ' text(1:2)='OK' text(3:62)=file_out(1:60) open(1,file='extract.mes',status='unknown') write(1,'(a)') text close(1) c call put_seisan_message(text) if (format_out(1:4).eq.'SEED') then write(6,'(1x,a,a)')' Output waveform file name is ', & seed_file_out else write(6,'(1x,a,a)')' Output waveform file name is ',file_out endif if(in.eq.1) goto 11 600 continue stop end subroutine signal_to_visual(ichan,max_points,option) c c reduce number of samples in single trace, so that total windows of all c traces will have max_points c c input: ichan - channel index for wav_out c' max_points - number of points for total window c option - 1: header only, 2: header and data c implicit none save include 'seidim.inc' include 'waveform.inc' integer ichan ! channel index integer max_points ! number of points in total time window integer option ! option real data(max_sample) ! data array real min,max ! min and max in interval delta double precision delta ! size of time window to take data double precision stop,now integer nsamp_save(max_trace) real rate_save(max_trace) integer i,k integer ind(2) ! index to sample of min and max c c option 1, only set wav_out_nsamp and wav_out_rate c if (option.eq.1) then c c get length of total output window c call wav_out_index_total delta = wav_out_total_time / max_points c write(*,*) ichan,wav_out_rate(ichan), c * wav_out_nsamp(ichan),wav_out_total_time nsamp_save(ichan)=wav_out_nsamp(ichan) rate_save(ichan)=wav_out_rate(ichan) wav_out_rate(ichan)=1./delta wav_out_nsamp(ichan)=wav_out_duration(ichan)/delta c write(*,*) ichan,wav_out_rate(ichan), c * wav_out_nsamp(ichan) return elseif (option.eq.2) then c c reset samples and rate c wav_out_rate(ichan)=rate_save(ichan) wav_out_nsamp(ichan)=nsamp_save(ichan) c c get length of total output window c delta = wav_out_total_time / max_points c c go out if max_points larger than number of samples c if (max_points.ge.wav_out_nsamp(ichan)) return c c get time interval for new sampling c stop = wav_abs_time(wav_first)+wav_out_start(ichan)+delta k=1 min=signal1(wav_out_first_sample(ichan)) ind(1)=wav_out_first_sample(ichan) max=signal1(wav_out_first_sample(ichan)) ind(2)=wav_out_first_sample(ichan) c write(17,*) 'max,min ',max,min c c loop through all samples c do i=wav_out_first_sample(ichan),wav_out_first_sample(ichan)+ * wav_out_nsamp(ichan)-1 c c time of sample c now=wav_abs_time(wav_first)+wav_out_start(ichan)+ * (i-wav_out_first_sample(ichan))/wav_out_rate(ichan) c write(17,*) now,stop,ichan,wav_out_first_sample(ichan), c * signal1(i) if (now.le.stop) then ! check if sample in interval if (signal1(i).ge.max) then max=signal1(i) ind(2)=i endif if (signal1(i).le.min) then min=signal1(i) ind(1)=i endif c c sort with time c if (ind(2).le.ind(1)) then data(k)=max data(k+1)=min else data(k+1)=max data(k)=min endif else ! new interval c write(17,*) k,data(k),data(k+1) k=k+2 stop=stop+2*delta max=signal1(i) min=signal1(i) data(k)=max data(k+1)=min endif enddo c c set number of samples and first sample in wav_out c wav_out_rate(ichan)=1./delta wav_out_nsamp(ichan)=wav_out_duration(ichan)/delta wav_out_first_sample(ichan)=1 c c write data to signal1 c do i=1,wav_out_nsamp(ichan) signal1(i)=data(i) enddo c c change number of samples and sample rate c endif return end subroutine get_argument_from_file(arguments,nargs) implicit none include 'libsei.inc' character*80 arguments(*) integer nargs integer read1 integer seiclen character*80 strline integer code logical b_flag integer i call sei open( old$+warn$, ! Open file & warn of erro. & ' ', ! Prompt (n/a). & arguments(2), ! Filename. & read1, ! Unit opened on. & b_flag, ! Flag. & code ) ! Returned condition. if( code .ne. e_ok$ ) then ! An error or does not exist. write(6,'(a,a)') ' File does not exist: ',arguments(2) stop ! stop if input from choise file ! endif do i=1,9999 read(read1,'(a)',end=9711) strline arguments(i) = strline(1:seiclen(strline)) enddo 9711 nargs = i - 1 close(read1) return end subroutine signal_to_rsam(ichan,option) c c compute 1 s p minute rsam data c c input: ichan - channel index for wav_out implicit none save include 'seidim.inc' include 'waveform.inc' integer ichan ! channel index real data(max_sample) ! data array double precision delta ! size of time window to take data double precision msec1,msec2,start integer i,k integer year,month,day,hour,min,doy real sec,secprev real rsam integer nrsam integer option if (option.eq.1) then msec1=wav_abs_time(ichan) msec2=wav_abs_time(ichan)+wav_nsamp(ichan)/wav_rate(ichan) k=int((msec2-msec1)/60.) call sectim(msec1,year,doy,month,day,hour,min,sec) if (sec.ne.0.) k=k+1 call sectim(msec2,year,doy,month,day,hour,min,sec) if (sec.ne.0.) k=k+1 wav_out_rate(ichan)=1./60. wav_out_nsamp(ichan)=k wav_out_first_sample(ichan)=1 else c c compute rsam data c secprev=-1 rsam=0. k=0 do i=1,wav_nsamp(ichan) c c get time of sample c msec1=wav_abs_time(ichan)+(i-1)/wav_rate(ichan) call sectim(msec1,year,doy,month,day,hour,min,sec) if (i.eq.1) then wav_out_sec(ichan)=0. endif if (sec.lt.secprev) then k=k+1 if (nrsam.ne.0) then data(k)=rsam/float(nrsam) else data(k)=0 endif nrsam=1 rsam=abs(signal1(i)) secprev=-1 else secprev=sec nrsam=nrsam+1 rsam=rsam+abs(signal1(i)) endif enddo k=k+1 data(k)=rsam/float(nrsam) c c write data to signal1 c do i=1,wav_out_nsamp(ichan) signal1(i)=data(i) enddo endif return end