C+ program st2sac C C Converts a st.XXX file (demultiplexed suds) to a set of sac C format files. C C 22 August 95 jas/vtso C- parameter (maxN=131074) real*4 signal(maxN) character contents*512,buffer*512,ymdhms*21,getstring*80 character filename*80,filen*10,net_sta(32)*4 character file_now*80,com*8/'TIME OK!'/ integer*2 nstruct,const1,const2,threshold,const3,const4 integer*2 sta(32),lta(32),abs_sta(32),abs_lta(32) integer*2 trig_value(32),num_triggers(32),i2block(256) real*8 beginttime,trig_time,real_time character net(32)*3,netwname*3,network*3,st_name(32)*4 character blank*1/' '/,t_sta_nam(32)*4,st_dummy*4 character null*1/0/,kevnm*16,component*1 integer*2 actual,trgnet,trval,decim,tstmin,tkeep integer*2 tavrag,toteta,enum,edenom,econst,etastn,knet integer*4 hour logical*1 polflp/.true./ data ymdhms(10:10)/' '/ C open(unit=2,name='st2sac.lst',status='unknown') CALL TIMDAT(2,'st2sac') nch = 0 open(unit=3,file='scratch.dat', 1 status='scratch',recl=512,access='direct') filename = getstring('Input st filespec') open(unit=1,file=filename,status='old', 1 form='unformatted',recl=512,access='direct') kevnm = getstring('Comment (up to 16 characters)') write(2,'(''Input file '',a)') filename(1:lenc(filename)) write(2,'(1x,a)') kevnm C C Get the TRIGSETTING structure (number 26) C nleft = 0 kb = 0 call get_suds_structure(1,contents,nhead,ndata, 1 buffer,nstruct,nleft,kb) read(contents(1:nhead-4),'(A3,1x,A8,5a2,2x,2a4)') netwname, 1 beginttime,const1,const2,threshold,const3,const4, 2 sweep,aperature ct write(2,*) nstruct,nhead,ndata,blank,netwname, ct 1 beginttime,const1,const2,threshold,const3,const4, ct 2 sweep,aperature,nleft C C Now for the rest of the prelims C m = 0 ntrig = 0 now = 0 nch = 0 do while (nleft .ge. 0) call get_suds_structure(1,contents,nhead,ndata, 1 buffer,nstruct,nleft,kb) ct write(2,*) 'In loop',nch, nstruct,nhead,ndata,nleft,kb if (nleft .ge. 0) then if (nstruct .eq. 25) then m = m + 1 read(contents(1:nhead),'(A3,1x,A4,1x,a1,2x,6a2,a8)') 1 net(m),net_sta(m),component,sta(m),lta(m),abs_sta(m), 2 abs_lta(m),trig_value(m),num_triggers(m),trig_time if (net_sta(m)(4:4) .eq. null) net_sta(m)(4:4) = blank if (net(m) .eq. netwname) then ntrig = ntrig + 1 t_sta_nam(ntrig) = net_sta(m) end if ct write(2,*) blank,network,blank,st_name, ct 1 blank,component,sta,lta,abs_sta,abs_lta,trig_value, ct 2 num_triggers,trig_time C else if (nstruct .eq. 5) then read(contents(1:10),'(A3,1x,A4,1x,a1)') network, 1 st_dummy,component if (nch .gt. 3 .and. component .eq. 'x') then nleft = -1 else nch = nch + 1 ct write(2,*) 'nstruct is 5',nch,blank,network, ct 1 blank,st_name,blank,component,blank,nleft end if C C nstruct=5 is stationcomp. I do not think there is anything we C want to keep from it. For the record, we probably want to tell C the system that the polarities are "r" rather then "n". I do C not think anything gets flipped automatically if you say it, though. C The station identifier is in the next one, so even that is not C needed here. C else if (nstruct .eq. 7) then read(contents(1:36),'(4x,A4,4X,a8,12X,a4)') st_name(nch), 1 real_time,rate_now if (st_name(nch)(4:4) .eq. null) st_name(nch)(4:4) = blank if (nch .eq. 1) then ipts4 = ndata/2 rate = rate_now call unixtime(real_time,jyear,jul,hour,min,sec) ct write(2,*) jyear,jul,hour,min,sec,' unixtime' call julian(jul,jyear,jmo,jday,ymdhms,.true.,.false.) write(ymdhms(11:21),'(i2.2,'':'',i2.2,'':'',F5.2)') 1 hour,min,sec if (ymdhms(17:17) .eq. ' ') ymdhms(17:17) = '0' ct write (2,'(1x,a21)') ymdhms C end if nblocks = ndata/512 nremain = ndata - 512*nblocks if (nremain .gt. 0) nblocks = nblocks + 1 nleft_in = nleft/2 nleft_out = nleft_in - nremain/2 if (nleft_out .lt. 0) nleft_out = nleft_out + 256 call move_suds_data(1,3,nleft_in,nleft_out,buffer, 1 nblocks,now,nch,kb) ct if (nch .lt. 3) write(2,*) 'Out of move_suds: now,kb',now,kb nleft = 2*nleft_out C C nstruct=7 is descriptrace. This has both a header and adata section. C This routine reads from unit 1 and writes to unit 3. C end if end if end do c c header time c actual = nch pretrg = 20.0 if (netwname .eq. 'gco' .or. netwname .eq. 'GCO') then trgnet = 1 else trgnet = 2 end if trval = 24 ! Guess decim = 3 tstmin = nint(aperature) tkeep = 60 tavrag = sweep toteta = 0 knet = ntrig etastn = 0 do k=1,m toteta = toteta + num_triggers(k) if (num_triggers(k) .gt. 0 .and. trig_value(k) .gt. 0) then etastn = etastn + 1 end if end do enum = const1 edenom = const2 econst = threshold write(2,'(''Initial: '',a21,10x,i2,'' channels''/)') ymdhms,nch call piktim(3,2,rate,ipts4,ymdhms,ncom,com,jul) tottim = (ipts4-1)/rate write(filen,'(a2,i3.3,2a2,''.'')') ymdhms(8:9),jul, 1 ymdhms(11:12),ymdhms(14:15) c write(2,'(''Stations:'',T13,13(a4,1x)/T13,13(a4,1x)/T13, 1 13(a4,1x))') (st_name(jj),jj=1,nch) write(2,2) YMDHMS,PRETRG write(2,3) RATE,IPTS4,TOTTIM,TOTETA,ETASTN write(2,4) TRGNET,(T_sta_nam(JJ),JJ=1,KNET) write(2,'(/,''net station sta Lta a_sta a_lta'', 1 '' eta num_trig'')') do k=1,m write(2,'(1x,a3,3x,a4,i8,i6,i8,i8,i8,2x,i5)') 1 net(k),net_sta(k),sta(k),lta(k),abs_sta(k),abs_lta(k), 2 trig_value(k),num_triggers(k) end do write(2,11) write(2,*) com write(2,*) nblocks,' blocks ',trval,' is trigger reset ', 1 decim,' is decimation factor' write(2,*) tstmin,' is minimum trigger time ',tkeep, 1 ' is time after last trigger' write(2,*) tavrag,' is TAVRAG ',enum,edenom,econst, 1 ' are num, denom and threshold' 2 FORMAT(/,'First point: ',a21,' Pretrigger time: ', . F5.2) 3 FORMAT('Digitizing rate: ',F6.2,' Points/channel: ',I5, . ' Total time: ',F6.2/'TOT Trigs =',I6,' ETASTN =',I3) 4 FORMAT('Trigger subnet # ',I1,T21,10(A4,1X)) 11 FORMAT('Polarity of data channels reversed to recover true', & ' polarity of ground motion.'/) read(ymdhms(11:12),*) hour read(ymdhms(14:15),*) minute read(ymdhms(17:18),*) nzsec read(ymdhms(20:21),*) nzmsec nzmsec = 10*nzmsec dt = 1.0/rate secs = 0.0 nfile = lenc(filename) filename = filename(1:nfile-4)//'.' nfile = nfile - 3 C call newhdr call setfhv('B',secs,nerr) call setnhv('NZYEAR',jyear,nerr) call setnhv('NZJDAY',jul,nerr) call setnhv('NZHOUR',hour,nerr) call setnhv('NZMIN',minute,nerr) call setnhv('NZSEC',nzsec,nerr) call setnhv('NZMSEC',nzmsec,nerr) call setkhv('knetwk','DIGSYS',nerr) call setkhv('kevnm',kevnm,nerr) call setnhv('NPTS',ipts4,nerr) call setfhv('DELTA',dt,nerr) k = 0 write(2,*) ' ' do j=1,actual jj = 0 do i=1,nblocks k = k + 1 read(3'k) i2block do kk=1,256 signal(jj+kk) = i2block(kk) end do jj = jj + 256 end do call setkhv('KSTNM',st_name(j),nerr) if (j .gt. 1) then if (st_name(j)(4:4) .eq. 'N') then call setkhv('KCMPNM','N',nerr) call setfhv('cmpinc',90.0,nerr) call setfhv('cmpaz',0.0,nerr) else if (st_name(j)(4:4) .eq. 'E') then call setkhv('KCMPNM','E',nerr) call setfhv('cmpinc',90.0,nerr) call setfhv('cmpaz',90.0,nerr) else call setkhv('KCMPNM','Z',nerr) call setfhv('cmpaz',-12345.,nerr) call setfhv('cmpinc',0.0,nerr) end if end if call maxmin(signal,ipts4,smin,smax,'R4',2,0) if (smax-smin .lt. 10) then write(2,*) 'No sac file written for channel ',st_name(j) else file_now = filen//st_name(j) call wsac0(file_now,dummy,signal,nerr) Write(2,*) ' Output to file ',file_now(1:nfile+4) end if end do call exit(0) end C+ subroutine piktim(Ludata,LUout,rate,ipts4,ymdhms,ncom,com,day) C subroutine piktim(Ludata,LUout,rate,ipts4,ymdhms,ncom,com,day) C C Adapted from J. Bollinger's PDP program to correct the rate and C the starting time. LUDATA is where the time channel lives C starting in block 2 for ipts4 points. RATE is 1/digitizing C interval, and YMDHMS is a 22 byte array with date and time. C Both RATE and YMDHMS may be changed for output. Things are C written to LUOUT. If there is an error, ncom is reset to 30 and C the last 8 characters is COM become "TIME BAD". C C JAB's deathless prose follows. C C James A. Bollinger, May 21,1985 C- REAL*4 TOTTIM,RATE,SEC,AVG INTEGER*2 clock(12800) character*21 ymdhms,ymdetc character*8 com INTEGER MIN,HRS,DAY LOGICAL*1 ERROR C npts = min0(ipts4,12800) nblocks = npts/256 if (256*nblocks .lt. npts) nblocks = nblocks + 1 kend = 0 do j=1,nblocks kst = kend + 1 kend = kend + 256 read (ludata'j) (clock(jj),jj=kst,kend) end do C C Call the subroutine that reads the time from the time signal. C CALL TIMSIG(CLOCK,NPTS,RATE,luout,YMDHMS,IY,DAY,HRS,MIN,SEC, & AVG,ERROR) C WRITE(luout,107)AVG !Write out digitizing rate write(*,107)AVG C C If there is no error, update the header. Time, date, and C digitizing rate and total seconds will be changed. C IF(.NOT.ERROR)THEN C C Convert from Julian to DD-MMM-YY for date. C CALL JULIAN(DAY,IY,IM,IDY,YMDHMS,.TRUE.,.FALSE.) C C Insert correct time into time-date array. C write(ymdhms(11:21),102) hrs,min,sec IF(YMDHMS(17:17).EQ.' ')YMDHMS(17:17)='0' WRITE(luout,103)YMDHMS,DAY !Write correct date & time write(*,103)YMDHMS,DAY RATE=AVG TOTTIM=FLOAT(IPTS4-1)/RATE ELSE !ERROR EXISTS ymdetc = ymdhms com = 'TIME BAD' IF(MIN.NE.99)THEN !We weren't totally lost, guess time. IF(DAY.NE.999) !Decode Julian day if sensible. & CALL JULIAN(DAY,IY,IM,IDY,YMDETC,.TRUE.,.FALSE.) C C Write estimate for time and date. C write(ymdetc(11:21),102) hrs,min,sec IF(YMDETC(17:17).EQ.' ')YMDETC(17:17)='0' WRITE(luout,105)YMDETC,DAY write(*,105)YMDETC,DAY ELSE !BLEW IT, NO USE GUESSING !! WRITE(luout,106) write(*,106) ENDIF ENDIF 999 return 102 FORMAT(I2.2,':',I2.2,':',F5.2) 103 FORMAT(1X,'The corrected time is: ',A21,5X,'Julian Day:',I4) 104 FORMAT(A8) 105 FORMAT(1X,'The calculated time is: ',A21,5X,'Julian Day:',I4, & ' -- HEADER NOT UPDATED') 106 FORMAT(1X,'UNABLE TO PICK THE TIME ', & ' -- HEADER NOT UPDATED') 107 FORMAT(1X,'The average digitizing rate over the entire series', & ' is: ',F10.3) END C+ SUBROUTINE TIMSIG(CLOCK,NPTS,RATE,luout,YMDHMS,IY,DAY,HRS,MIN, & SEC,AVG,ERROR) C C James A. Bollinger, May 15,1985 C C This subroutine attempts to figure out the absolute UCT time C for the first point in a VTF file based on the digitized signal C from the satellite clock. C input arguments are CLOCK, the time signal C RATE - the target digitizing rate C luout - the unit for writing report C YMDHMS - the trigger time, according to C system clock C Outputs are: TIME - DAY,HRS,MIN,SEC (IY=YEAR) C AVG - the true average digitizing rate. C ERROR - logical error flag C C If ERROR is true, the part of the time that could not be C accurately determined is set = nines in its significant C digits, (i.e. 99 for minutes or 99.99 for seconds). C C- INTEGER*2 CLOCK(*) integer PEAK(130,2),PEKPTR,HINDX,MINPTR,FORTIM INTEGER MIN,HRS,DAY,MINNDX,LOTRIG,HITRIG,I,NPTS,luout,IY LOGICAL*1 DAYS,HOURS,MINUTS,ERROR,SWITCH,INVERT,NEWDAY,LEAPYR character*21 YMDHMS REAL*4 AVG,SEC,RATE C C Initiallize the peak trigger (1600) and the trough trigger (400) C Initiallize the pointers. C SWITCH=.FALSE. INVERT=.FALSE. ERROR=.FALSE. HINDX = NPTS HITRIG=1600 LOTRIG=400 5 PEKPTR=0 C Find all the peaks in the time signal. Make an error check. If C no peaks are found or too many are found -- EXIT! Something is C wrong, and the rest of the routine can't deal with it! C CALL FINDPK(CLOCK,HINDX,RATE,PEAK,PEKPTR,LOTRIG,HITRIG,ERROR) 6 IF(ERROR)THEN CALL MISSED(DAY,HRS,MIN,SEC,ERROR) GO TO 999 ENDIF C C Break the peaks into three groups based on their width in time. C Code: SHORT=20 samples=0; MEDIUM=50 samples=1; LONG=80 =2. C CALL BINS(PEAK(1,2),PEKPTR) C C Check to see whether or not the time signal is inverted. If it is, C reverse it and begin again. SWITCH keeps this process (reversing C and retrying) from going on forever. C CALL CKINVT(PEAK(1,2),PEKPTR,CLOCK,HINDX,INVERT) IF((.NOT.SWITCH).AND.INVERT)THEN SWITCH=.TRUE. WRITE(luout,101) write(*,101) GO TO 5 ELSEIF(SWITCH.AND.INVERT)THEN ERROR=.TRUE. WRITE(luout,102) write(*,102) GO TO 6 ENDIF C C Figure out the true digitizing rate. C CALL GETRAT(PEAK,PEKPTR,AVG) C C Try to find the "double long" pulse which marks the beginning of C the minute. If I can't find the minute mark, exit. C MINPTR=0 DO 10 I=2,PEKPTR IF((PEAK(I,2).EQ.2).AND.(PEAK(I-1,2).EQ.2))THEN MINPTR=I GO TO 20 ENDIF 10 CONTINUE CALL MISSED(DAYS,HRS,MIN,SEC,ERROR) GO TO 999 C C Find out how many seconds are left in the CLOCK signal to the C right of the minute mark. 20 CALL LKAHED(PEKPTR,MINPTR,MINUTS,HOURS,DAYS) C Find out the minutes if possible, otherwise set error flag. C CALL GETMIN(PEAK,MINPTR,MINUTS,MIN,ERROR) C C Find the hours, if possible. C CALL GETHRS(PEAK,MINPTR,HOURS,HRS,MIN,ERROR) C C Figure out whether this minute mark is the beginning of a new day C NEWDAY=(HRS.EQ.0).AND.(MIN.EQ.0) C C Now find the days. C CALL GETDAY(PEAK,MINPTR,DAYS,DAY,NEWDAY,ERROR) C C Try to figure out from the system trigger time what the year is. C CALL GETYRS(YMDHMS,NEWDAY,IY) C C Extrapolate back in time to the first point in the file. C MINNDX=PEAK(MINPTR,1) SECDIF=(FLOAT(MINNDX-1))/AVG CAUG88 IF (SECDIF .GT. 60.0) THEN SECDIF = SECDIF - 60.0 MIN= MIN - 1 END IF CAUG88 C C Try to handle all the marginal time situations. This would C include the first minute in the hour, first hour in the day, C first day in the year, if its a leap year, etc. C CALL MARGIN(IY,DAY,HRS,MIN,SEC,SECDIF,DAYS,NEWDAY,ERROR) 999 RETURN 101 FORMAT(1X,'Time signal is inverted.') 102 FORMAT(/1X,'Inversion yields inconsistent signal - EXITING.') END C+ SUBROUTINE FINDPK(CLOCK,HINDX,RATE,PEAK,PEKPTR,LOTRIG,HITRIG, & ERROR) C C This subroutine finds the number of "peaks" in a time series C and their width. This information is recorded in a two- C dimensional array where the first value in each pair is the C time-series index of the beginning of the peak, and the second C value in the pair is the width of the peak in time counts. C C CLOCK is the time series array -- input parameter C HINDX is the length of the time series in time counts -- input C PEAK is the array containing the output info outlined above C PEKPTR is the number of peaks found in the time series C HITRIG is the trigger level for the transition from low to high C LOTRIG is the trigger level for the transition from high to low C ERROR is the flag set if we have problems. C- INTEGER*2 CLOCK(*) ! Input time series integer HITRIG,LOTRIG INTEGER PEAK(130,2),PEKPTR,TIMPTR,PEKBEG,WIDTH,HINDX LOGICAL*1 ERROR REAL*4 RATE C C Initiallize local pointer. C TIMPTR=1 C C Since getting the tail end of a peak doesn't help, throw away C everything up to the first low signal, then begin processing. 5 IF(CLOCK(TIMPTR).LT.LOTRIG)THEN !SIGNAL BEGINS LOW GO TO 10 !BEGIN PROCESSING NORMALLY ELSE !THROW AWAY UNTIL LOW TIMPTR=TIMPTR+1 IF(TIMPTR.GT.HINDX)GO TO 99 !QUIT IF RUN OUT OF DATA GO TO 5 ENDIF C C Trigger on a peak (Now locally set 1600) 10 TIMPTR=TIMPTR+1 !AT THIS POINT - LOW! IF(TIMPTR.GT.HINDX)GO TO 99 !RAN OUT OF DATA WHILE LOW IF(CLOCK(TIMPTR).GT.HITRIG)THEN !Found a peak PEKPTR=PEKPTR+1 !Increment peak counter PEKBEG=TIMPTR !Set ptr to begin. of peak C C Count the length of the peak. 20 TIMPTR=TIMPTR+1 C Note that partial peaks left over at the end of data are C discarded. C IF(TIMPTR.GT.HINDX)THEN !RAN OUT OF DATA WHILE HI PEKPTR=PEKPTR-1 !DISCARD ERRONEOUS PEAK GO TO 99 !EXIT ENDIF IF(CLOCK(TIMPTR).LT.LOTRIG)THEN !SIGNAL TRIGGERS LOW C C Calculate width of peak, record the info about this peak. C WIDTH=TIMPTR-PEKBEG PEAK(PEKPTR,1)=PEKBEG PEAK(PEKPTR,2)=WIDTH GO TO 10 ELSE !SIGNAL IS STILL HIGH GO TO 20 ENDIF ENDIF GO TO 10 C C Before exit, make sure there aren't too many peaks or too few. C Note error if present. 99 IF((PEKPTR.LT.HINDX/NINT(1.1*RATE)).OR. & (PEKPTR.GT.HINDX/NINT(0.9*RATE))) ERROR=.TRUE. RETURN END C+ SUBROUTINE BINS(WIDTH,PEKPTR) C C BINS takes information from the array PEAKS and breaks C the peaks into three groups by "width" (length in time) - C LONG, MEDIUM, SHORT C C Choosing the 0 code value for short and 1 for medium allows C one to calculate the time sequence more easily since short C pulses are off(=0) and medium indicates on(=1). C C This information is copied back on top of the width to C save space. C- integer WIDTH(*),PEKPTR,I C DO 10 I=1,PEKPTR IF((WIDTH(I).GE.10).AND.(WIDTH(I).LE.35))THEN WIDTH(I)=0 ELSEIF((WIDTH(I).GT.35).AND.(WIDTH(I).LE.65))THEN WIDTH(I)=1 ELSEIF(WIDTH(I).GT.65)THEN WIDTH(I)=2 ELSE WIDTH(I)=0 ENDIF 10 CONTINUE RETURN END C+ SUBROUTINE GETRAT(PEAK,PEKPTR,AVGRAT) C C This subroutine uses information from the array PEAKS, which C is created in the subroutine FINDPK to calculate the digitizing C rate by ONE methods. Method 1 is the average rate, which is the C total samples divided by the total time (AVGRAT). C- integer PEAK(130,2),PEKPTR REAL*4 AVGRAT,DIGINC,DELTAT C C Compute the average digitizing rate. C DELTAT=FLOAT(PEKPTR-1) DIGINC=FLOAT(PEAK(PEKPTR,1)-PEAK(1,1)) AVGRAT=DIGINC/DELTAT RETURN END C+ SUBROUTINE CKINVT(PEAK,PEKPTR,CLOCK,HINDX,INVERT) C C This subroutine checks to see if the signal is inverted C by looking for three wide "peaks" in row. C- LOGICAL*1 INVERT integer PEAK(*),PEKPTR,I,HINDX integer*2 CLOCK(*) C DO 10 I=3,PEKPTR IF((PEAK(I).EQ.2).AND.(PEAK(I-1).EQ.2).AND. & (PEAK(I-2).EQ.2))THEN INVERT=.TRUE. GO TO 20 ENDIF 10 CONTINUE INVERT=.FALSE. 20 IF(INVERT)THEN DO 30 I=1,HINDX CLOCK(I)=2048-CLOCK(I) 30 CONTINUE ENDIF RETURN END C+ SUBROUTINE MISSED(DAYS,HRS,MIN,SEC,ERROR) C C MISSED sets the date and time to values that are out-of-range C to flag error conditions. C- integer DAYS,HRS,MIN REAL*4 SEC LOGICAL*1 ERROR C ERROR=.TRUE. DAYS=999 HRS=99 MIN=99 SEC=99.99 RETURN END C+ SUBROUTINE LKAHED(PEKPTR,MINPTR,MINUTS,HOURS,DAYS) C C The next few subroutines calculate the Julian day and time from C the time code. To perform this task, these subroutines need C information about how many seconds there are between the minute C mark and the end of the signal in order to decide whether to try C to calculate the time from code ahead or behind the minute mark C in time. LKAHED looks ahead of the minute mark and decides which C information can be calculated after the minute mark. C- integer PEKPTR,MINPTR LOGICAL*1 MINUTS,HOURS,DAYS C FORTIM=PEKPTR-MINPTR MINUTS=.FALSE. HOURS=.FALSE. DAYS=.FALSE. IF((FORTIM.GE.17).AND.(FORTIM.LT.26))THEN MINUTS=.TRUE. ELSEIF((FORTIM.GE.26).AND.(FORTIM.LT.41))THEN MINUTS=.TRUE. HOURS=.TRUE. ELSEIF(FORTIM.GE.41)THEN MINUTS=.TRUE. HOURS=.TRUE. DAYS=.TRUE. ENDIF RETURN END C+ SUBROUTINE GETMIN(PEAK,MINPTR,MINUTS,MIN,ERROR) C C GETMIN tries to calculate the minute from the time code. C- integer MIN,PEAK(130,2),MINPTR LOGICAL*1 MINUTS,ERROR C IF(MINUTS)THEN MIN=PEAK(MINPTR+10,2)+PEAK(MINPTR+11,2)*2+PEAK(MINPTR+12,2)*4 & +PEAK(MINPTR+13,2)*8+PEAK(MINPTR+15,2)*10+PEAK(MINPTR+16,2) 2 *20+PEAK(MINPTR+17,2)*40 ELSEIF(MINPTR-50.GE.1)THEN MIN=PEAK(MINPTR-50,2)+PEAK(MINPTR-49,2)*2+PEAK(MINPTR-48,2)*4 & +PEAK(MINPTR-47,2)*8+PEAK(MINPTR-45,2)*10+PEAK(MINPTR-44,2) 2 *20+PEAK(MINPTR-43,2)*40+1 ELSE ERROR=.TRUE. MIN=99 GO TO 99 ENDIF MIN=MOD(MIN,60) 99 RETURN END C+ SUBROUTINE GETHRS(PEAK,MINPTR,HOURS,HRS,MIN,ERROR) C C GETHRS computes the hour from the time code if possible. C- integer HRS,PEAK(130,2),MINPTR,MIN LOGICAL*1 ERROR,HOURS C IF(HOURS)THEN HRS=PEAK(MINPTR+20,2)+PEAK(MINPTR+21,2)*2+PEAK(MINPTR+22,2)*4 & +PEAK(MINPTR+23,2)*8+PEAK(MINPTR+25,2)*10+PEAK(MINPTR+26,2)*20 ELSEIF((MINPTR-40.GE.1).AND.(MIN.EQ.0))THEN HRS=PEAK(MINPTR-40,2)+PEAK(MINPTR-39,2)*2+PEAK(MINPTR-38,2)*4 & +PEAK(MINPTR-37,2)*8+PEAK(MINPTR-35,2)*10+PEAK(MINPTR-34,2)*20 2 +1 ELSEIF(MINPTR-40.GE.1)THEN HRS=PEAK(MINPTR-40,2)+PEAK(MINPTR-39,2)*2+PEAK(MINPTR-38,2)*4 & +PEAK(MINPTR-37,2)*8+PEAK(MINPTR-35,2)*10+PEAK(MINPTR-34,2)*20 ELSE ERROR=.TRUE. HRS=99 GO TO 99 ENDIF HRS=MOD(HRS,24) 99 RETURN END C+ SUBROUTINE GETDAY(PEAK,MINPTR,DAYS,DAY,NEWDAY,ERROR) C C GETDAY calculates the day from the time code if possible. C- integer PEAK(130,2),MINPTR,DAY LOGICAL*1 DAYS,ERROR,NEWDAY C IF(DAYS)THEN DAY=PEAK(MINPTR+30,2)+PEAK(MINPTR+31,2)*2+PEAK(MINPTR+32,2)*4 & +PEAK(MINPTR+33,2)*8+PEAK(MINPTR+35,2)*10+PEAK(MINPTR+36,2) 2 *20+PEAK(MINPTR+37,2)*40+PEAK(MINPTR+38,2)*80+PEAK(MINPTR+40 3 ,2)*100+PEAK(MINPTR+41,2)*200 ELSEIF((MINPTR-30.GE.1).AND.NEWDAY)THEN DAY=PEAK(MINPTR-30,2)+PEAK(MINPTR-29,2)*2+PEAK(MINPTR-28,2)*4 & +PEAK(MINPTR-27,2)*8+PEAK(MINPTR-25,2)*10+PEAK(MINPTR-24,2) 2 *20+PEAK(MINPTR-23,2)*40+PEAK(MINPTR-22,2)*80+PEAK(MINPTR-20 3 ,2)*100+PEAK(MINPTR-19,2)*200+1 ELSEIF(MINPTR-30.GE.1)THEN DAY=PEAK(MINPTR-30,2)+PEAK(MINPTR-29,2)*2+PEAK(MINPTR-28,2)*4 & +PEAK(MINPTR-27,2)*8+PEAK(MINPTR-25,2)*10+PEAK(MINPTR-24,2) 2 *20+PEAK(MINPTR-23,2)*40+PEAK(MINPTR-22,2)*80+PEAK(MINPTR-20 3 ,2)*100+PEAK(MINPTR-19,2)*200 ELSE ERROR=.TRUE. DAY=999 GO TO 99 ENDIF 99 RETURN END C+ SUBROUTINE GETYRS(YMDHMS,NEWDAY,YEAR) C C This subroutine looks at the year returned in the header and C only messes around with it if it's New Years Day. C- character*21 YMDHMS LOGICAL*1 NEWDAY integer YEAR CHARACTER*6 SYSDAY C read(ymdhms(1:10),'(a6,1x,i2)') sysday,year IF((SYSDAY.EQ.'01-JAN').AND.NEWDAY)YEAR=YEAR-1 RETURN END C+ SUBROUTINE MARGIN(YEAR,DAY,HRS,MIN,SEC,SECDIF,DAYS,NEWDAY,ERROR) C C This subroutine is supposed to handle marginal cases in time, C like the first minute in an hour, first minute in day, first C minute in year, etc. I have such great confidence in this C algorithm that I put it in a subroutine so that it would be C easier to change. C- REAL*4 SEC,SECDIF integer YEAR,DAY,HRS,MIN LOGICAL*1 NEWDAY,DAYS,ERROR,LEAPYR C LEAPYR=MOD(IY+1900,4).EQ.0 IF(.NOT.ERROR)THEN IF(NEWDAY)THEN IF(DAYS.AND.(DAY.EQ.1))THEN IF(LEAPYR)THEN DAY=366 ELSE DAY=365 ENDIF ELSE DAY=DAY-1 ENDIF HRS=23 MIN=59 ELSEIF(MIN.EQ.0)THEN HRS=HRS-1 MIN=59 ELSE MIN=MIN-1 ENDIF SEC=60.0-SECDIF ELSEIF((MIN.NE.99).AND.(HRS.NE.99))THEN IF(MIN.EQ.0)THEN HRS=HRS-1 MIN=59 ELSE MIN=MIN-1 ENDIF SEC=60.0-SECDIF ELSE SEC=99.99 ENDIF RETURN END C+ subroutine get_suds_structure(nunit,contents,nhead,ndata, 1 buffer,nstruct,nleft,kb) c c Reads in next suds structure header from a 512-byte/block unix C file. nhead is the length of the structure header, ndata is C the number of byte is the data part (the total actually C takes 12 more bytes than the sum in original storage), nstruct C is the SUDS structure number, contents is a byte array of C dimension nhead. nleft is the number of bytes in the current C block left over from the current structure. If nleft is zero, C starts by reading a new record=block. C buffer is the contents of the current block. It is used by the C calling program if ndata is greater than zero. C nunit is the logical unit number of the input file. C- character contents*512, buffer*512,residual*12 integer*2 nstruct integer*4 nhead,ndata C if (nleft .eq. 0) then kb = kb + 1 read(nunit'kb) buffer read(buffer(3:12),'(a2,2a4)') nstruct, nhead, ndata ct write(2,*) nstruct,nhead,ndata,nleft,' first one' if (nstruct .le. 0) then nleft = -1 return end if nleft = 500 else if (nleft .ge. 12) then ct write(2,*) 'About to read from',515-nleft,524-nleft read(buffer(515-nleft:524-nleft),'(a2,2a4)') 1 nstruct, nhead, ndata ct write(2,*) nstruct,nhead,ndata,nleft,' second one' if (nstruct .le. 0) then nleft = -1 return end if nleft = nleft - 12 else residual(1:nleft) = buffer(513-nleft:512) kb = kb + 1 read(nunit'kb,iostat=ierr) buffer if (ierr .ne. 0) then nleft = -1 return end if residual(nleft+1:12) = buffer(1:12-nleft) read(residual(3:12),'(a2,2a4)') nstruct, nhead, ndata ct write(2,*) nstruct,nhead,ndata,nleft,' third one' if (nstruct .le. 0) then nleft = -1 return end if nleft = 500 + nleft end if c now = 0 do while (now .lt. nhead) if (nhead-now .le. nleft) then na = 512 - nleft do j=now+1,nhead na = na + 1 contents(j:j) = buffer(na:na) end do nleft = 512 - na now = nhead else ct write(2,*) nhead,nleft,now,' Just before contents in get_suds' contents(now+1:now+nleft) = buffer(513-nleft:512) kb = kb + 1 read(nunit'kb,iostat=ierr) buffer if (ierr .ne. 0) then nleft = -1 return end if now = now + nleft nleft = 512 end if end do return end C+ subroutine move_suds_data(LUIN,LUOUT,nleft_in,nleft_out,buffer, 1 nblocks,now,nch,kb) C C Input from a st.*** file on unit LUIN, data section starting at C 256-nleft_in (I*2) for the next 256*nblocks I*2 words ending at C position nleft_out. The initial C bunch may be already read, hence passing buffer. Writes to C LUOUT. now is the current block number in LUOUT. C All channels but #1 (time) and #5 (blab) are flipped in polarity C- integer*2 i2data(256),i2d(256) character*512 buffer ct if (nch .lt. 3) ct 1 write(2,*) 'nleft in, out nch,kb,now',nleft_in,nleft_out,nch,kb,now if (nleft_in .eq. 0 .and. nleft_out .eq. 0) then do j=1,nblocks kb = kb + 1 read(luin'kb) i2data do k=1,256 if (.not.(nch.eq.5 .or. nch.eq.1)) then if (i2data(k) .eq. -32768) i2data(k) = -32767 i2data(k) = -i2data(k) end if end do now = now + 1 write(luout'now) i2data end do else nst = 257 - nleft_in nend = nst - 1 ntop = 256 ct if (nch .lt. 3) write(2,*) 'nst,nend',nst,nend do j=1,nblocks if (j .eq. nblocks) then if (nleft_out .lt. nleft_in) then ntop = 256 - nleft_out else nend = 256 - nleft_out end if end if nn = ntop - nst + 1 ct if (nch .lt. 3) write(2,*) 'j,nst,nn,ntop,',j,nst,nn,ntop read(luin'kb) i2d do k=1,nn i2data(k) = i2d(nst-1+k) end do ct read(buffer(2*nst-1:2*ntop),'(256a2)') (i2data(k),k=1,nn) if ((j .lt. nblocks) .or. (nleft_out .ge. nleft_in)) then kb = kb + 1 read(luin'kb) buffer read(luin'kb) i2d do k=1,nend i2data(nleft_in+k) = i2d(k) end do ct read(buffer(1:2*nend),'(256a2)') (i2data(k),k=nst,ntop) end if do k=1,256 if (.not.(nch.eq.5 .or. nch.eq.1)) then if (i2data(k) .eq. -32768) i2data(k) = -32767 i2data(k) = -i2data(k) end if end do now = now + 1 write(luout'now) i2data end do end if return end