c ----------------- Main Routine -----------------------
c
c copyright 2001, 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
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
c itr[]is an integer buffer used to store a USP trace element
c SZLNHD is a "big number" that is "guaranteed to provide more than
c enough space to store a trace (including both trace header and
c trace data) or a line header"; it is defined in lhdrsz.h.
integer itr ( SZLNHD )
integer nsamp, nsi, ntrc, nrec, nreco, iform, obytes
integer luin , luout, lbytes, nbytes, lbyout
integer ist, iend, irs, ire, argis, jerr
real UnitSc, dt
character ntap*255, otap*255, name*4
logical verbos
c Program Specific _ dynamic memory variables
integer RecordSize, HeaderSize, errcd1, errcd2, errcd3, abort
integer Headers
real Record, Record_WorkSpace
pointer (mem_Record, Record(2))
pointer (mem_Space, Record_WorkSpace(2))
pointer (mem_Headers, Headers(2))
c Program Specific _ static memory variables
integer ifmt_StaCor,l_StaCor,ln_StaCor, StaCor
integer hdr_index, tr_index, JJ, KK, lenth, length
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 generally
c not useful for Fortran programmers.) If for your application you'll
c need to do a lot of fancy text processing or low-level (non-USP)
c I/O, consider writing your program in C instead of FORTRAN. (SIS
c I/O itself, for example, is 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 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.
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
length = lenth(ntap)
write(LERR,*)'PRGM: no line header on input dataset',
: ntap(1:length)
write(LER,*)' '
write(LER,*)'PRGM: '
write(LER,*)' no line header on input dataset',ntap(1:length)
write(LER,*)'FATAL'
write(LER,*)' '
stop
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
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
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
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)
c number output bytes
obytes = SZTRHD + SZSMPD * nsamp
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 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 used later
c in the program.
call savelu ( 'StaCor', ifmt_StaCor, l_StaCor, ln_StaCor,
: TRACEHEADER )
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 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 )
c reset array load points for this trace
tr_index = 1 - nsamp
hdr_index = 1 - ITRWRD
c write output data
DO KK = 1, ntrc
tr_index = tr_index + nsamp
hdr_index = hdr_index + ITRWRD
call vmov ( Record(tr_index), 1, itr(ITHWP1), 1, nsamp )
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[] -- process start time (ms) (1)'
write(LER,*)'-e[] -- process 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