c Checks for 'CHK2' changed to checks for 'CHK2 ', because it occurred that c compressed waveform data started with 'CHK2...' on the first column of a c line. Urs Kradolfer, 3 May 2000 c c Set hdr_instr = '_' (instead of ' ') in subr. init_hdrvars. 3Dec99 K.Koch c c Output-Lines in GSE-Format have always a length of 80 characters, also the c last line of the DAT2 section. Urs Kradolfer, 6.6.99 c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::== C MODIFICATION FOR SESAME HV CODE by a.tento 04-05/02: c === write(6, --> write(*, c === problem read(a(1:1),*) bla c === GSEIN : if (hdr_nsamp .gt. c_sigsize) c ==> exit routine without readind data and with ierr = 0 c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::== c c Version 98.002 - 27.11.98 c============================================================================== c Karl Koch, 08.05.98 c c============================================================================== c c----- subroutine gsein(iin,iout,cbuf,iy,ichecksum,ierr) c------------------------------------------------------------------------------ c Input/Output parameters: c iin ......... logical unit with input data c iout ........ logical unit for output data c cbuf ........ character array to store compressed data (see c_bufsize) c iy .......... integer array to hold raw waveform data (see c_sigsize) c ichecksum ... checksum value for data stored in iy according to format c ierr ........ error code c c in common block c hdr_ifmtnum . GSE format identifier of input data c hdr_nsamp ... number of samples in array iy c c disabled c iformat_in .. GSE format identifier (ie. =2 for GSE2.0) (output) c nmax ........ number of samples in array iy (output) c lin ......... character array to hold header lines c ilin ........ number of header lines (output) c------------------------------------------------------------------------------ c implicit none c -- include header variables (hdr_...) and constants (c_...) include 'codeco_common.f' c -- parameters integer*4 iin, iout character*1 cbuf(1) integer*4 iy(1) integer*4 ichecksum, ierr c -- local variables integer*4 icomp, i, j, ierror, nchecksum integer*4 numout, nchar, ii integer*4 trimlen integer*4 icheck_sum_gse1,icheck_sum_gse2 real fsec ! seconds and milliseconds character line*132, fmt*6 character firstline*132 ! store first line for error message character buff*60 c -- for month conversions: integer mnths(12) data mnths/31,28,31,30,31,30,31,31,30,31,30,31/ c -- executable code hdr_coosys = ' ' if (hdr_debug .gt. 2) write (*,*) '-- entering gsein' ierr = 0 c -- loop all lines in GSE file until WID is found 1 read(iin,'(a)',end=999) line if (line(1:4) .ne. 'WID1' .and. line(1:4) .ne. 'WID2') then c -- copy all lines before WID-line if (hdr_ofmtnum .gt. 0) & write(iout,'(a)') line(1:trimlen(line)) goto 1 else if (line(1:4).eq.'WID1') hdr_ifmtnum = 1 if (line(1:4).eq.'WID2') hdr_ifmtnum = 2 endif firstline = line if (hdr_debug .gt. 3) then write (*,*) '-- found hdr line: ' write (*,'(a)') firstline endif c if (hdr_ifmtnum .eq. 2) goto 20 c c -- input data are in GSE1.0 format c if (hdr_debug .gt. 0) write(*,*) & '-- reading GSE1.0 format data' c -- parse first header line of GSE1.0 100 format(i4,i3,3(1x,i2),1x,i3,1x,i8,1x,a6,1x,a8,1x, & a2,1x,f11.7,1x,a6,1x,a4,1x,i1) read( line(7:), 100 ) hdr_year, hdr_jday, hdr_hour, hdr_min, & hdr_sec, hdr_msec, hdr_nsamp, hdr_station, hdr_stadescr, & hdr_chan, hdr_smprate, hdr_instr, hdr_ifmtname, hdr_idiff c -- read and parse second header line read( iin, '(a)' ) line 101 format(f9.6,i1,f7.4,1x,4(f9.4,1x),2(f7.2,1x),f6.1,1x,f6.1) read( line, 101 ) hdr_calfac, hdr_calunit, hdr_calper, & hdr_stalat, hdr_stalon, hdr_staelev, hdr_stadepth, & hdr_beamaz, hdr_beamslo, hdr_hang if (hdr_nsamp .gt. c_sigsize) goto 998 c -- compute day and month from julian day cks if(ifix(float(hdr_year)/4.0)*4 .eq. hdr_year) mnths(2)=29 cks hdr_month = 0 cks do i=1,12 cks hdr_month = hdr_month + mnths(i) cks if (hdr_month .ge. hdr_jday) goto 60 cks enddo cks60 hdr_month = hdr_month - mnths(i) cks hdr_day = hdr_jday - hdr_month cks hdr_month = i call julien( hdr_day, hdr_month, hdr_year, hdr_jday, 0 ) c -- adjust channel name (insert a 'H' in the middle) hdr_chan(3:3) = hdr_chan(2:2) hdr_chan(2:2) = 'H' call casefold( hdr_chan ) c -- read away all lines until DAT1 if (hdr_debug .gt. 3) write (*,*) '-- find DAT1' 2 read( iin, '(a)' ) line if (line(1:4) .ne. 'DAT1') goto 2 c -- Now read in data according to the format: if (hdr_debug .gt. 2) write (*,*) '-- read in samples' if (hdr_ifmtname .eq. 'INTV') then if (hdr_debug .gt. 3) write (*,*) '-- read INTV data' read(iin,*) (iy(i),i=1,hdr_nsamp) read(iin,'(4x,6x,i15)') ichecksum if (hdr_debug .gt. 3) write (*,*) '-- done' endif if(hdr_ifmtname(1:3).eq.'INT'.and.hdr_ifmtname.ne.'INTV')then if (hdr_debug .gt. 3) write (*,*) '-- read INTx data' buff=hdr_ifmtname(4:4) read(buff,*) ii read(iin,'(a)') line i=(trimlen(line)+ii-1)/ii write(fmt,'(''('',i2.2,''i'',i1.1,'')'')') i,ii backspace(3) read(iin,fmt) (iy(i),i=1,hdr_nsamp) read(iin,'(4x,6x,i15)') ichecksum if (hdr_debug .gt. 3) write (*,*) '-- done' endif c if(hdr_ifmtname(1:3).eq.'CMP')then if (hdr_debug .gt. 3) write (*,*) '-- read CMPx data' buff=hdr_ifmtname(4:4) read(buff,*) icomp j=0 c -- loop all data lines 10 continue read(iin,'(a)') line if(line(1:4).eq.'CHK1') goto 11 nchar=trimlen(line) do i=1,nchar j=j+1 cbuf(j)=line(i:i) enddo goto 10 11 continue read(line(5:),*) ichecksum if (hdr_debug .gt. 3) write (*,*) '-- done' nchar=j numout=c_sigsize ! on input, give max. size of array iout (=iy) if (hdr_debug .gt. 3) write (*,*) '-- decompress', icomp if(icomp.eq.6) call DCOMP6(nchar,cbuf,numout,iy,ierror) if(icomp.eq.7) call DCOMP7(nchar,cbuf,numout,iy,ierror) if(icomp.eq.8) call DCOMP8(nchar,cbuf,numout,iy,ierror) if(ierror.eq.-1) write(*,*)'Integer*4 array not large enough !' if(numout.ne.hdr_nsamp)then write(*,*)'Error in retrieving correct number of samples !' write(*,*)'Expected ',hdr_nsamp,', received ',numout endif if (hdr_debug .gt. 3) write (*,*) '-- done' endif c -- remove first or second differences if necessary if (hdr_debug .gt. 3) write (*,*) '-- remove differences diff:', & hdr_idiff, 'no. samples: ', hdr_nsamp do i=1,hdr_idiff call REMDIF1(iy,hdr_nsamp) enddo if (hdr_debug .gt. 3) write (*,*) '-- done' c if (hdr_ifmtnum.eq.1) goto 50 write(*,*) 'There is good news and bad news at this point:' write(*,*) 'First the good news: no more CPU time will be used!!!' write(*,*) 'and now the bad news: You have screwed up !!!' stop c c -- The input data are GSE2.0 c 20 if (hdr_debug .gt. 0) write(*,*) '-- reading GSE2.0 format data' c -- parse WID2 line 200 format(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f6.3,1x,a5,1x,a3,1x, & a4,1x,a3,1x,i8,1x,f11.6,1x,e10.2,1x,f7.3,1x,a6,1x,f5.1,1x,f4.1) read( line(6:), 200 ) hdr_year, hdr_month, hdr_day, hdr_hour, & hdr_min, fsec, hdr_station, hdr_chan, hdr_stadescr, & hdr_ifmtname, hdr_nsamp, hdr_smprate, hdr_calfac, hdr_calper, & hdr_instr, hdr_hang, hdr_vang hdr_sec = ifix( fsec ) hdr_msec = nint( (fsec-float(hdr_sec)) * 1000.0 ) if (hdr_debug .gt. 3) write (*,*) '-- got start time: ', & hdr_year, hdr_month, hdr_day, hdr_hour, hdr_min, hdr_sec, & hdr_msec if (hdr_nsamp .gt. c_sigsize) goto 998 c -- compute julian day cks if (ifix(float(hdr_year)/4.0)*4 .eq. hdr_year) mnths(2)=29 cks hdr_jday = 0 cks do i=1,hdr_month-1 cks hdr_jday = hdr_jday + mnths(i) cks enddo cks hdr_jday = hdr_jday + hdr_day call julien( hdr_day, hdr_month, hdr_year, hdr_jday, 1 ) c -- read away all lines until DAT2, read possibly existing STA2-line if (hdr_debug .gt. 3) write (*,*) '-- find DAT2' 22 read( iin, '(a)' ) line if (line(1:4) .eq. 'STA2') then read( line, & '(5x,a9,1x,f9.5,1x,f10.5,1x,a12,1x,f5.3,1x,f5.3)' ) & hdr_network, hdr_stalat, hdr_stalon, hdr_coosys, & hdr_staelev, hdr_stadepth hdr_staelev = hdr_staelev * 1000.0 hdr_stadepth = hdr_stadepth * 1000.0 if (hdr_debug .gt. 3) write (*,*) '-- found STA2 ', & hdr_network, hdr_stalat, hdr_stalon, hdr_coosys, & hdr_staelev, hdr_stadepth endif if (line(1:4) .ne. 'DAT2') goto 22 if (hdr_debug .gt. 3) write (*,*) '-- done' c -- Now read in data according to the format: if (hdr_debug .gt. 2) write (*,*) '-- reading samples' if (hdr_ifmtname .eq. 'INT') then if (hdr_debug .gt. 3) write (*,*) '-- reading INT data' read(iin,*) (iy(i),i=1,hdr_nsamp) 30 read(iin,'(a)') line if(line(1:5).eq.'CHK2 ') goto 31 goto 30 31 read(line(5:),*) ichecksum if (hdr_debug .gt. 3) write (*,*) '-- done' endif c if (hdr_ifmtname(1:2) .eq. 'CM') then buff=hdr_ifmtname(3:3) read(buff,*) icomp if (hdr_debug .gt. 3) write (*,*) '-- reading CM', icomp j=0 c -- reading all compressed lines 32 continue read(iin,'(a)') line if(line(1:5).eq.'CHK2 ') goto 34 nchar=trimlen(line) do i=1,nchar j=j+1 cbuf(j)=line(i:i) enddo goto 32 34 continue read(line(5:),*) ichecksum if (hdr_debug .gt. 3) write (*,*) '-- done' nchar=j numout=c_sigsize ! on input, give max. size of array iout (=iy) if (hdr_debug .gt. 3) write (*,*) '-- decompress data' if(icomp.eq.6) call DCOMP6(nchar,cbuf,numout,iy,ierror) if(icomp.eq.7) call DCOMP7(nchar,cbuf,numout,iy,ierror) if(icomp.eq.8) call DCOMP8(nchar,cbuf,numout,iy,ierror) if(ierror.eq.-1) write(*,*)'Integer*4 array not large enough !' if (hdr_debug .gt. 3) write (*,*) '-- done' if(numout.ne.hdr_nsamp)then write(*,*)'Error in retrieving correct number of samples !' write(*,*)'Expected ',hdr_nsamp,', received ',numout endif c -- remove first and second differences if (hdr_debug .gt. 3) write (*,*) '-- remove 2 differences' call REMDIF1(iy,hdr_nsamp) call REMDIF1(iy,hdr_nsamp) if (hdr_debug .gt. 3) write (*,*) '-- done' endif c 50 continue if (hdr_debug .gt. 3) write (*,*) '-- check checksums' if (hdr_ifmtnum .eq. 1) then c -- calculate checksum for GSE1.0 format nchecksum=icheck_sum_gse1(iy,hdr_nsamp) elseif (hdr_ifmtnum .eq. 2) then c -- calculate checksum for GSE2.0 format nchecksum=icheck_sum_gse2(iy,hdr_nsamp) else print *, 'Unknown GSE-Format ',hdr_ifmtnum,'.x' stop endif if(nchecksum.ne.ichecksum)then write(*,*)'Error in checksum in the following waveform' write(*,'(a)') firstline(1:trimlen(firstline)) write(*,*)'Read checksum =',ichecksum write(*,*)'Calculated checksum =',nchecksum write(*,*) ichecksum=nchecksum ! ... use the correct checksum !!! endif if (hdr_debug .gt. 2) write (*,*) '-- leaving gsein' return c c998 print *,' >>>> WARNING: Program dimensions too small for ', c & ' the following waveform :' c print *, line c goto 1 998 return 999 ierr = 999 return end