As an example of SEPlib90 data access, below is presented a short program called Filter which performs smoothing via triangle filtering. Though the code is written in Fortran90, it uses almost exclusively Fortran77 functions, and but for a few lines could be compiled with f77. The only features unfamiliar to Fortran77 programmers should be the array declarations of variables trace and grid, the presence of an optional subroutine argument, and the use of a module.
Aside from variable declarations, there are really just two parts into this program. The first is a series of error checking statements, each of which calls the program-halting subroutine seperr if some condition is not met. SEPlib90-specific subroutines are identifiable in the code by their naming convention, using multi-word phrases beginning with "sep". The first of these, sep_get_number_grid_axes, returns a zero value if successful, and some non-zero integer value if it fails. The value returned can be used to provide some diagnosis of the problem if one exists, but in this case the program simply aborts any time that the value is non-zero, with the message ``Trouble getting grid axes.'' Assuming that no problems occur, program Filter stores the number of grid axes (obtained by sep_get_number_grid_axes from the grid format file) in the variable num_axes. There is then a further check that aborts the program if num_axes falls outside the range (2,3,4,5,6). A quick glance at the main loop structure reveals why this is done.
Notice that nowhere are the familiar SEPlib hetch or from history statements that are used to get information on the input data's dimensionality. This job is instead carried out by the subroutines sep_get_data_axis_par and sep_get_grid_axis_par, which appear next in the code. The length and sampling of the time axis is determined by the former. Similar information for whatever spatial (or other) axes exist in the data are determined by the latter. This assumes that the data are regular according to the gridding coordinate system determined in the grid format file. This would be a foolhardy assumption for most processes, where spatial information is important and would be gotten from the headers. However, Filter works solely on the time axis, and is a nonphysical operation, such that precise information about source and receiver coordinates, etc. is unnecessary.
The final error checking is for the existence the header key, data_record_number. Because SEPlib90 header and data are stored separately, and may be manipulated independently, the possibility exists that the headers will be ordered differently on disk than the data (and there may even be fewer headers than traces). The header key data_record_number serves as a pointer from a given header to a trace. Program Filter checks for the existence of a header called data_record_number, and if it finds such a header key, it assumes that headers and data are not stored in the same order. To deal with this, the header key index of data_record_number is stored as the variable data_key, and the logical variable in_order is set to false.
The remainder of the program is fairly simple. Two arrays are allocated. One, named trace exists to temporarily store each input trace for filtering. The other, named grid, is an integer vector containing a subset of the header record numbers in the input data. Also, nwind, jwind, and fwind, parameters to the subroutine that gets header record numbers from the grid (sep_get_grid_window), are set to default values. Finally, each axis of the grid is looped over. Header record numbers are grabbed a vector at a time from the grid and stored in the array grid. The array grid is then looped over, each trace being input, filtered, and output.
program Filter
use filtmods
implicit none
real ,dimension(:),pointer :: trace
integer,dimension(:),pointer :: grid
integer :: n(6)
real :: d(6),o(6)
integer :: i1,i2,i3,i4,i5,i6
integer :: sep_get_number_grid_axes
integer :: sep_get_grid_axis_par
integer :: getch,sep_get_data_axis_par,fetch
integer :: sep_get_key_index,putch,esize
integer :: sep_get_grid_window
integer :: same_record_number,num_axes
integer :: length,data_key
integer :: fwind(6),jwind(6),nwind(6),ierr
logical :: in_order
character(len=128) :: label(6)
call initpar()
if(0==fetch("esize","d",esize)) then
esize=4
end if
if(0.ne.putch("esize","d",esize)) then
call seperr("trouble writing to output history file")
end if
if(esize .ne. 4) then
call seperr("this program can only deal with real data")
end if
if(0 .ne. sep_get_number_grid_axes("in",num_axes)) then
call seperr("Trouble getting grid axes")
end if
if(num_axes >6 .or. num_axes <2 ) then
call seperr("Not equipped to deal with given number of axes")
end if
n=1
if(0 .ne. sep_get_data_axis_par("in",1,n(1),o(1),d(1),label(1))) then
call seperr("Trouble getting data axis information")
end if
do i1=2,num_axes
if(0 .ne. sep_get_grid_axis_par("in",i1,n(i1),o(i1),d(i1), &
label(i1))) then
call seperr("Trouble getting grid axis information")
end if
end do
if(0==getch("length","d",length)) then
length=4
end if
if(0==fetch("same_record_number","d",same_record_number)) then
same_record_number=1
end if
if(same_record_number == 0) then
in_order=.false.
ierr=sep_get_key_index("in","data_record_number",data_key)
if(ierr .ne. 0) then
call seperr("Data is out of order and cannot find &
data record number key")
end if
else
in_order=.true.
end if
allocate(grid(1:n(2)))
allocate(trace(1:n(1)))
jwind=1
nwind(2:num_axes-1)=1
fwind(1) = 0
nwind(1)=n(2)
do i6=1,n(6)
fwind(5) = i6 -1
do i5=1,n(5)
fwind(4) = i5 -1
do i4=1,n(4)
fwind(3) = i4 -1
do i3=1,n(3)
fwind(2) = i3 -1
call sep_get_grid_window("in",num_axes,n(2),nwind,fwind,jwind,grid)
do i2=1,n(2)
if(grid(i2) .ne. -1 ) then
call get_trace("in",in_order,n(1),esize,grid(i2), &
trace,data_key=data_key)
call filter_trace(trace,n(1),length)
call put_trace("out",n(1),esize,grid(i2),trace)
end if
end do
end do
end do
end do
end do
!if(0 .ne. putch("out",same_record_number,1)) then
if(0 .ne. putch("same_record_number",'d',1)) then
same_record_number=1
end if
end program Filter
Four subroutines are called within the inner loop. Subroutine filter_trace is simply convolution and will be skipped over here. Subroutine get_trace is a short routine which checks to see if the data and headers are stored in the same order, and acts accordingly to call sseek and sreed. Subroutine put_trace handles the output in a straightforward manner. Notice that even in the case that the input headers and data appear in different orders, the output will be in the same order.
Finally, subroutine sep_get_grid_window is the least obvious in its usage. It is complicated by the fact that the one axis (time axis) is implicit in a SEPlib90 grid. A data set gridded according to source x and source y header values only is actually described as having three grid axes. A look into the grid format file would show that source x and source y were the two and three axes, respectively. The one axis is not described in the grid format file. What this means is that the num_axes argument to the sep_get_grid_window subroutine is an integer whose value is one greater than the length of the vectors nwind, jwind, fwind, and n.
|
blocks
Figure 1 Input to the program Filter.x |  |
|
temp
Figure 2 Output of the program Filter.x |  |
module filtmods
implicit none
contains
subroutine get_trace(tag,in_order,n,esize,record_number,trace,data_key)
implicit none
character(len=*) :: tag
logical :: in_order
integer :: n
integer :: esize
real :: trace(n)
integer, optional :: data_key
integer :: record_number,sep_get_val_by_index
integer :: data_record_number,ierr
integer :: sreed,sseek
integer*8 :: sseek_pos
real :: tempa(512)
if(.not. in_order) then
if (.not. present(data_key)) then
call seperr("must supply a data key if the traces &
are out of order")
end if
end if
if(in_order) then
sseek_pos = (record_number -1) * esize * n
else
if(0 .ne. sep_get_val_by_index(tag,record_number,data_key, &
1, data_record_number)) then
call seperr("trouble getting data_record_number ")
end if
sseek_pos = (data_record_number -1) * esize * n
end if
ierr=sseek(tag,sseek_pos,0)
ierr=sreed(tag,trace,n*esize)
end subroutine get_trace
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine put_trace(tag,n,esize,record_number,trace)
character(len=*) :: tag
integer :: n
integer :: esize
integer :: record_number,sep_get_val_by_index
real :: trace(n)
integer :: data_record_number
integer*8 :: sseek_pos
sseek_pos = (record_number -1) * esize * n
call sseek(tag,sseek_pos,0)
call srite(tag,trace,n*esize)
end subroutine put_trace
end module filtmods