FreeUSP: Algorithm Implementation and Testing
Record 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 record processing in USP.
c Author[s]: P.Gutowski, J.Wade, P.Garossino, J.Dellinger
c Program Changes:
c - original written: June 26, 2001
c Sept 19, 2006: Added dynamic memory allocation for initial itr[]
c buffer. Added nsampo usage, increased size of
c ntap[] and otap[]
c Garossino
c
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, nrec, nreco, iform, obytes
integer luin, luout, lbytes, nbytes, lbyout
integer ist, iend, irs, ire, argis, jerr
real UnitSc, dt
character ntap*512, otap*512, name*4
logical verbos
c Program Specific dynamic memory variables
integer RecordSize, HeaderSize, errcd, errcd1, abort, alloc_size
integer Headers
integer itr
real Record, Record_WorkSpace
pointer (mem_Record, Record(2))
pointer (mem_Space, Record_WorkSpace(2))
pointer (mem_Headers, Headers(2))
pointer (mem_itr, itr(2))
c Program Specific static memory variables
integer ifmt_StaCor,l_StaCor,ln_StaCor, StaCor
integer hdr_index, tr_index, JJ, KK
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 process id number: PRGM.nnnnn
c this should be unique even for multiple occurences of the same
c process in a pipeline.
#include <f77/open.h>
c Scan the command line for arguments (the routine gcmdln is not
c a library routine; it is defined near the end of this source file).
call cmdln ( ntap, otap, 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
c or low-level (non-USP) I/O, consider writing your program in
c C instead of FORTRAN. (SIS I/O itself, for example, is
c primarily 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 read of line header
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 print HLH to printout file and add prgm line to HLH
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
c this parameterization assumes input 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 also
c have to update TmMsFS on output of the new line header.
ist = nint ( float(ist) / float(nsi) ) + 1
if ( ist .le. 0 ) ist = 1
if ( iend .eq. -99999 ) then
iend = nsamp
else
iend = nint ( float(iend) / float(nsi) ) + 1
if ( iend .le. 0 .or. iend .gt. nsamp ) iend = nsamp
endif
nreco = ire - irs + 1
nsampo = iend - ist + 1
c modify line header to reflect actual record configuration output
c NOTE: in this case the trace and sample limits are used to
c limit processing only. All data within the selected record
c range are actually passed.
call savew(itr, 'NumRec', nreco, 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, 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 ( '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'
stop
else
write(LERR,*)'Allocating workspace:'
write(LERR,*) alloc_size,' bytes requested '
endif
call vclr ( itr, 1, nsamp )
c dynamic memory allocation:
RecordSize = ntrc * nsamp
HeaderSize = ntrc * ITRWRD
call galloc (mem_Record, RecordSize * SZSMPD, errcd1, abort)
call galloc (mem_Space, RecordSize * SZSMPD, errcd2, abort)
call galloc (mem_Headers, HeaderSize * SZSMPD, errcd3, abort)
if ( errcd1 .ne. 0 .or.
: errcd2 .ne. 0 .or.
: errcd3 .ne. 0 )then
write(LERR,*)' '
write(LERR,*)'Unable to allocate workspace:'
write(LERR,*) 2*RecordSize+HeaderSize* SZSMPD, ' bytes'
write(LERR,*)' '
write(LER,*)' '
write(LER,*)'Unable to allocate workspace:'
write(LER,*) 2*RecordSize+HeaderSize* SZSMPD, ' bytes'
write(LER,*)' '
go to 999
else
write(LERR,*)' '
write(LERR,*)'Allocating workspace:'
write(LERR,*) 2*RecordSize+HeaderSize* SZSMPD, ' bytes'
write(LERR,*)' '
endif
c initialize memory
call vclr ( Record, 1, RecordSize )
call vclr ( Record_WorkSpace, 1, RecordSize )
call vclr ( Headers, 1, HeaderSize )
c BEGIN PROCESSING
c skip unwanted input 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 load record to memory
tr_index = 1 - nsamp
hdr_index = 1 - ITRWRD
DO KK = 1, ntrc
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 set array load points for this trace
tr_index = tr_index + nsamp
hdr_index = hdr_index + ITRWRD
c process only live traces and zero out dead traces. By convention, a "Station Correction"
c of 30000 is used as a dead trace flag!
call saver2 ( itr, ifmt_StaCor, l_StaCor, ln_StaCor, StaCor,
: TRACEHEADER )
if ( StaCor .ne. 30000 ) then
c load trace to array Record[]
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, Record(tr_index), 1, nsamp )
else
call vclr ( Record(tr_index), 1, nsamp )
endif
c load trace header to array Headers[]
call vmov ( itr, 1, Headers(hdr_index), 1, ITRWRD )
ENDDO
c Put your subroutine here [remember to declare any arguments you need
c over and above those already declared]
call subs ( Record, Headers, Record_WorkSpace, nsamp, ntrc,
: ist, iend, nsampo )
c reset array load points for this trace
tr_index = 1 - nsampo
hdr_index = 1 - ITRWRD
c write output data
DO KK = 1, ntrc
tr_index = tr_index + nsampo
hdr_index = hdr_index + ITRWRD
call vmov ( Record(tr_index), 1, itr(ITHWP1), 1, nsampo )
call vmov ( Headers(hdr_index), 1, itr(1), 1, ITRWRD )
call wrtape (luout, itr, obytes)
ENDDO
ENDDO
c close data files
call lbclos ( luin )
call lbclos ( luout )
write(LERR,*)'prgm: Normal Termination'
write(LER,*)'prgm: Normal Termination'
stop
999 continue
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
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, xman or xuspman'
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,*)'-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[] -rs[] -re[] -V'
write(LER,*)' '
write(LER,*)'===================================================='
return
end
c ----------------- Subroutine -----------------------
c pick up command line arguments
subroutine cmdln ( ntap, otap, irs, ire, ist, iend,
: name, verbos )
implicit none
#include
integer ist, iend, irs, ire, argis
character ntap*(*), otap*(*), name*(*)
logical verbos
call argi4 ( '-e', iend, -99999, -99999 )
call argstr ( '-N', ntap, ' ', ' ' )
call argstr ( '-O', otap, ' ', ' ' )
call argi4 ( '-re', ire, 0, 0 )
call argi4 ( '-rs', irs, 0, 0 )
call argi4 ( '-s', ist, 1, 1 )
verbos = (argis('-V') .gt. 0)
c check for extraneous arguments and abort if found to
c catch all manner of user typo's
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, verbos)
#include
integer nsamp, ntrc, iform, ist, iend, irs, ire, 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 = ', ire
write(LERR,*) ' processing sample start = ', ist
write(LERR,*) ' processing sample end = ', iend
if ( verbos ) write(LERR,*) ' verbose printout requested'
write(LERR,*)' '
write(LERR,*)'========================================== '
write(LERR,*)' '
return
end
