CTITLESAZM06B-- PERFORMS DEPTH MIGRATION OPERATION FOR ZMIG (DIP.GT.65) 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA 00020000 CA AUTHOR FRANCIS SHERRILL--AOGC R&D 00030000 CA DESIGNER FRANCIS SHERRILL--AOGC R&D 00040000 CA LANGUAGE FORTRAN 00050000 CA SYSTEM IBM (SEE CRAY) 00060000 C WRITTEN MARCH 1980 00070000 C REVISED NOV,1980 BY L.B. LIN 00080000 C REVISED DEC,1980 BY F.G.SHERRILL 00090000 C REVISED APR,1981 BY R.D.KNIGHT 00100000 C REVISED APR,1986 BY R.D.KNIGHT; MOVED TO CRAY. 00110000 C REVISED APR,1986 BY R.D.KNIGHT; MOVED TO CRAY. 00120000 C REVISED OCT,1986 BY F.A.COADY ; MOVED VEL FIELD FROM 1ST WORK 00130000 C FILE TO BLANK COM ARRAY (V). 00140000 C REVISED JAN,1988 BY JAMES SUN 00141004 C (1) ADDED 80 DEGREE OPTION 00142004 C (2) CORRECTED MISTAKE IN THE C-ARRAY 00143004 CA 00150000 CA CALL SAZM06B(BETA,BUF1,BUF2, 00160000 CA DOM,DT,DX,DZ,DZ2,DZ3,ETA,IFLAG,IPASS,IPR, 00170000 CA ISTO,IUP,IVCHK,IVDATM, V ,KPWRK2,NFB,NFCPB, 00180000 CA NOM,NOMCHK,NOMLO,NOMLST,NTB,NTPB,NTT,NX, 00190000 CA NZ,NZE,NZS,REF,S,VEL,WAVI,WAVR, WAVE ) 00200000 CA 00210000 CA IN/OUT ARG TYPE DESCRIPTION 00220000 CA 00230000 CA IN BETA R4 VARIABLE USED IN SECOND DERIVATIVE APPROX. 00240000 CA SCR BUF1 R4 ARRAY USED FOR BUFFERING DATA IN FROM DISK 00250000 CA SCR BUF2 R4 ARRAY USED FOR BUFFERING DATA IN FROM DISK 00260000 CA IN DOM R4 RADIAL FREQUENCY SPACING 00270000 CA IN DT R4 SAMPLE RATE OF DATA IN MILS 00280000 CA IN DX R4 TRACE SPACING IN FEET 00290000 CA IN DZ R4 DEPTH INCREMENT IN FEET 00300000 CA IN DZ2 R4 MAXIMUM DEPTH 00310000 CA IN DZ3 R4 ORIGINAL DEPTH INCREMENT USED IN DATUMING 00320000 CA (EXCEPT FOR FINAL Z-STEP) 00330000 CA IN ETA R4 DIP FILTERING COEFFICIENT 00340000 CA IN IFLAG I4 SWITCH TO INDICATE MIGRATION OR DATUMING 00350000 CA IN IPASS I4 USED IN DATUMING OPTION WHEN THE UPWARD 00360000 CA EXTRAPOLATION IS SPECIFIED. EQUALS 1 ON 00370000 CA A DOWNWARD MIGRATION STEP, AND 2 ON THE 00380000 CA MODELING STEP 00390000 CA IN IPR I4 UNIT NUMBER OF OUTPUT PRINT FILE 00400000 CA IN ISTO I4 ARRAY CONTAINING RANGES WITH UNIFORM VELOCITY 00410000 CA IN IUP I4 SWITCH TO INDICATE WHETHER WAVEFIELD IS TO BE 00420000 CA CARRIED BACK TO SURFACE DURING DATUMING 00430000 CA IN IVCHK I4 ARRAY INDICATING ADJACENT LAYERS WITH EQUAL 00440000 CA VELOCITY 00450000 CA IN IVDATM I4 REPLACEMENT VELOCITY FOR USE IN DATUMING 00460000 CA IN V I4 VELOCITY FIELD ARRAY DIMENSIONED (NX X 3*NZ) 00470000 CA IN KPWRK2 I4 DCB ADDRESS OF WORKFILE 2 00480000 CA IN NFB I4 NUMBER OF FREQUENCY BLOCKS 00490000 CA IN NFCPB I4 NUMBER OF FREQUENCY COMPONENTS PER BLOCK 00500000 CA IN NOM I4 TOTAL NUMBER OF FREQUENCY COMPONENTS TO BE 00510000 CA DOWNWARD CONTINUED 00520000 CA IN NOMCHK I4 ARRAY OF INDICATORS FOR DOWNWARD CONTINUATION 00530000 CA OF SPECTRAL FREQUENCIES 00540000 CA IN NOMLO I4 LOWEST FREQUENCY COMPONENT TO MIGRATE 00550000 CA IN NOMLST I4 LAST FREQUENCY COMPONENT TO MIGRATE 00560000 CA IN NTB I4 NUMBER OF TRACE BLOCKS 00570000 CA IN NTPB I4 NUMBER OF TRACES PER BLOCK 00580000 CA IN NTT I4 NUMBER OF SAMPLES IN FFT-AUGMENTED TRACELENGTH00590000 CA IN NX I4 TOTAL TRACES TO MIGRATE 00600000 CA IN NZ I4 NUMBER OF DOWNWARD CONTINUATION DEPTH STEPS 00610000 CA IN NZE I4 ENDING DEPTH STEP NUMBER FOR CURRENT PASS 00620000 CA IN NZS I4 START DEPTH STEP NUMBER FOR CURRENT PASS 00630000 CA SCR REF R4 DEPTH DOMAIN MIGRATED WAVEFIELD BUFFER 00640000 CA SCR S R4 CONTIGUOUS SCRATCH ARRAYS FROM BLANK COMMON 00650000 CA IN VEL R4 VELOCITY FIELD 00660000 CA SCR WAVI R4 ARRAY CONTAINING IMAGINARY PART OF WAVEFIELD 00670000 CA SCR WAVR R4 ARRAY CONTAINING REAL PART OF WAVEFIELD 00680000 CA SCR WAVE R4 DUMMY ARRAY -- NOT REFERENCED IN IBM VERSION 00690000 CA 00700000 CA 00710000 CA THIS SUBROUTINE PERFORMS THE DOWNWARD CONTINUATION OF ARRAY WAVE. 00720000 CA FREQUENCY COMPONENTS ARE CARRIED DOWN IN A MULTI-STEP FASHION. 00730000 CA THIS ALLOWS CERTAIN COMPUTATIONAL ADVANTAGES IF THE VELOCITY IS 00740000 CA CONSTANT OVER ADJACENT DEPTH LEVELS. THIS IS A DEPTH MIGRATION 00750000 CA ALGORITHM IN THAT THE DOWNARD EXTRAPOLATION IS IN INCREMENTS OF 00760000 CA DEPTH. THE INPUT WAVEFIELD MUST BE IN THE SPACE-FREQUENCY DOMAIN. 00770000 CA 00780000 CA 00790000 CMEMORY 00800000 SUBROUTINE SAZM06B(BETA,BUF1,BUF2, 00810002 * DOM,DT,DX,DZ,DZ2,DZ3,ETA,IFLAG,IPASS,IPR, 00820000 * ISTO,IUP,IVCHK,IVDATM, V ,KPWRK2,NFB,NFCPB, 00830000 * NOM,NOMCHK,NOMLO,NOMLST,NTB,NTPB,NTT,NX, 00840000 * NZ,NZE,NZS,REF,S,VEL,WAVI,WAVR, WAVE ) 00850000 C 00860000 C COMPLEX VARIABLES--LOCAL AND EXTERNAL 00870000 C 00880000 COMPLEX BBNX 00890000 COMPLEX BB1 00900000 COMPLEX BNX 00910000 COMPLEX BCS 00920000 COMPLEX B1 00930000 COMPLEX CMPLX 00940000 COMPLEX VIS 00950000 COMPLEX Z1 00960000 COMPLEX Z2 00970000 COMPLEX Z3 00980000 COMPLEX Z4 00990000 COMPLEX Z5 01000000 C 01010000 C REAL ARRAYS -- PARAMETER LIST 01020000 C 01030000 REAL BUF1 ( 1) 01040000 REAL BUF2 ( 1) 01050000 REAL REF ( NX) 01060000 REAL S ( 1) 01070000 REAL VEL ( NX, 2) 01080000 REAL WAVI ( NX, NFCPB) 01090000 REAL WAVR ( NX, NFCPB) 01100000 REAL WAVE ( 1) 01110000 CMEMORY 01120000 REAL V ( NX, 3*NZ) 01130000 C 01140000 C INTEGER ARRAYS--PARAMETER LIST 01150000 C 01160000 DIMENSION ISTO ( NX) 01170000 DIMENSION IVCHK ( NZ) 01180000 DIMENSION NOMCHK( NOM) 01190000 C 01200000 C INTEGER VARIABLES--LOCAL 01210000 C 01220000 INTEGER DAWRK 01230000 INTEGER DAWRK2 01240000 C 01250000 C INTEGER VARIABLES--LOCAL 01260000 C 01270000 INTEGER AAI 01280000 INTEGER AAR 01290000 INTEGER AI 01300000 INTEGER AR 01310000 INTEGER BBI 01320000 INTEGER BBR 01330000 INTEGER BI 01340000 INTEGER BR 01350000 INTEGER DENDEN 01360000 INTEGER DENI 01370000 INTEGER DENR 01380000 INTEGER DI 01390000 INTEGER DR 01400000 INTEGER EI 01410000 INTEGER ER 01420000 INTEGER FI 01430000 INTEGER FR 01440000 INTEGER SHIFTI 01450000 INTEGER SHIFTR 01460000 INTEGER TI 01470000 INTEGER TR 01480000 C 01490000 C CONSTANTS--LOCAL 01500000 C 01510000 DATA ID/'D'/ 01520000 DATA IM/'M'/ 01530000 DATA Q1INA/0.9143/ 01551000 DATA Q2INA/0.8740/ 01552000 DATA Q1INB/0.6800/ 01553000 DATA Q2INB/0.2227/ 01554000 DATA DW/0.00001/ 01560000 C 01570000 C ******************************************************************** 01580000 C GENERATE PARAMETERS AND SET CONSTANTS: 01590000 C RR1,RR3,----REAL AND COMPLEX COEFFICIENTS USED IN COMPUTING THE 01600000 C DIAGONAL ARRAYS IN THE TRIDIAGONAL MATRIX 01610000 C DIPFLT----- I*ETA 01620000 C ******************************************************************** 01630000 C 01640000 WRITE(IPR,3) Q1INA,Q2INA,Q1INB,Q2INB 01640103 3 FORMAT(' *** 80 DEGREE OPTION ENTERED ***',/, 01640203 + ' Q1INA,Q2INA,Q1INB,Q2INB =',4F7.4) 01640303 C 01641003 DXQ=.25 * DX 01650000 N1=NX-1 01660000 C 01661000 Q1A=Q1INA-Q2INA 01670000 Q2A=2.*Q2INA 01680000 Q3A=2.*Q1A 01690000 RR11A = Q2A/DZ 01700000 RR31A = 2.*DX*DX/DZ 01710000 C 01710100 Q1B=Q1INB-Q2INB 01711000 Q2B=2.*Q2INB 01712000 Q3B=2.*Q1B 01713000 RR11B = Q2B/DZ 01714000 RR31B = 2.*DX*DX/DZ 01715000 C 01720000 SCF = 2.0/FLOAT(NTT) 01730000 C 01740000 AAI = 0 01750000 AAR = AAI + NX 01760000 AI = AAR + NX 01770000 AR = AI + NX 01780000 BBI = AR + NX 01790000 BBR = BBI + NX 01800000 BI = BBR + NX 01810000 BR = BI + NX 01820000 DENDEN = BR + NX 01830000 DENI = DENDEN+ NX 01840000 DENR = DENI + NX 01850000 DI = DENR + NX 01860000 DR = DI + NX 01870000 EI = DR + NX 01880000 ER = EI + NX 01890000 FI = ER + NX 01900000 FR = FI + NX 01910000 SHIFTI = FR + NX 01920000 SHIFTR = SHIFTI+ NX 01930000 TI = SHIFTR+ NX 01940000 TR = TI + NX 01950000 C 01960000 IF(IFLAG.EQ.IM) GO TO 10 01970000 RR12A=Q2A/DZ2 01980000 RR32A = 2.* DX * DX/DZ2 01990000 RR12B=Q2B/DZ2 01991000 RR32B = 2.* DX * DX/DZ2 01992000 IF(IPASS.EQ.2) DZ =DZ3 02000000 C 02010000 C ****************************************************************** 02020000 C THE DC COMPONENT IS NOT DOWNWARD CONTINUED 02030000 C ****************************************************************** 02040000 C 02050000 10 ISTART=1 02060000 IBLKP =0 02070000 IF(NOMLO.GT.1) GO TO 20 02080000 ISTART=2 02090000 NOMLO=2 02100000 C 02110000 C ******************************************************************* 02120000 C TAKE ONE FREQUENCY DOWN NZE-NZS+1 Z-STEPS 02130000 C ******************************************************************* 02140000 C 02150000 20 DO 280 IOM = ISTART,NOM 02160000 IBLK = (IOM+NFCPB-1)/NFCPB 02170000 IDIOM = MOD(IOM-1,NFCPB) + 1 02180000 C 02190000 IF(IBLK.EQ.IBLKP) GO TO 50 02200000 IBLKP = IBLK 02210000 DAWRK2 = IBLK 02220000 I = 1 02230000 C 02240000 DO 30 J=1,NTB 02250000 CALL FORDSD(KPWRK2,DAWRK2,BUF1(I)) 02260000 DAWRK2 = DAWRK2-1+NFB 02270000 CALL FORDSD(KPWRK2,DAWRK2,BUF2(I)) 02280000 I = I+NFCPB*NTPB 02290000 30 DAWRK2 = DAWRK2-1+NFB 02300000 C 02310000 CALL MTRNSP(BUF1,WAVR,NFCPB,NX) 02320000 CALL MTRNSP(BUF2,WAVI,NFCPB,NX) 02330000 C 02340000 IST = IOM 02350000 IF(IST.EQ.ISTART) IST=1 02360000 N = MIN0(NOM,IOM+NFCPB-1) 02370000 C 02380000 DO 40 J=IST,N 02390000 IF(NOMCHK(J).NE.0) GO TO 40 02400000 CALL ARSET(WAVR(1,J-IST+1),NX,DW) 02410000 CALL ARSET(WAVI(1,J-IST+1),NX,DW) 02420000 40 CONTINUE 02430000 C 02440000 C ******************************************************************* 02450000 C DEFINE CONSTANTS WHICH ARE FUNCTIONS OF OMEGA ONLY 02460000 C ******************************************************************* 02470000 C 02480000 50 IF(NOMCHK(IOM).NE.1) GO TO 250 02490000 OM =(NOMLO+IOM-2)*DOM 02500000 DZ=DZ3 02510000 RR1A = RR11A 02520000 RR3A = RR31A 02530000 RR1B = RR11B 02531000 RR3B = RR31B 02532000 OMDZ =OM*DZ 02540000 VIS = 1./CMPLX(OM,ETA) 02550000 RVIS = REAL(VIS) 02560000 CVIS = AIMAG(VIS) 02570000 C 02580000 C ******************************************************************* 02590000 C DOWNWARD CONTINUE THIS FREQUENCY COMPONENT 02600000 C ******************************************************************* 02610000 C 02620000 DO 240 IZ = NZS,NZE 02630000 C 02640000 C ******************************************************************* 02650000 C 02660000 C THE WAVEFIELD T(IX) IS DOWNWARD CONTINUED IN A TWO STEP PROCESS: 02670000 C 02680000 C (1) APPLY PHASE SHIFT : THIS IS THE THIN LENSE TERM OR THE 02690000 C SNELL'S LAW APPROXIMATION AND THE LINE OF CODE IS 02700000 C T(IX)= SHIFT * WAVE(IX,IOM) 02710000 C 02720000 C (2) DIFFRACT STEP 02730000 C 02740000 C THE DIFFRACT STEP IS: 02750000 C 02760000 C T(IX-1)*A(IX-1) + T(IX)*B(IX) + T(IX+1)*A(IX+1) 02770000 C = T+(IX-1)*AA(IX-1) + T+(IX)*BB(IX) + T+(IX+1)*AA(IX+1) 02780000 C WHERE 02790000 C T(IX) IS THE WAVEFIELD AT THE IZ LEVEL 02800000 C T+(IX) IS THE WAVEFIELD AT THE IZ+1 LEVEL 02810000 C 02820000 C --------------------------------------- 02830000 C | | | | 02840000 C | AA(IX-1) | BB(IX) | AA(IX+1) | IZ LEVEL 02850000 C | | | | 02860000 C |___________|____________|_____________| 02870000 C | | | | 02880000 C | A(IX-1) | B(IX) | A(IX+1) | IZ+1 LEVEL 02890000 C | | | | 02900000 C |___________|____________|_____________| 02910000 C 02920000 C IN MATRIX FORM THIS IS: 02930000 C ___ ___ __ __ __ __ 02940000 C | B(1) A(2) . | | T+(1) | | D(1) | 02950000 C | A(1) B(2) . | | T+(2) | | D(2) | 02960000 C | A(2) . . | | . | | . | 02970000 C | . . | | . | = | . | 02980000 C | . . | | . | | . | 02990000 C | . A(NX) | | . | | . | 03000000 C | A(NX-1) B(NX) | | T+(NX)| | D(NX)| 03010000 C |__ ___| |__ __| |__ __| 03020000 C 03030000 C WHERE THE D MATRIX IS DEFINED BY : 03040000 C 03050000 C ____ _____ ___ __ ___ ___ 03060000 C | BB(1) AA(2) . | | T(1) | | D(1) | 03070000 C | AA(1) BB(2) . | | T(2) | | D(2) | 03080000 C | AA(2) . . | | . | | . | 03090000 C | . . | | . | = | . | 03100000 C | . . | | . | | . | 03110000 C | . AA(NX) | | . | | . | 03120000 C | AA(NX-1)BB(NX) | | T(NX)| | D(NX)| 03130000 C |____ ____| |__ __| |__ __| 03140000 C 03150000 C DZ=DZ3 03160000 C RR1 = RR11 03170000 C RR3 = RR31 03180000 IF(IFLAG.EQ.IM .OR. IZ.NE.NZE) GO TO 60 03190000 DZ=DZ2 03200000 RR1A=RR12A 03210000 RR3A=RR32A 03220000 RR1B=RR12B 03221000 RR3B=RR32B 03222000 OMDZ=OM*DZ 03230000 C IF(DZ.LT.0.5) GO TO 240 03240000 C 03250000 C==================================================================== 03251000 C PASS A 03252000 C==================================================================== 03253000 CJCS 03254000 60 CONTINUE 03260000 IF(Q1INA.EQ.0.) GO TO 6666 03260100 CJCS60 IF((IVCHK(IZ).EQ.1).AND.(IZ.NE.NZS)) GO TO 160 03261000 C 03270000 C********************************************************************* 03280000 C FIND THE TRI-D ELEMENTS FOR THIS LAYER 03290000 C********************************************************************* 03300000 C 03310000 IBGN = 1 03320000 C 03330000 IF(IZ.EQ.NZS) GO TO 70 03340000 DAWRK = 3*(IZ-1)-2 03350000 CMEMORY CALL FORDSD(KPWRKD,DAWRK,VEL(1,1)) 03360000 CALL ARMVE( V(1,DAWRK), VEL(1,1), NX) 03370000 70 DAWRK = 3* IZ -2 03380000 CMEMORY CALL FORDSD(KPWRKD,DAWRK,VEL(1,2)) 03390000 CALL ARMVE( V(1,DAWRK), VEL(1,2), NX) 03400000 CMEMORY CALL FORDSD(KPWRKD,DAWRK, ISTO) 03410000 CALL ARMVE( V(1,DAWRK+1), ISTO, NX) 03420000 C 03430000 DO 90 KK = 1 ,NX 03440000 KEND =ISTO(KK) 03450000 VCT = VEL(KEND, 2) 03460000 RRSET =RR1A*VCT 03470000 CC3 = RR3A*OM/VCT 03480000 COMP =(RVIS * RRSET) + (CC3 *BETA) 03490000 RREAL = -(RRSET * CVIS) 03500000 ZAAR = Q1A-RREAL 03510000 ZAAI = - COMP 03520000 ZAR = ZAAR-Q3A 03530000 CJCS ZBBR = 2.*RREAL-1. 03541000 ZBBR = 2.*RREAL-Q3A 03542000 ZBBI = 2. * COMP -CC3 03550000 CJCS ZBR = ZBBR + 2. 03560000 ZBR = ZBBR + 2. * Q3A 03561000 ZSFTR= COS(OMDZ/VCT) 03570000 ZSFTI= SIN(OMDZ/VCT) 03580000 C 03590000 CXIR$ IVDEP 03600000 DO 80 IK = IBGN,KEND 03610000 S(AAR+IK) = ZAAR 03620000 S(AAI+IK) = ZAAI 03630000 S(AR+IK) = ZAR 03640000 S(AI+IK) = ZAAI 03650000 S(BBR+IK) = ZBBR 03660000 S(BBI+IK) = ZBBI 03670000 S(BR+IK) = ZBR 03680000 S(BI+IK) = ZBBI 03690000 S(SHIFTR+IK) = ZSFTR 03700000 S(SHIFTI+IK) = ZSFTI 03710000 80 CONTINUE 03720000 C 03730000 IF(KEND.EQ.NX) GO TO 100 03740000 90 IBGN =KEND+1 03750000 C 03760000 C ******************************************************************** 03770000 C ABSORBING SIDE CONDITIONS : FOR EXPLANATION SEE JEFF 03780000 C THORSON'S ARTICLE IN SEP REPORT #20 03790000 C ******************************************************************** 03800000 C 03810000 C LEFT BOUNDARIES 03820000 CJCS 03821000 100 CONTINUE 03830000 C 100 IF(IZ.EQ.NZS) GO TO 110 03831000 C IF(VEL(1, 2).EQ.VEL(1,1)) GO TO 120 03840000 110 BAB1= OM/VEL(1,2) * DXQ 03850000 RA1D = 1. + BAB1*BAB1 03860000 RA1R = (1 - BAB1*BAB1)/RA1D 03870000 RA1I = (BAB1+BAB1)/RA1D 03880000 BCR = S(AR+1)*RA1R - S(AI+1)*RA1I 03890000 BCI = S(AR+1)*RA1I + S(AI+1)*RA1R 03900000 BBCR = S(AAR+1)*RA1R - S(AAI+1)*RA1I 03910000 BBCI = S(AAR+1)*RA1I + S(AAI+1)*RA1R 03920000 BR1 = S(BR+1) + BCR 03930000 BI1 = S(BI+1) + BCI 03940000 BBR1= S(BBR+1) + BBCR 03950000 BBI1= S(BBI+1) + BBCI 03960000 C 03970000 C RIGHT BOUNDARIES 03980000 120 CONTINUE 03990000 C 120 IF(IZ.EQ.NZS) GO TO 130 03991000 C IF(VEL(NX, 2).EQ.VEL(NX,1)) GO TO 140 04000000 130 BAB1= OM/VEL(NX,2) * DXQ 04010000 RA1D = 1. + BAB1*BAB1 04020000 RA1R = (1 - BAB1*BAB1)/RA1D 04030000 RA1I = (BAB1+BAB1)/RA1D 04040000 BCR = S(AR+NX)*RA1R - S(AI+NX)*RA1I 04050000 BCI = S(AR+NX)*RA1I + S(AI+NX)*RA1R 04060000 BBCR = S(AAR+NX)*RA1R - S(AAI+NX)*RA1I 04070000 BBCI = S(AAR+NX)*RA1I + S(AAI+NX)*RA1R 04080000 BRNX = S(BR+NX) + BCR 04090000 BINX = S(BI+NX) + BCI 04100000 BBRNX= S(BBR+NX) + BBCR 04110000 BBINX= S(BBI+NX) + BBCI 04120000 C 04130000 C ****************************************************************** 04140000 C THE FOLLOWING NINE LINES ARE PART OF THE CRANK NICHOLON MATRIX 04150000 C SOLUTION,BUT DO NOT CHANGE IF THE VELOCITIES ON ONE Z-LEVEL 04160000 C IS CONSTANT 04170000 C ****************************************************************** 04180000 C 04190000 140 DB = BR1 * BR1 + BI1 * BI1 04200000 S(ER+1) =-(BR1 * S(AR+2) + BI1* S(AI+2))/DB 04210000 S(EI+1) = (BI1 * S(AR+2) - BR1* S(AI+2))/DB 04220000 C 04230000 CXIR$ IVDEP 04240000 DO 150 I =2 ,N1 04250000 S(DENR+I)=S(BR+I)+S(AR+I-1)*S(ER+I-1)-S(AI+I-1)*S(EI 04260000 * +I-1) 04270000 S(DENI+I) = S(BI+I) + S(AI+I-1)*S(ER+I-1)+S(AR+I-1)* 04280000 * S(EI+I-1) 04290000 S(DENDEN+I) = S(DENR+I) * S(DENR+I) + S(DENI+I)*S(DE 04300000 * NI+I) 04310000 S(ER+I) =-(S(AR+I+1)*S(DENR+I)+S(AI+I+1)*S(DENI+I))/ 04320000 * S(DENDEN+I) 04330000 S(EI+I) = (S(AR+I+1)*S(DENI+I)-S(AI+I+1)*S(DENR+I))/ 04340000 * S(DENDEN+I) 04350000 150 CONTINUE 04360000 C 04370000 C ***************************************************************** 04380000 C COMPUTE T AND D ARRAYS 04390000 C ***************************************************************** 04400000 C 04410000 CXIR$ IVDEP 04420000 160 DO 170 IX=1,NX 04430000 S(TR+IX)=S(SHIFTR+IX)*(WAVR(IX,IDIOM))-S(SHIFTI+IX)* 04440000 * (WAVI(IX,IDIOM)) 04450000 S(TI+IX)=S(SHIFTI+IX)*(WAVR(IX,IDIOM))+S(SHIFTR+IX)* 04460000 * (WAVI(IX,IDIOM)) 04470000 170 CONTINUE 04480000 C 04490000 S(DR+1)=BBR1*S(TR+1)-BBI1*S(TI+1)+S(AAR+2)*S(TR+2)-S(A 04500000 * AI+2)*S(TI+2) 04510000 S(DI+1)=BBR1*S(TI+1)+BBI1*S(TR+1)+S(AAR+2)*S(TI+2)+S(A 04520000 * AI+2)*S(TR+2) 04530000 S(DR+NX)=BBRNX*S(TR+NX)-BBINX*S(TI+NX)+S(AAR+N1)*S(TR+ 04540000 * N1)-S(AAI+N1)*S(TI+N1) 04550000 S(DI+NX)=BBRNX*S(TI+NX)+BBINX*S(TR+NX)+S(AAR+N1)*S(TI+ 04560000 * N1)+S(AAI+N1)*S(TR+N1) 04570000 C 04580000 CXIR$ IVDEP 04590000 DO 180 IX=2,N1 04600000 S(DR+IX)=S(BBR+IX)*S(TR+IX)-S(BBI+IX)*S(TI+IX)+S(AAR 04610000 * +IX-1)*S(TR+IX-1)-S(AAI+IX-1)*S(TI+IX-1)+ S(AAR+IX+1 04620000 * )*S(TR+IX+1)-S(AAI+IX+1)*S(TI+IX+1) 04630000 S(DI+IX)=S(BBR+IX)*S(TI+IX)+S(BBI+IX)*S(TR+IX)+S(AAR 04640000 * +IX-1)*S(TI+IX-1)+S(AAI+IX-1)*S(TR+IX-1)+ S(AAR+IX+1 04650000 * )*S(TI+IX+1)+S(AAI+IX+1)*S(TR+IX+1) 04660000 180 CONTINUE 04670000 C 04680000 C ****************************************************************** 04690000 C SOLVE TRI-DIAGONAL SYSTEM OF EQUATIONS BY THE CRANK 04700000 C NICHOLSON METHOD 04710000 C ****************************************************************** 04720000 C 04730000 S(FR+1) =(BR1*S(DR+1)+BI1*S(DI+1))/DB 04740000 S(FI+1) =(BR1*S(DI+1)-BI1*S(DR+1))/DB 04750000 C 04760000 CDIR$ NEXTSCALAR 04770000 DO 190 ITRI = 2,N1 04780000 FTR=S(DR+ITRI)-S(AR+ITRI-1)*S(FR+ITRI-1)+S(AI+ITRI-1 04790000 * )*S(FI+ITRI-1) 04800000 FTI=S(DI+ITRI)-S(AR+ITRI-1)*S(FI+ITRI-1)-S(AI+ITRI-1 04810000 * )*S(FR+ITRI-1) 04820000 S(FR+ITRI)=(FTR*S(DENR+ITRI)+FTI*S(DENI+ITRI))/S(DEN 04830000 * DEN+ITRI) 04840000 S(FI+ITRI)=(FTI*S(DENR+ITRI)-FTR*S(DENI+ITRI))/S(DEN 04850000 * DEN+ITRI) 04860000 190 CONTINUE 04870000 C 04880000 TTNR=S(DR+NX)-S(AR+N1)*S(FR+N1)+S(AI+N1)*S(FI+N1) 04890000 TTNI=S(DI+NX)-S(AR+N1)*S(FI+N1)-S(AI+N1)*S(FR+N1) 04900000 TTDR=BRNX+S(AR+N1)*S(ER+N1)-S(AI+N1)*S(EI+N1) 04910000 TTDI=BINX+S(AR+N1)*S(EI+N1)+S(AI+N1)*S(ER+N1) 04920000 TTD2 = TTDR*TTDR + TTDI*TTDI 04930000 S(TR+NX)=(TTNR*TTDR+TTNI*TTDI)/TTD2 04940000 S(TI+NX)=(TTNI*TTDR-TTNR*TTDI)/TTD2 04950000 C 04960000 CDIR$ NEXTSCALAR 04970000 DO 200 JTRI =1,N1 04980000 IRI=NX-JTRI 04990000 S(TR+IRI)=S(ER+IRI)*S(TR+IRI+1)-S(EI+IRI)*S(TI+IRI+1 05000000 * ) +S( FR+IRI) 05010000 S(TI+IRI)=S(ER+IRI)*S(TI+IRI+1)+S(EI+IRI)*S(TR+IRI+1 05020000 * ) +S( FI+IRI) 05030000 200 CONTINUE 05040000 C 05050000 C==================================================================== 05051000 C PASS B 05052000 C==================================================================== 05053000 C 05053100 6666 IF(Q1INB.EQ.0.) GO TO 6667 05053200 C 05058000 C********************************************************************* 05059000 C FIND THE TRI-D ELEMENTS FOR THIS LAYER 05059100 C********************************************************************* 05059200 C 05059300 IBGN = 1 05059400 C 05059500 C IF(IZ.EQ.NZS) GO TO 70 05059600 C DAWRK = 3*(IZ-1)-2 05059700 CMEMORY CALL FORDSD(KPWRKD,DAWRK,VEL(1,1)) 05059800 C CALL ARMVE( V(1,DAWRK), VEL(1,1), NX) 05059900 C 70 DAWRK = 3* IZ -2 05060000 CMEMORY CALL FORDSD(KPWRKD,DAWRK,VEL(1,2)) 05060100 C CALL ARMVE( V(1,DAWRK), VEL(1,2), NX) 05060200 CMEMORY CALL FORDSD(KPWRKD,DAWRK, ISTO) 05060300 C CALL ARMVE( V(1,DAWRK+1), ISTO, NX) 05060400 C 05060500 DO 1090 KK = 1 ,NX 05060600 KEND =ISTO(KK) 05060700 VCT = VEL(KEND, 2) 05060800 RRSET =RR1B*VCT 05060900 CC3 = RR3B*OM/VCT 05061000 COMP =(RVIS * RRSET) + (CC3 *BETA) 05061100 RREAL = -(RRSET * CVIS) 05061200 ZAAR = Q1B-RREAL 05061300 ZAAI = - COMP 05061400 ZAR = ZAAR-Q3B 05061500 CJCS ZBBR = 2.*RREAL-1. 05061600 ZBBR = 2.*RREAL-Q3B 05061700 ZBBI = 2. * COMP -CC3 05061800 CJCS ZBR = ZBBR + 2. 05061900 ZBR = ZBBR + 2. * Q3B 05062000 ZSFTR= COS(OMDZ/VCT) 05062100 ZSFTI= SIN(OMDZ/VCT) 05062200 C 05062300 CXIR$ IVDEP 05062400 DO 1080 IK = IBGN,KEND 05062500 S(AAR+IK) = ZAAR 05062600 S(AAI+IK) = ZAAI 05062700 S(AR+IK) = ZAR 05062800 S(AI+IK) = ZAAI 05062900 S(BBR+IK) = ZBBR 05063000 S(BBI+IK) = ZBBI 05063100 S(BR+IK) = ZBR 05063200 S(BI+IK) = ZBBI 05063300 S(SHIFTR+IK) = ZSFTR 05063400 S(SHIFTI+IK) = ZSFTI 05063500 1080 CONTINUE 05063600 C 05063700 IF(KEND.EQ.NX) GO TO 1100 05063800 1090 IBGN =KEND+1 05063900 C 05064000 C ******************************************************************** 05064100 C ABSORBING SIDE CONDITIONS : FOR EXPLANATION SEE JEFF 05064200 C THORSON'S ARTICLE IN SEP REPORT #20 05064300 C ******************************************************************** 05064400 C 05064500 C LEFT BOUNDARIES 05064600 CJCS 05064700 1100 CONTINUE 05064801 C 100 IF(IZ.EQ.NZS) GO TO 110 05064900 C IF(VEL(1, 2).EQ.VEL(1,1)) GO TO 120 05065000 1110 BAB1= OM/VEL(1,2) * DXQ 05065100 RA1D = 1. + BAB1*BAB1 05065200 RA1R = (1 - BAB1*BAB1)/RA1D 05065300 RA1I = (BAB1+BAB1)/RA1D 05065400 BCR = S(AR+1)*RA1R - S(AI+1)*RA1I 05065500 BCI = S(AR+1)*RA1I + S(AI+1)*RA1R 05065600 BBCR = S(AAR+1)*RA1R - S(AAI+1)*RA1I 05065700 BBCI = S(AAR+1)*RA1I + S(AAI+1)*RA1R 05065800 BR1 = S(BR+1) + BCR 05065900 BI1 = S(BI+1) + BCI 05066000 BBR1= S(BBR+1) + BBCR 05066100 BBI1= S(BBI+1) + BBCI 05066200 C 05066300 C RIGHT BOUNDARIES 05066400 1120 CONTINUE 05066500 C 120 IF(IZ.EQ.NZS) GO TO 130 05066600 C IF(VEL(NX, 2).EQ.VEL(NX,1)) GO TO 140 05066700 1130 BAB1= OM/VEL(NX,2) * DXQ 05066800 RA1D = 1. + BAB1*BAB1 05066900 RA1R = (1 - BAB1*BAB1)/RA1D 05067000 RA1I = (BAB1+BAB1)/RA1D 05067100 BCR = S(AR+NX)*RA1R - S(AI+NX)*RA1I 05067200 BCI = S(AR+NX)*RA1I + S(AI+NX)*RA1R 05067300 BBCR = S(AAR+NX)*RA1R - S(AAI+NX)*RA1I 05067400 BBCI = S(AAR+NX)*RA1I + S(AAI+NX)*RA1R 05067500 BRNX = S(BR+NX) + BCR 05067600 BINX = S(BI+NX) + BCI 05067700 BBRNX= S(BBR+NX) + BBCR 05067800 BBINX= S(BBI+NX) + BBCI 05067900 C 05068000 C ****************************************************************** 05068100 C THE FOLLOWING NINE LINES ARE PART OF THE CRANK NICHOLON MATRIX 05068200 C SOLUTION,BUT DO NOT CHANGE IF THE VELOCITIES ON ONE Z-LEVEL 05068300 C IS CONSTANT 05068400 C ****************************************************************** 05068500 C 05068600 1140 DB = BR1 * BR1 + BI1 * BI1 05068700 S(ER+1) =-(BR1 * S(AR+2) + BI1* S(AI+2))/DB 05068800 S(EI+1) = (BI1 * S(AR+2) - BR1* S(AI+2))/DB 05068900 C 05069000 CXIR$ IVDEP 05069100 DO 1150 I =2 ,N1 05069200 S(DENR+I)=S(BR+I)+S(AR+I-1)*S(ER+I-1)-S(AI+I-1)*S(EI 05069300 * +I-1) 05069400 S(DENI+I) = S(BI+I) + S(AI+I-1)*S(ER+I-1)+S(AR+I-1)* 05069500 * S(EI+I-1) 05069600 S(DENDEN+I) = S(DENR+I) * S(DENR+I) + S(DENI+I)*S(DE 05069700 * NI+I) 05069800 S(ER+I) =-(S(AR+I+1)*S(DENR+I)+S(AI+I+1)*S(DENI+I))/ 05069900 * S(DENDEN+I) 05070000 S(EI+I) = (S(AR+I+1)*S(DENI+I)-S(AI+I+1)*S(DENR+I))/ 05070100 * S(DENDEN+I) 05070200 1150 CONTINUE 05070300 C 05070400 C ***************************************************************** 05070500 C COMPUTE T AND D ARRAYS 05070600 C ***************************************************************** 05070700 C 05070800 C 05071600 S(DR+1)=BBR1*S(TR+1)-BBI1*S(TI+1)+S(AAR+2)*S(TR+2)-S(A 05071700 * AI+2)*S(TI+2) 05071800 S(DI+1)=BBR1*S(TI+1)+BBI1*S(TR+1)+S(AAR+2)*S(TI+2)+S(A 05071900 * AI+2)*S(TR+2) 05072000 S(DR+NX)=BBRNX*S(TR+NX)-BBINX*S(TI+NX)+S(AAR+N1)*S(TR+ 05072100 * N1)-S(AAI+N1)*S(TI+N1) 05072200 S(DI+NX)=BBRNX*S(TI+NX)+BBINX*S(TR+NX)+S(AAR+N1)*S(TI+ 05072300 * N1)+S(AAI+N1)*S(TR+N1) 05072400 C 05072500 CXIR$ IVDEP 05072600 DO 1180 IX=2,N1 05072700 S(DR+IX)=S(BBR+IX)*S(TR+IX)-S(BBI+IX)*S(TI+IX)+S(AAR 05072800 * +IX-1)*S(TR+IX-1)-S(AAI+IX-1)*S(TI+IX-1)+ S(AAR+IX+1 05072900 * )*S(TR+IX+1)-S(AAI+IX+1)*S(TI+IX+1) 05073000 S(DI+IX)=S(BBR+IX)*S(TI+IX)+S(BBI+IX)*S(TR+IX)+S(AAR 05073100 * +IX-1)*S(TI+IX-1)+S(AAI+IX-1)*S(TR+IX-1)+ S(AAR+IX+1 05073200 * )*S(TI+IX+1)+S(AAI+IX+1)*S(TR+IX+1) 05073300 1180 CONTINUE 05073400 C 05073500 C ****************************************************************** 05073600 C SOLVE TRI-DIAGONAL SYSTEM OF EQUATIONS BY THE CRANK 05073700 C NICHOLSON METHOD 05073800 C ****************************************************************** 05073900 C 05074000 S(FR+1) =(BR1*S(DR+1)+BI1*S(DI+1))/DB 05074100 S(FI+1) =(BR1*S(DI+1)-BI1*S(DR+1))/DB 05074200 C 05074300 CDIR$ NEXTSCALAR 05074400 DO 1190 ITRI = 2,N1 05074500 FTR=S(DR+ITRI)-S(AR+ITRI-1)*S(FR+ITRI-1)+S(AI+ITRI-1 05074600 * )*S(FI+ITRI-1) 05074700 FTI=S(DI+ITRI)-S(AR+ITRI-1)*S(FI+ITRI-1)-S(AI+ITRI-1 05074800 * )*S(FR+ITRI-1) 05074900 S(FR+ITRI)=(FTR*S(DENR+ITRI)+FTI*S(DENI+ITRI))/S(DEN 05075000 * DEN+ITRI) 05075100 S(FI+ITRI)=(FTI*S(DENR+ITRI)-FTR*S(DENI+ITRI))/S(DEN 05075200 * DEN+ITRI) 05075300 1190 CONTINUE 05075400 C 05075500 TTNR=S(DR+NX)-S(AR+N1)*S(FR+N1)+S(AI+N1)*S(FI+N1) 05075600 TTNI=S(DI+NX)-S(AR+N1)*S(FI+N1)-S(AI+N1)*S(FR+N1) 05075700 TTDR=BRNX+S(AR+N1)*S(ER+N1)-S(AI+N1)*S(EI+N1) 05075800 TTDI=BINX+S(AR+N1)*S(EI+N1)+S(AI+N1)*S(ER+N1) 05075900 TTD2 = TTDR*TTDR + TTDI*TTDI 05076000 S(TR+NX)=(TTNR*TTDR+TTNI*TTDI)/TTD2 05076100 S(TI+NX)=(TTNI*TTDR-TTNR*TTDI)/TTD2 05076200 C 05076300 CDIR$ NEXTSCALAR 05076400 DO 1200 JTRI =1,N1 05076500 IRI=NX-JTRI 05076600 S(TR+IRI)=S(ER+IRI)*S(TR+IRI+1)-S(EI+IRI)*S(TI+IRI+1 05076700 * ) +S( FR+IRI) 05076800 S(TI+IRI)=S(ER+IRI)*S(TI+IRI+1)+S(EI+IRI)*S(TR+IRI+1 05076900 * ) +S( FI+IRI) 05077000 1200 CONTINUE 05077100 C 05077200 C ****************************************************************** 05077300 C IMAGE THE DATA BY SUMMING OVER THE REAL PART OF THE 05078000 C FREQUENCIES TO GET T(0). 05080000 C NOTE: NO IMAGING IS DONE WHEN DATUMING 05090000 C ****************************************************************** 05100000 C 05110000 6667 CONTINUE 05111000 210 DO 220 IX = 1,NX 05120000 WAVR(IX,IDIOM) = S(TR+IX) 05130000 WAVI(IX,IDIOM) = S(TI+IX) 05140000 220 CONTINUE 05150000 C 05160000 IF(IFLAG.EQ.ID) GO TO 240 05170000 C 05180000 DAWRK = 3*IZ 05190000 CMEMORY CALL FORDSD(KPWRKD,DAWRK,REF) 05200000 CALL ARMVE( V(1,DAWRK), REF, NX) 05210000 C 05220000 DO 230 IX=1,NX 05230000 230 REF(IX) = REF(IX) + S(TR+IX) 05240000 C 05250000 IF(IOM.EQ.NOMLST) CALL ARMPFC(REF,REF,SCF,NX) 05260000 C 05270000 CMEMORY DAWRK = DAWRK-1 05280000 CMEMORY CALL FOWDSD(KPWRKD,DAWRK,REF) 05290000 CALL ARMVE( REF, V(1,DAWRK), NX) 05300000 C 05310000 240 CONTINUE 05320000 C 05330000 250 IF(IOM.EQ.NOM) GO TO 260 05340000 IF(IDIOM.NE.NFCPB) GO TO 280 05350000 C 05360000 260 CALL MTRNSP(WAVR,BUF1,NX,NFCPB) 05370000 CALL MTRNSP(WAVI,BUF2,NX,NFCPB) 05380000 C 05390000 DAWRK2 = IBLK 05400000 I=1 05410000 C 05420000 DO 270 J=1,NTB 05430000 CALL FOWDSD(KPWRK2,DAWRK2,BUF1(I)) 05440000 DAWRK2 = DAWRK2-1+NFB 05450000 CALL FOWDSD(KPWRK2,DAWRK2,BUF2(I)) 05460000 I=I+NFCPB*NTPB 05470000 270 DAWRK2=DAWRK2-1+NFB 05480000 C 05490000 280 CONTINUE 05500000 C 05510000 IF(IPASS.EQ.2) GO TO 320 05520000 IF(IFLAG.NE.ID) GO TO 320 05530000 C 05540000 C ******************************************************************** 05550000 C INSERT REPLACEMENT VELOCITY IN FILE FOR DATUMING 05560000 C ******************************************************************** 05570000 C 05580000 290 IF(IUP.EQ.1) GO TO 320 05590000 C 05600000 IPASS=2 05610000 NZE1=NZE+1 05620000 COADY 05630000 VTMP = -0.5*IVDATM 05640000 CALL ARSET(VEL(1,1),NX,VTMP) 05650000 C 05660000 DO 300 IZ=1,NZ 05670000 DAWRK = IZ*3 - 2 05680000 CMEMORY CALL FOWDSD(KPWRKD,DAWRK,VEL(1,1)) 05690000 CALL ARMVE (VEL(1,1), V(1,DAWRK), NX) 05700000 300 CONTINUE 05710000 C 05720000 DO 310 IZ=1,NZE1 05730000 310 IVCHK(IZ)=1 05740000 C 05750000 ETA=-ETA 05760000 C 05770000 320 RETURN 05780000 END 05790000