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.
Figure 1 Input to the program Filter.x
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