@PROCESS DIRECTIVE('*VDIR:') C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CTITLESAMPFKB -- PERFORM THE MIGRATION CA AUTHOR BRUCE VER WEST CA DESIGNER BRUCE VER WEST CA LANGUAGE VS FORTRAN CA SYSTEM IBM (SEE CRAY) CA WRITTEN 07/01/86 C REVISED 12-04-87 BJV INCREASED BUFFER BLOCKING ARRAY C DIMENSIONS TO 39 AND 19. C REVISED 03-03-88 BJV CHANGE UNBUFFERED I/O TO USE UWOPEN C AND CALLS DIRECTLY TO READWA, WRITEWA C AND WUNIT INSTEAD OF USING BENCHLIB. C REVISED 08-25-88 ESN MINOR MODIFICATIONS OF 300 DO-LOOP FOR C IBM VECTORIZATION. C REVISED 02-20-89 JJC FOR SPARC PRODUCTION. C REVISED 03-24-89 JJC CHANGED IWPAN TO IWP. C REVISED 10-26-89 LWC INT FUNCTIONS FOR CFT77. C REVISED 05-09-91 CLJ BREAK UP DATA FILE INTO EIGHT PARTS C INSTEAD OF FOUR C CHANGE TO IBM ONLY BECAUSE LOOP 470 HAS C VECTOR DIRECTIVE (PREFER SCALAR) C BECAUSE OF COMPILER BUG WHEN VECTORIZED C CA SUBROUTINE SAMPFKB(VELD2,W,AKX,AKH,AKZSQ,IFAXX,TRIGSX, CA + WR,WI,P,Q,COEF,CNORM,LR,CSHIFT, CA + INDEX2,INDEX3,INDEX4,INDEX5,INDEX6, CA + INDEX7,INDEX8,INDEX9,INDEXA,ITHRSH,IPR,IABORT) C CA CA THIS SUBROUTINE PERFORMS THE DOWNWARD EXTRAPOLATION AND THE CA IMAGING FOR THE MIGRATION PROCESS USING THE INPUT VELOCITY. CA A SINC COMPLEX INTERPOLATION IS USED IN EVALUATING THE COMPLEX CA TEMPORAL FREQUENCY WAVE FIELD BEFORE INVERSE TRANSFORMING FROM CA TEMPORAL FREQUENCY TO TIME. CA CA********************************************************************** CA CA ARGUMENTS: CA CA VELD2 INPUT VELOCITY TO BE USED IN THE WAVE FIELD CA EXTRAPOLATION FOR ALL MID-POINTS TO CA BE PROCESSED CA CA W INPUT TEMPORAL RADIAL FREQUENCIES RESULTING CA FROM THE FOURIER TRANSFORM OF THE CA INPUT TIME TRACE CA CA AKX INPUT SPATIAL RADIAL FREQUENCIES RESULTING CA FROM THE FOURIER TRANSFORM OF THE CA MID-POINT COORDINATE COMPLEX FIELD CA CA AKH INPUT SPATIAL RADIAL FREQUENCIES RESULTING CA FROM THE FOURIER TRANSFORM OF THE CA OFFSET COORDINATE COMPLEX FIELD CA CA AKZSQ OUTPUT AN ARRAY USED TO STORE V**2/W**2 CA AS A FUNCTION OF TEMPORAL FREQUENCY CA W FOR THE EXTRAPOLATION (OR WEIGHTING CA FACTOR) FOR TIME MIGRATION CA CA IFAXX INPUT ARRAY WHICH CONTAINS THE PRIME FACTORS CA FOR THE GIVEN LENGTH OF THE SPATIAL CA FFT AS REQUIRED BY CFFT99. THE LENGTH CA IS NKX. CA CA TRIGSX INPUT TRIG TABLES FOR THE GIVEN LENGTH OF CA SPATIAL FFT AS REQUIRED BY CFFT99 CA CA WR OUTPUT REAL PART OF THE COMPLEX FIELD CA P(W,KX,KH) FOR A GIVEN OFFSET WAVE CA NUMBER KH AS A FUNCTION OF TEMPORAL CA FREQUENCY W. USED IN THE DOWNWARD CA CONTINUATION AND IMAGING. CA CA WI OUTPUT IMAGINARY PART OF THE COMPLEX FIELD CA P(W,KX,KH) FOR A GIVEN OFFSET WAVE CA NUMBER KH AS A FUNCTION OF TEMPORAL CA FREQUENCY W. USED IN THE DOWNWARD CA CONTINUATION AND IMAGING. CA CA P INPUT/ COMPLEX ARRAY WHICH CONTAINS THE CA OUTPUT INPUT/OUTPUT TO THE PARALLEL FFT CA ROUTINE CFFT99 FOR TRANSFORMING FROM CA MID-POINT COORDINATE TO FREQUENCY KM CA CA Q OUTPUT COMPLEX WORK ARRAY OF SAME LENGTH AS CA P WHICH IS REQUIRED BY CFFT99 CA CA COEF INPUT REAL COEFFICIENT ARRAY FOR THE SIX CA POINT SINC INTERPOLATOR FUNCTION FOR CA INTERPOLATING THE FIELD OVER TEMPORAL CA FREQUENCY. THIS IS A COMPLEX INTER- CA POLATION CA CA CNORM INPUT FREQUENCY DOMAIN COMPLEX SINC INTER- CA POLATION FILTER NORMALIZATION FACTOR. CA CA LR INPUT THE NUMBER OF POINTS IN THE SINC CA INTERPOLATION FILTER. CA CA CSHIFT FREQUENCY DOMAIN COMPLEX SINC INTER- CA POLATION FILTER TIME SHIFT FOR DE- CA MODULATION/MODULATION. CA CA INDEX2 INPUT DATA SET NAME WHICH CONTAINS THE CA COMPLEX IMAGE FIELD (INTEGRATED CA OVER KH) AS A FUNCTION OF TEMPORAL CA FREQUENCY W AND MID-POINT SPATIAL CA WAVE NUMBER KX; P(W,KX) CA CA INDEX3 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEX4 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEX5 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEX6 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEX7 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEX8 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEX9 INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA INDEXA INPUT DATA SET FILE NAME WHICH CONTAINS CA A BLOCK OF TRACES FOR INPUT TO THE CA FFT ROUTINE CFFT99 AS A FUNCTION CA OF OFFSET FREQUENCY; P(W,KX,KH) CA CA ITHRSH INPUT HARDWIRED VALUE OF THREE WHICH SERVES CA AS A LOWER THRESHOLD FOR WEIGHTING THE CA IMAGE BY THE ACTUAL FOLD. CA CA IPR INPUT SPARC LOGICAL UNIT NUMBER FOR PRINT CA CA IABORT INPUT/ FLAG WHICH IS RETURNED TO SPARC TO CA OUTPUT DETERMINE IF PROCESSING IS TO CONTIN- CA UE CA CA********************************************************************** C SUBROUTINE SAMPFKB(VELD2,W,AKX,AKH,AKZSQ,IFAXX,TRIGSX, + WR,WI,P,Q,COEF,CNORM,LR,CSHIFT, + INDEX2,INDEX3,INDEX4,INDEX5,INDEX6, + INDEX7,INDEX8,INDEX9,INDEXA,ITHRSH,IPR,IABORT) C COMMON/CMPFKC/IF1,LNT,LW,LW2,IW1,ALPHA,SCALE,CMIN,VMUTE,AFFR,IKHHI COMMON/CMPFKT/NT,DT,NW,DW,NWD2,NWD21,NWP2 COMMON/CMPFKX/NX,DX,NKX,DKX,NKXD2,NKXD21,NKX2,NKXP2 COMMON/CMPFKH/NH,DH,NKH,DKH,NKHD2,NKHD21,NKH2,IHBEG COMMON/CMPFKB/IKXBF(39),IWBF(19),MKXBF(39),MWBF(19), + NBF,MBF,NKXBF,NWBF C COMPLEX P(LNT/2,NKXBF,1),Q(LNT/2,1),CPHASE,CNORM(101),CC,CSHIFT(*) DIMENSION W(1),AKH(1),AKZSQ(1),WR(1),WI(1),COEF(LR,101) DIMENSION AKX(1),IFAXX(1),TRIGSX(1) C C C C LR = 6 ONLY !!!!!!!!!!!!!!!!!!!!! C CALL WUNIT(INDEX3) CALL WUNIT(INDEX4) CALL WUNIT(INDEX5) CALL WUNIT(INDEX6) CLJ1 CALL WUNIT(INDEX7) CALL WUNIT(INDEX8) CALL WUNIT(INDEX9) CALL WUNIT(INDEXA) CLJ2 C IKHHI4=(IKHHI+3)/4 TPI = 2.*3.1415926536 ALW1 = FLOAT(LW+1) ALW = FLOAT(LW) WLMAX= ALW1*DW CALL ARSET(WR(1),LW+15,0.0) CALL ARSET(WI(1),LW+15,0.0) C DO 100 IWP=2,LW AKZL=W(IWP)/VELD2 100 AKZSQ(IWP)=1./(AKZL*AKZL) C C DETERMINE HIGHEST KH VALUE TO PROCESS. PFACT IS THE FRACTION OF C MAXIMUM TIME DIP 1/V WHICH IS ALLOWED FOR AT MAXIMUM FREQUENCY. C USUALLY THE MAXIMUM TIME DIP IN A GATHER IS MUCH LESS THAT 1/V C DUE TO MAXIMUM OFFSET AND MUTING. C PFACT = 0.70*SQRT(2.0/(1.0+(2.0*VELD2/CMIN)**2)) IF(IKHHI.EQ.1) THEN IKHMAX=1 ELSE C IKHMAX = 1+INT(2.0*PFACT*DW*LW/(2.*VELD2*DKH)) IKHMAX = 1+INT((2.0*PFACT*DW*LW + .0001) / (2.*VELD2*DKH)) IKHMAX = MIN0(IKHMAX,IKHHI) ENDIF VEL = 2.0*VELD2 WRITE (IPR,8112) VEL,IKHMAX,PFACT 8112 FORMAT(/,' VELOCITY =',F7.1,' IHKMAX =',I6,' PFACT =',F6.3) C J2=0 DO 450 JBF=1,NBF KSEQDA=J2*LNT+1 ISEQDA=J2*LNT+1 J1=J2+1 J2=MIN0(IKXBF(JBF),NKX) JF1=LNT*MKXBF(JBF) C C TRIPLE BUFFERING OF INPUT FOR IMAGING STEP C C ICBF = CURRENT BUFFER (1,2, 3) MPC = CURRENT FILE (0,1,2,3) C INBF = NEXT BUFFER (1,2, 3) MP = NEXT FILE (0,1,2,3) C INNB = NEXT+ BUFFER (1,2, 3) MPP = NEXT+ FILE (0,1,2,3) C C READ FIRST AND SECOND BUFFER C LP = 0 MP = 0 INBF = 1 JSEQDA = LP*IF1+ISEQDA CALL READWA(INDEX3,P(1,1,INBF),JSEQDA,JF1,1) C LP = 0 MPP = 1 INNB = 2 IF (IKHMAX .GT. 1) THEN JSEQDA = LP*IF1+ISEQDA CALL READWA(INDEX4,P(1,1,INNB),JSEQDA,JF1,1) ENDIF C DO 500 IKH=1,IKHMAX AKHL=AKH(IKH) BKHL=ABS(AKHL) AKH2L=AKHL*AKHL ICBF = INBF MPC = MP INBF = INNB MP = MPP INNB = INNB + 1 IF (INNB .GT. 3) THEN INNB = 1 ENDIF C IF (IKH .LT. IKHMAX-1) THEN C C START READ OF NEXT BUFFER C IKHP = IKH+2 CCLJ* LP = (IKHP-1)/4 CCLJ* MPP = MOD(IKHP-1,4) LP = (IKHP-1)/8 MPP = MOD(IKHP-1,8) JSEQDA = LP*IF1+ISEQDA C IF(MPP .EQ. 0) THEN CALL READWA(INDEX3,P(1,1,INNB),JSEQDA,JF1,1) ELSE IF(MPP .EQ. 1) THEN CALL READWA(INDEX4,P(1,1,INNB),JSEQDA,JF1,1) ELSE IF(MPP .EQ. 2) THEN CALL READWA(INDEX5,P(1,1,INNB),JSEQDA,JF1,1) CLJ1 CCLJ* ELSE CCLJ* CALL READWA(INDEX6,P(1,1,INNB),JSEQDA,JF1,1) ELSE IF(MPP .EQ. 3) THEN CALL READWA(INDEX6,P(1,1,INNB),JSEQDA,JF1,1) ELSE IF(MPP .EQ. 4) THEN CALL READWA(INDEX7,P(1,1,INNB),JSEQDA,JF1,1) ELSE IF(MPP .EQ. 5) THEN CALL READWA(INDEX8,P(1,1,INNB),JSEQDA,JF1,1) ELSE IF(MPP .EQ. 6) THEN CALL READWA(INDEX9,P(1,1,INNB),JSEQDA,JF1,1) ELSE CALL READWA(INDEXA,P(1,1,INNB),JSEQDA,JF1,1) CLJ1 ENDIF ENDIF C C IF(IKH.EQ.1) THEN CALL WUNIT(INDEX2) CALL ARSET(Q,JF1,0) ENDIF C C CHECK IF READ OF CURRENT BUFFER IS COMPLETE C IF(MPC .EQ. 0) THEN CALL WUNIT(INDEX3) ELSE IF(MPC .EQ. 1) THEN CALL WUNIT(INDEX4) ELSE IF(MPC .EQ. 2) THEN CALL WUNIT(INDEX5) CLJ1 CCLJ* ELSE CCLJ* CALL WUNIT(INDEX6) ELSE IF(MPC .EQ. 3) THEN CALL WUNIT(INDEX6) ELSE IF(MPC .EQ. 4) THEN CALL WUNIT(INDEX7) ELSE IF(MPC .EQ. 5) THEN CALL WUNIT(INDEX8) ELSE IF(MPC .EQ. 6) THEN CALL WUNIT(INDEX9) ELSE CALL WUNIT(INDEXA) ENDIF C JI=0 DO 400 IKX=J1,J2 JI=JI+1 AKXL=AKX(IKX) BKXL=ABS(AKXL) AKX2L=AKXL*AKXL C AKHX2=(VELD2*VELD2*AKHL*AKXL)**2 C IND = LR/2-1 DO 200 IW=1,LW WR(IND+IW)= REAL(P(IW,JI,ICBF)) 200 WI(IND+IW)=AIMAG(P(IW,JI,ICBF)) C C IWPLO=INT(VELD2*SQRT(BKXL*BKHL)/DW)+2 IWPLO = INT((VELD2*SQRT(BKXL*BKHL) + .0001) / DW) + 2 C C IF POST-STACK (NH .EQ. 1) THEN USE P -> P TRANSFORM C IF PRE-STACK (NH .GT. 1) THEN USE P -> M TRANSFORM C IF (NH .EQ. 1) THEN DO 290 IWP=IWPLO,LW WL=W(IWP)*SQRT((1.+AKX2L*AKZSQ(IWP))*(1.+AKH2L*AKZSQ(IWP))) C SINC INTERPOLATION WLM = AMIN1(WL,WLMAX) C JW=INT(WLM/DW) JW = INT((WLM + .0001) / DW) IND = INT((WLM/DW-JW)*100.+1.5) C --- P --- WGHT=W(IWP)*(1.0-AKHX2/W(IWP)**4)/WL C SINC INTERPOLATION N=6 AWR= (COEF( 1,IND)*WR(JW+1) + COEF( 2,IND)*WR(JW+2) * +COEF( 3,IND)*WR(JW+3) + COEF( 4,IND)*WR(JW+4) * +COEF( 5,IND)*WR(JW+5) + COEF( 6,IND)*WR(JW+6))*WGHT AWI= (COEF( 1,IND)*WI(JW+1) + COEF( 2,IND)*WI(JW+2) * +COEF( 3,IND)*WI(JW+3) + COEF( 4,IND)*WI(JW+4) * +COEF( 5,IND)*WI(JW+5) + COEF( 6,IND)*WI(JW+6))*WGHT Q(IWP,JI)=Q(IWP,JI) * +CNORM(IND)*CSHIFT(JW+1)*CMPLX(AWR,AWI) 290 CONTINUE ELSE CESN C CHANGED IWP IWPA FOR LOOP 300 CESN DO 300 IWPA=IWPLO,LW WL=W(IWPA)*SQRT((1.+AKX2L*AKZSQ(IWPA))*(1.+AKH2L*AKZSQ(IWPA))) C SINC INTERPOLATION CESN WLM = AMIN1(WL,WLMAX) WLM = WL IF(WL.GT.WLMAX) WLM = WLMAX CESN C JW=INT(WLM/DW) JW = INT((WLM + .0001) / DW) IND = (WLM/DW-JW)*100.+1.5 C --- M --- C WGHT= (1.0-AKHX2/W(IWPA)**4)/SQRT(2.0*WL) WGHT= W(IWPA) * (1.0-AKHX2/W(IWPA)**4) / SQRT(2.0*WL) C WGHT=INT(W(IWPA)*(1.0-AKHX2/W(IWPA)**4)/SQRT(2.0*WL)) C SINC INTERPOLATION N=6 C AWR= (COEF( 1,IND)*WR(JW+1) + COEF( 2,IND)*WR(JW+2) C * +COEF( 3,IND)*WR(JW+3) + COEF( 4,IND)*WR(JW+4) C * +COEF( 5,IND)*WR(JW+5) + COEF( 6,IND)*WR(JW+6))*WGHT C AWI= (COEF( 1,IND)*WI(JW+1) + COEF( 2,IND)*WI(JW+2) C * +COEF( 3,IND)*WI(JW+3) + COEF( 4,IND)*WI(JW+4) C * +COEF( 5,IND)*WI(JW+5) + COEF( 6,IND)*WI(JW+6))*WGHT AWR= COEF( 1,IND)*WR(JW+1) + COEF( 2,IND)*WR(JW+2) AWI= COEF( 1,IND)*WI(JW+1) + COEF( 2,IND)*WI(JW+2) AWR=AWR +COEF( 3,IND)*WR(JW+3) + COEF( 4,IND)*WR(JW+4) AWI=AWI +COEF( 3,IND)*WI(JW+3) + COEF( 4,IND)*WI(JW+4) AWR=AWR +COEF( 5,IND)*WR(JW+5) + COEF( 6,IND)*WR(JW+6) AWI=AWI +COEF( 5,IND)*WI(JW+5) + COEF( 6,IND)*WI(JW+6) Q(IWPA,JI)=Q(IWPA,JI) * +CNORM(IND)*CSHIFT(JW+1) * *CMPLX(WGHT*(AWR+AWI),WGHT*(AWI-AWR)) 300 CONTINUE ENDIF C 400 CONTINUE 500 CONTINUE C C WEIGHT BY ACTUAL FOLD C THRSH = FLOAT(ITHRSH) IF (NH .GT. 1) THEN SNORM2 = 1.0/(VELD2*DKH*SQRT((2.0*VELD2/VMUTE)**2-1.0)) JI = 0 DO 480 IKX=J1,J2 JI = JI+1 BKXL = ABS(AKX(IKX)) IF (BKXL .EQ. 0.0) THEN SNORM = FLOAT(IKHMAX)/DW**2*1.01 ELSE SNORM = 1.0/(VELD2**2*DKH*BKXL) ENDIF C*VDIR: PREFER SCALAR DO 470 IWP = 2,LW ANP = AINT(1.0+SNORM*W(IWP)**2) AN2 = AINT(1.0+SNORM2*W(IWP)) ANP = AMIN1(ANP,AN2,FLOAT(IKHMAX)) ANP = AMAX1(ANP,THRSH) WGHT = 1.0/(ANP*2.0*VELD2) Q(IWP,JI) = Q(IWP,JI)*WGHT 470 CONTINUE 480 CONTINUE ENDIF C IF(NBF.GT.1) THEN CALL WUNIT(INDEX2) CALL WRITEWA(INDEX2,Q,KSEQDA,JF1,1) ENDIF 450 CONTINUE C IF(NX.EQ.1) RETURN C IF(NBF.EQ.1) THEN CALL CFFT99(Q,P,TRIGSX,IFAXX,LNT/2,1,NKX,LW,+1) ELSE CALL SAMPFKG(IFAXX,TRIGSX,P,Q,INDEX2,IPR) ENDIF C RETURN END