The following code fragments express in Fortran 90 a Kirchhoff prestack time-migration algorithm. In this example the data are assumed to be sorted in common-offset sections and the traveltimes are computed analytically assuming an RMS slowness function. Furthermore, for the sake of simplicity, I did not include the constant phase shift that should be applied to the results of Kirchoff migration.
First I show how the desired data layout is specified to the Fortran compiler.
C C arrays declarations C REAL DATAL(N1,N2) REAL DATAR(N1,N2) REAL MIGRES(NTAU,N2)
C C arrays layout C C :SERIAL = local axis C :NEWS = parallel axis C
CMF$ LAYOUT DATAL(:SERIAL,:NEWS) CMF$ LAYOUT DATAR(:SERIAL,:NEWS) CMF$ LAYOUT MIGRES(:SERIAL,:NEWS)
Then, I show the main loop of the algorithm.
C C initialize buffers for right-side half of operator C DATAR=DATAL
C C loop over operator lag C DO IOPER=1,LOPER
C C compute constants on the front-end C DELMIDPOINT=(IOPER)*D2 DELMH2=(DELMIDPOINT-HOFFSET)**2 DELPH2=(DELMIDPOINT+HOFFSET)**2
C C shift input data to the left and to the right C DATAL=CSHIFT(DATAL,DIM=2,SHIFT=1) DATAR=CSHIFT(DATAR,DIM=2,SHIFT=-1)
C C loop over output times C DO ITAU=3,NTAU
C C compute traveltimes on the front-end C TIME=SLOWTAU(ITAU) * 1 (SQRT(DEPTH2(ITAU)+(DELMH2)) + 2 SQRT(DEPTH2(ITAU)+(DELPH2))) TIME=((TIME-O1)*INVD1)+1 ITIME=INT(TIME)
C C compute amplitudes (your favorite amplitude function can go here) C AMP=1./TIME
IF (ITIME .LE. (N1-3)) THEN
C C compute interpolation samples and weights on the front-end, C using the interpolator table TABINT C ITAB=INT((TIME-ITIME)*LTABINT)+1
INDEX1=ITIME INDEX2=ITIME+1 INDEX3=ITIME-1 INDEX4=ITIME+2 INDEX5=ITIME-2 INDEX6=ITIME+3
INTER1=AMP*TABINT(0,ITAB) INTER2=AMP*TABINT(-1,ITAB) iNTER3=AMP*TABINT(1,ITAB) INTER4=AMP*TABINT(-2,ITAB) INTER5=AMP*TABINT(2,ITAB) INTER6=AMP*TABINT(-3,ITAB)
C C compute contributions to the output and sum into output array C
MIGRES(ITAU,:)=MIGRES(ITAU,:) 1 + INTER1*DATAL(INDEX1,:) 2 + INTER1*DATAR(INDEX1,:) 3 + INTER2*DATAL(INDEX2,:) 4 + INTER2*DATAR(INDEX2,:) 5 + INTER3*DATAL(INDEX3,:) 6 + INTER3*DATAR(INDEX3,:) 7 + INTER4*DATAL(INDEX4,:) 8 + INTER4*DATAR(INDEX4,:) 9 + INTER5*DATAL(INDEX5,:) 1 + INTER5*DATAR(INDEX5,:) 2 + INTER6*DATAL(INDEX6,:) 3 + INTER6*DATAR(INDEX6,:)
In this simple case each processor does exactly the same operation because the traveltimes are independent of the midpoint location. The same would be true if dip moveout were implemented instead of migration. On the contrary, in a depth migration algorithm the impulse response changes with the midpoint location, and thus each processor needs to perform a slightly different operation. More precisely, the indices of input time samples used in the summation differ from processor to processor. This requires the use of the indirect addressing capabilities of the processing hardware, and it usually leads to an increase in run time. Nevertheless, Kirchhoff depth migration can be efficiently run on a massively parallel computer (Van Trier, 1991). If the data are not evenly sampled, as in 3-D surveys, implementing a Kirchhoff migration algorithm would become even more complicated. However, also in this case, the basic idea of exploiting the spatial continuity of the impulse response can be used for deriving efficient algorithms.