FreeUSP: Algorithm Implementation and Testing
Trace Template
c ----------------- Main Routine -----------------------
c
c copyright 2001-6, Amoco Production Company
c All rights reserved
c an affiliate of BP America Inc.
c
c ----------------- ------------ -----------------------
c Program Description:
c FORTRAN 77 template for trace processing in USP.
c Author[s]: P.Gutowski, J.Wade, P.Garossino, J.Dellinger
c Program Changes:
c - original written: June 26, 2001
implicit none
c Define various machine-dependent parameters.
c These include files are inserted into the code by the compiler;
c they can be found under "~usp/include".
#include <f77/iounit.h>
#include <f77/lhdrsz.h>
#include <f77/sisdef.h>
#include <save_defs.h>
c dimension standard USP variables
integer nsamp, nsampo, nsi, ntrc, ntrco, nrec, nreco, iform
integer luin, luout, lbytes, nbytes, lbyout, obytes
integer ist, iend, irs, ire, ns, ne, argis, jerr
integer JJ, KK
real UnitSc, dt
character ntap*512 , otap*512, name*4
logical verbos
c Program Specific _ dynamic memory variables
integer errcd, errcd1, alloc_size, abort
integer itr
real Trace_WorkSpace, tri
c The two in parenthesis in the pointer statements
c indicates the number of elements that the debugger will be able to
c see in the main routine while debugging the code. It has nothing
c to do with the actual amount of memory allocated or requested.
c If you want to look at lots of the array then change the two
c to something larger, like 100000.
pointer ( mem_Trace_WorkSpace, Trace_WorkSpace(2))
pointer ( mem_itr, itr(2))
pointer ( mem_tri, tri(2))
c Program Specific _ static memory variables
integer ifmt_RecNum,l_RecNum,ln_RecNum, RecNum
integer ifmt_TrcNum,l_TrcNum,ln_TrcNum, TrcNum
integer ifmt_DstSgn,l_DstSgn,ln_DstSgn, DstSgn
integer ifmt_StaCor,l_StaCor,ln_StaCor, StaCor
c Initialize variables
data abort/0/
data name/"PRGM"/
c give command line help if requested
if ( argis ( '-?' ) .gt. 0 .or.
: argis ( '-h' ) .gt. 0 .or.
: argis ( '-help' ) .gt. 0 ) then
call help()
stop
endif
c open printout file in preparation
c for processing. The output file consists of a 2-character
c identifier, the program name followed by the process id
c number: PRGM.nnnnn. The file will contain information
c important to debugging any problems that may crop up
c during execution.
c
#include <f77/open.h>
c Scan the command line for arguments.
call cmdln ( ntap, otap, ns, ne, irs, ire, ist, iend, name,
: verbos )
c-----
c Get logical unit numbers for input and output of seismic data
c using the USP library routine "getln".
c
c Input values are strings ntap and otap (names for the input and
c output files, respectively), the read 'r' or write 'w'
c designations of these files, and the default logical unit numbers
c to use if ntap and/or otap are blank strings. Usually you will
c want the default unit numbers to be:
c
c 0 = stdin for reading
c 1 = stdout for writing
c
c The output value (getln's first argument) is the SIS I/O
c "logical unit number" to use for accessing the file. If the
c returned number is less than 0 it means there was a fatal error
c opening the corresponding file. (Usually getln will not get the
c chance to return in case of error, however, because the lower-level
c SIS I/O routines will have aborted the entire program the instant
c the error occurred.)
c
c Note that SIS I/O "logical unit numbers" are NOT the same things
c as Fortran I/O units, nor are they the same as C file descriptors.
c A logical unit number returned by getln can ONLY be used in
c conjunction with SIS I/O calls. If you need to do raw Fortran I/O,
c use the standard Fortran library routines to open, read, write, and
c close the associated file. In particular, use the FORTRAN routine
c "alloclun" to come up with an unused unit number (since you can't
c easily find out what unit numbers SIS I/O might be using). (The USP
c library routine "lucfd" can be used to find the C file descriptor
c associated with an SIS I/O logical unit number, but this is
c generally not useful for Fortran programmers.) If for your
c application you'll need to do a lot of fancy text processing or
c low-level (non-USP) I/O, consider writing your program in C
c instead of FORTRAN. (SIS I/O itself, for example, is primarily
c written in C.)
c
c In the getln call below, note that if the user did not specify
c -N or -O on the command line, ntap and otap will be blanks, and
c the call will set up stdin and stdout for input and output
c (respectively).
c
c getln can actually use any of several read/write modes:
c 'r' opens an already existing file for reading only;
c 'r+' can be used to open an already existing file for both
c reading and writing;
c 'w' opens a file for writing only (destroying any file of
c that name that may already have been there);
c 'w+' opens a file for reading and writing (destroying any
c file of that name that may already have been there);
c 'a' opens an already existing file for both reading and writing,
c with reading starting at the beginning and writing starting
c at the end (i.e., writing can only append);
c 'a+' is the same as 'a' but with both reading and writing
c starting at the end.
call getln(luin , ntap,'r', 0)
call getln(luout, otap,'w', 1)
c allocate memory for line header buffer itr[]
alloc_size = SZLNHD
errcd = 0
call galloc(mem_itr,alloc_size * SZSMPD,errcd,abort)
if (errcd .ne. 0) then
write(LERR,*) 'ERROR: Unable to allocate workspace '
write(LERR,*) ' ',alloc_size * SZSMPD,' bytes requested '
write(LERR,*) ' FATAL'
write(LER,*) 'PRGM: Unable to allocate workspace '
write(LER,*) ' ',alloc_size * SZSMPD,' bytes requested '
write(LER,*) 'FATAL'
goto 999
else
write(LERR,*)'Allocating workspace:'
write(LERR,*) alloc_size * SZSMPD,' bytes requested '
endif
call vclr ( itr, 1, alloc_size )
c Read the input dataset's line header, which comprises the
c very first record in a USP dataset.
c Among other useful things, the line header tells us the
c dimensions of the dataset we're about to process.
c
c One call to the standard input subroutine rtape reads one
c record of data. For USP datasets, this is either a line header
c or a trace header and trace. Here rtape reads data into the
c array itr; lbytes is the number of bytes actually read.
lbytes = 0
call rtape(luin,itr,lbytes)
c----
c
c About output streams
c----------------------------------------------------------------
c
c There are three output streams, LOT, LER, and LERR.
c
c LOT is standard out, probably where the output usp data file
c is going. You certainly don't want to send any error messages
c there!
c
c LER is standard error; that's where messages should go if you
c want the user to see them on their terminal. (Standard error is
c NOT buffered: messages appear on the user's screen as soon as
c the call to write occurs.)
c
c LERR is the USP log file. Note that output to the log file
c is buffered, so there's no guarantee that messages written
c there will actually go out at the time of your call to
c write. They may not get written until the program exits,
c and if the program exits abnormally they may never get written
c at all!
c
c Remember to always include the program name in any error
c message, so the user knows which program it came from!
c
c----
c
c Complain if we can't read the line header.
if(lbytes.eq.0)then
write(LERR,*)'PRGM: no line header on input dataset',
: ntap
write(LER,*)' '
write(LER,*)'PRGM: '
write(LER,*)' no line header on input dataset',ntap
write(LER,*)'FATAL'
write(LER,*)' '
goto 999
endif
c Here we use saver to find the basic dimensions of the dataset
c as specified in the line header. (Examples of using other
c header-access library routines can be found further along in the
c code.)
c
c LINHED and TRACEHEADER (and more not used in this program)
c are defined in the include file <save_defs.h>.
c
c Make sure to specify keywords in single quotes. FORTRAN
c converts strings in double quotes to upper case; saver
c and related subroutines won't recognize those!
call saver(itr, 'NumSmp', nsamp, LINHED)
call saver(itr, 'SmpInt', nsi , LINHED)
call saver(itr, 'NumTrc', ntrc , LINHED)
call saver(itr, 'NumRec', nrec , LINHED)
call saver(itr, 'Format', iform, LINHED)
call saver(itr, 'UnitSc', UnitSc, LINHED)
c read historical line header and print to printout file this routine
c also write the program name to the historical line header
call hlhprt (itr, lbytes, name, 4, LERR)
c POLICEMAN: check header for units scaling. Using UnitSc, remember
c that UnitSc default is milliseconds [i.e. 0.001] and UnitSc
c is a floating point variable. A UnitSc entry of 1.0 would
c mean units are in seconds. A UnitSc entry of 0 indicates that
c the unit was not defined. In this case milliseconds are
c assumed and loaded to the header for further processing.
if ( UnitSc .eq. 0.0 ) then
write(LERR,*)'********************************************'
write(LERR,*)'WARNING: sample unit scaler in LH = ',UnitSc
write(LERR,*)' will set to .001 (millisec default)'
write(LERR,*)'********************************************'
UnitSc = 0.001
call savew ( itr, 'UnitSc', UnitSc, LINHED)
endif
c compute delta T in seconds
dt = real (nsi) * UnitSc
c check user supplied boundary conditions and set defaults
if ( irs .eq. 0 ) irs = 1
if ( ire .eq. 0 .or. ire .gt. nrec ) ire = nrec
if ( ns .eq. 0 ) ns = 1
if ( ne .eq. 0 .or. ne .gt. ntrc ) ne = ntrc
c this logic assumes parameterization in units of the dataset. It is
c also assumed here that sample one is time zero. In your code you may
c want to consider the line header entry TmMsFS, if it is important
c that this parameter be used in terms of time based on the original
c field recording. If the data was previously windowed, then TmMsFS
c will be non-zero and may need to be considered here. You will
c also have to update TmMsFS on output of the new line header.
c convert start and end times to nearest sample values
ist = nint ( float(ist) / float(nsi) ) + 1
if ( ist .eq. 0 ) ist = 1
if ( iend .eq. -99999 ) then
iend = nsamp
else
iend = nint ( float(iend) / float(nsi) ) + 1
if ( iend .eq. 0 .or. iend .gt. nsamp ) iend = nsamp
endif
nreco = ire - irs + 1
ntrco = ne - ns + 1
nsampo = iend - ist + 1
c modify line header to reflect actual record configuration output
c NOTE: in this case the sample limits ist and iend are used to
c limit processing only. All data samples are actually passed.
call savew(itr, 'NumRec', nreco, LINHED)
call savew(itr, 'NumTrc', ntrco, LINHED)
call savew(itr, 'NumSmp, nsampo, LINHED)
c number output bytes
obytes = SZTRHD + SZSMPD * nsampo
c Now inject the current command line into the historical line
c header as well. For historical reasons savhlh outputs
c the modified line-header length separately in the 3rd argument.
call savhlh ( itr, lbytes, lbyout )
c We've finished modifying the line header; write it out.
c (Write to unit luout the lbyout bytes in the array itr.)
call wrtape ( luout, itr, lbyout )
c verbose output of all pertinent information before processing begins
call verbal( ntap, otap, nsamp, nsi, ntrc, nrec, iform, ist,
: iend, irs, ire, ns, ne, verbos)
c Here is an example of how to use savelu to find out the
c format, offset and length of various trace header parameters to be
c used later in the program.
call savelu('RecNum',ifmt_RecNum,l_RecNum,ln_RecNum,TRACEHEADER)
call savelu('TrcNum',ifmt_TrcNum,l_TrcNum,ln_TrcNum,TRACEHEADER)
call savelu('DstSgn',ifmt_DstSgn,l_DstSgn,ln_DstSgn,TRACEHEADER)
call savelu('StaCor',ifmt_StaCor,l_StaCor,ln_StaCor,TRACEHEADER)
c reallocate buffer memory to be big enough to include the trace header
c and trace time series
alloc_size = obytes
if (nsamp .gt. nsampo) alloc_size = SZTRHD + SZSMPD * nsamp
errcd = 0
call grealloc(mem_itr,alloc_size,errcd,abort)
if (errcd .ne. 0) then
write(LERR,*) 'ERROR: Unable to allocate trace workspace '
write(LERR,*) ' ',alloc_size,' bytes requested '
write(LERR,*) ' FATAL'
write(LER,*) 'PRGM: Unable to allocate trace workspace '
write(LER,*) ' ',alloc_size,' bytes requested '
write(LER,*) 'FATAL'
goto 999
else
write(LERR,*)'Allocating workspace:'
write(LERR,*) alloc_size,' bytes requested '
endif
call vclr ( itr, 1, nsamp )
c dynamic memory allocation:
call galloc (mem_tri, nsamp * SZSMPD, errcd, abort)
call galloc (mem_Trace_WorkSpace, nsamp * SZSMPD, errcd1, abort)
if ( errcd .ne. 0 .or.
: errcd1 .ne. 0 )then
write(LERR,*)' '
write(LERR,*)'Unable to allocate workspace:'
write(LERR,*) 2 * nsamp * SZSMPD, ' bytes'
write(LERR,*)' '
write(LER,*)' '
write(LER,*)'Unable to allocate workspace:'
write(LER,*) 2 * nsamp * SZSMPD, ' bytes'
write(LER,*)' '
go to 999
else
write(LERR,*)' '
write(LERR,*)'Allocating workspace:'
write(LERR,*) 2 * nsamp * SZSMPD, ' bytes'
write(LERR,*)' '
endif
c initialize memory
call vclr ( tri, 1, nsamp )
call vclr ( Trace_WorkSpace, 1, nsamp )
c BEGIN PROCESSING
c Skip over unwanted initial records.
c recskp reads along a trace at a time for a precalculated
c number of traces (it doesn't actually look at the record
c number to verify that it has stopped on the right one).
c Records are assumed to start with number 1, not 0. The
c first argument specifies the current record, the second
c the record to stop on. Although only the difference between
c these is actually used, if either the starting or ending
c record numbers are 0 recskp silently aborts!
c The third argument is the logical unit to read from;
c the fourth is the number of traces in a record.
c Space sufficient to store a trace in (here itr) must be
c passed to recskp in the final argument for use as
c temporary scratch space. (Upon return it may or may not
c contain a copy of the last skipped-over trace.)
call recskp ( 1, irs-1, luin, ntrc, itr )
DO JJ = irs, ire
c We're at the start of the correct record ... now
c skip to the desired starting trace within that record.
c
c In this call JJ is the current record number. (Depending on
c whether trcskp thinks the input is a pipe or a file it may or may
c not rely on this information, so be careful when debugging!)
c ns-1 is the trace number within that record to stop on,
c luin is the logical unit,
c ntrc is the number of traces in a record, and
c itr is again temporary scratch space sufficient
c for storing a trace.
call trcskp ( JJ, 1, ns-1, luin, ntrc, itr )
DO KK = ns, ne
c read a trace
nbytes = 0
call rtape( luin, itr, nbytes)
c if end of data encountered (nbytes=0) then bail out
if(nbytes .eq. 0) then
write(LERR,*)'Premature EOF on input at:'
write(LERR,*)' rec= ',JJ,' trace= ',KK
go to 999
endif
c Use previously derived pointers to get trace header values
c efficiently.
call saver2( itr, ifmt_RecNum, l_RecNum, ln_RecNum,
: RecNum, TRACEHEADER )
call saver2( itr, ifmt_TrcNum, l_TrcNum, ln_TrcNum,
: TrcNum, TRACEHEADER)
call saver2( itr, ifmt_DstSgn, l_DstSgn, ln_DstSgn,
: DstSgn, TRACEHEADER )
call saver2( itr, ifmt_StaCor, l_StaCor, ln_StaCor,
: StaCor, TRACEHEADER )
c process only live traces. By convention, a "Station Correction"
c of 30000 is used as a dead trace flag!
if ( StaCor .ne. 30000) then
c Copy the data part of the trace across to tri;
c
c CALL VMOV (A,IA,C,IC,N)
c where,
c A Real input vector.
c IA Integer input stride for vector A.
c C Real output vector.
c IC Integer input stride for vector C.
c N Integer input element count.
call vmov ( itr(ITHWP1), 1, tri, 1, nsamp )
c Put your subroutine here [remember to declare any arguments you need
c over and above those already declared]
c load trace portion of itr[] to real array tri[]
call subs( tri, Trace_WorkSpace, nsamp, ist, iend,
: nsampo )
c Copy the data part of the trace back to the output buffer itr
call vmov ( tri, 1, itr(ITHWP1), 1, nsampo )
endif
c write output data
call wrtape (luout, itr, obytes)
ENDDO
c skip to end of record
call trcskp ( JJ, ne+1, ntrc, luin, ntrc, itr )
ENDDO
c close data files and notify user of normal termination
c This flushes any data left in the output
c buffer. If the output buffer is NOT closed you can sometimes
c end up with missing data, because data left in the buffer
c at program exit may fail to be written.
call lbclos ( luin )
call lbclos ( luout )
write(LERR,*)'prgm: Normal Termination'
write(LER,*)'prgm: Normal Termination'
stop
999 continue
c close data files and warn user of abnormal termination.
call lbclos ( luin )
call lbclos ( luout )
write(LERR,*)'prgm: ABNORMAL Termination'
write(LER,*)'prgm: ABNORMAL Termination'
stop
end
c ----------------- Subroutine -----------------------
subroutine help()
c provide terse online help [detailed help goes in man page]
#include <f77/iounit.h>
write(LER,*)' '
write(LER,*)'===================================================='
write(LER,*)' '
write(LER,*)' Command Line Arguments for prgm: USP template'
write(LER,*)' '
write(LER,*)' For a more detailed description of these parameters'
write(LER,*)' see the online man page using uman or see the USP'
write(LER,*)' intranet site '
write(LER,*)' '
write(LER,*)'Input......................................... (def)'
write(LER,*)' '
write(LER,*)'-N[] -- input data set (stdin)'
write(LER,*)'-O[] -- output data set (stdout)'
write(LER,*)'-s[] -- start time (ms) (0)'
write(LER,*)'-e[] -- end time (ms) (last sample)'
write(LER,*)'-ns[] -- start trace number (1)'
write(LER,*)'-ne[] -- end trace number (last trace)'
write(LER,*)'-rs[] -- start record (1)'
write(LER,*)'-re[] -- end record (last record)'
write(LER,*)'-V -- verbos printout'
write(LER,*)' '
write(LER,*)'Usage:'
write(LER,*)' prgm -N[] -O[] -s[] -e[] -ns[] -ne[] -rs[]'
write(LER,*)' -re[] -V'
write(LER,*)' '
write(LER,*)'===================================================='
return
end
c ----------------- Subroutine -----------------------
c pick up command line arguments
subroutine cmdln ( ntap, otap, ns, ne, irs, ire, ist, iend,
: name, verbos )
implicit none
#include <f77/iounit.h>
integer ist, iend, ns, ne, irs, ire, argis
character ntap*(*), otap*(*), name*(*)
logical verbos
call argi4 ( '-e', iend, -99999, -99999 )
call argi4 ( '-ne', ne, 0, 0 )
call argi4 ( '-ns', ns, 0, 0 )
call argstr ( '-N', ntap, ' ', ' ' )
call argstr ( '-O', otap, ' ', ' ' )
call argi4 ( '-re', ire, 0, 0 )
call argi4 ( '-rs', irs, 0, 0 )
call argi4 ( '-s', ist, 0, 0 )
verbos = (argis('-V') .gt. 0)
c check for extraneous arguments and abort if found to
c catch all manner of user typo's (eat,die)
call xtrarg ( name, ler, .FALSE., .FALSE. )
call xtrarg ( name, lerr, .FALSE., .TRUE. )
return
end
c ----------------- Subroutine -----------------------
c verbal printout of pertinent program particulars
subroutine verbal ( ntap, otap, nsamp, nsi, ntrc, nrec, iform,
: ist, iend, irs, ire, ns, ne, verbos)
#include <f77/iounit.h>
integer nsamp, ntrc, iform, ist, iend, irs, ire, ns, ne, nsi
character ntap*(*), otap*(*)
logical verbos
write(LERR,*)' '
write(LERR,*)' Input Line Header Parameters'
write(LERR,*)' '
write(LERR,*) ' input data set name = ', ntap
write(LERR,*) ' samples per trace = ', nsamp
write(LERR,*) ' traces per record = ', ntrc
write(LERR,*) ' number of records = ', nrec
write(LERR,*) ' data format = ', iform
write(LERR,*) ' sample interval = ', nsi
write(LERR,*)' '
write(LERR,*)' '
write(LERR,*)' Command Line Parameters '
write(LERR,*)' '
write(LERR,*) ' output data set name = ', otap
write(LERR,*) ' start record = ', irs
write(LERR,*) ' end record = ', irs
write(LERR,*) ' start trace = ', ns
write(LERR,*) ' end trace = ', ne
write(LERR,*) ' sample start = ', ist
write(LERR,*) ' sample end = ', iend
if( verbos ) write(LERR,*) ' verbose printout requested'
write(LERR,*)' '
write(LERR,*)'========================================== '
write(LERR,*)' '
return
end
