CTITLEMNOQD3 -- REMOVE CONSTANT AND QUAD. LEAST SQUARES COMPONENTS 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D.D. THOMPSON 00000200 CA DESIGNER D.D. THOMPSON 00000300 CA LANGUAGE FORTRAN H 00000400 CA SYSTEM IBM OR CRAY 00000500 CA WRITTEN SEPTEMBER, 1979 00000600 C REVISED MO-DA-YR BY PROGRAMMER FOR REASON 00000700 C REVISED 07-07-86 JMP. DUAL IBM/CRAY VERSION. 00000800 C REVISED 05-05-88 LWC. DECLARE SYSTEM VARIABLES INTEGER. 00000800 C REVISED 11-13-89 RDK. REMOVE EXTERNAL FOR S1ATP. 00000800 C 00000900 CA 00001000 CA 00001100 CA CALL MNOQD3 (Q, ICDP, KCDP, LINE, LD, NLINE, D, XCDP, XRNMO, 00001200 CA ZA, ZB, ZC, NUMC, NUML) 00001300 CA 00001400 CA IN/OUT ARGUMENT TYPE DESCRIPTION 00001500 CA 00001600 CA IN/OUT Q R8 LAG VALUES FOR EACH TRACE. LENGTH LINE(NLINE)00001700 CA (REJECTED TRACES ARE INDICATED BY Q > 1.0E6) 00001800 CA IN ICDP I4 ARRAY OF INDEX POSITIONS IN IND CORRESPONDING00001900 CA TO THE LAST TRACE IN EACH CDP. LENGTH < OR = 00002000 CA NLINE*LD 00002100 CA IN KCDP I4 ARRAY OF STARTING CDP NUMBERS FOR EACH LINE. 00002200 CA LENGTH = NLINE 00002300 CA IN LINE I4 POINTER ARRAY INDICATING LAST ELEMENT IN Q 00002400 CA FOR EACH LINE. LENGTH = NLINE 00002500 CA IN LD I4 MAXIMUM CDP # WHICH OCCURS ON ANY LINE. 00002600 CA IN NLINE I4 NUMBER OF LINES. 00002700 CA IN D R4 OFFSET SQUARED FOR EACH CDP. LENGTH = 00002800 CA LINE(NLINE) 00002900 CA OUT XCDP R8 OUTPUT CDP AVERAGE FOR EACH CDP. ARRANGED 00003000 CA BY LINES WITH LD ENTRIES/LINE. NOTE SOME 00003100 CA OF THESE CDP ARE NOT REPRESENTED BY TRACES. 00003200 CA LENGTH = LD*NLINE 00003300 CA OUT XRNMO R8 OUTPUT RNMO FOR EACH CDP ARRANGED AS XCDP. 00003400 CA LENGTH = LD*NLINE 00003500 CA IN ZA R8 SCRATCH ARRAY OF LENGTH = NLINE*LD 00003600 CA IN ZB R8 SCRATCH ARRAY OF LENGTH = NLINE*LD 00003700 CA IN ZC R8 SCRATCH ARRAY OF LENGTH = NLINE*LD 00003800 CA IN NUMC I4 # CDP ON EACH SIDE OF CENTER IN RNMO AVERAGE 00003900 CA IN NUML I4 #LINES ON EACH SIDE OF CENTER IN RNMO AVERAGE00004000 CA 00004100 CA THIS ROUTINE REMOVES THE CONSTANT AND QUADRATIC COMPONENTS 00004200 CA (WITH RESPECT TO OFFSET) FROM THE INPUT LAG VALUES. 00004300 CA THE QUADRATIC COMPONENT IS AVERAGED OVER A SPECIFIED 00004400 CA RECTANGLE OF ADJACENT CDP FOR 3-D DATA. THE COEFFICIENTS OF 00004500 CA THE REMOVED COMPONENTS ARE ALSO RETURNED. 00004600 CA 00004700 CA NOTE: AT ENDS AND CORNERS OF 3-D GRID THE RNMO AVERAGES 00004800 CA ARE COMPUTED OVER WHATEVER TRACES ARE AVAILABLE 00004900 CA ON THE AVERAGING RECTANGLE. THUS THE AVERAGES ARE 00005000 CA NOT SYMMETRIC ABOUT THE CDP'S BEING EVALUATED. 00005100 CA 00005200 C *** SEE THE FIRST PART OF APPENDIX A IN THE TECHNICAL NOTES 00005300 C FOR THE PROCESS 'TRAX'. COMMENTS HERE WILL REFERENCE THE 00005400 C STEP NUMBERS AND VARIABLE NAMES DESCRIBED THERE. 00005500 C 00005600 C 00005700 SUBROUTINE MNOQD3 (Q, ICDP, KCDP, LINE, LD, NLINE, D, XCDP, XRNMO,00005800 * ZA, ZB, ZC, NUMC, NUML) 00005900 C 00006000 C 00990000 C *********** WARNING ************************************* 01000000 C *** NO IMPLICIT INTEGER STATEMENT **************** 01010000 C ********************************************************* 01020000 C 01030000 C EXTERNAL S1ATP 00006100 C 00006200 C INTEGER ARRAYS IN PARAMETER LIST 00006300 C 00006400 INTEGER ICDP (1) 00006500 INTEGER KCDP (1) 00006600 INTEGER LINE (1) 00006700 INTEGER SYSTEM 01330000 INTEGER SYBYPW 01340000 INTEGER SYLOCF 01350000 INTEGER JAPNMS 01360000 C 00006800 C REAL ARRAYS IN PARAMETER LIST 00006900 C 00007000 REAL D (1) 00007100 DOUBLE PRECISION Q (1) 00007200 DOUBLE PRECISION XCDP (1) 00007300 DOUBLE PRECISION XRNMO (1) 00007400 DOUBLE PRECISION ZA (1) 00007500 DOUBLE PRECISION ZB (1) 00007600 DOUBLE PRECISION ZC (1) 00007700 C 00007800 C REAL VARIABLES--LOCAL 00007900 C 00008000 DOUBLE PRECISION B 00008100 DOUBLE PRECISION Q2 00008200 DOUBLE PRECISION SD 00008300 DOUBLE PRECISION SW 00008400 DOUBLE PRECISION SX 00008500 DOUBLE PRECISION SX2 00008600 DOUBLE PRECISION SY 00008700 DOUBLE PRECISION SYD 00008800 C 00008900 COMMON /SYSTEM/ SYSTEM, SYBYPW, SYLOCF, JAPNMS 00009000 C IF (1.EQ.2) CALL S1ATP C 00010000 IF (SYBYPW .EQ. 8) THEN 00020000 NDPW = 1 00030000 ELSE 00040000 NDPW = 2 00050000 ENDIF 00050100 C 00050200 C DETERMINE THE MAXIMUM VALUE OF DEPTH POINTS 00050300 C 00050400 LDD = LD * NLINE 00050500 LDD2 = NDPW * LDD 00050600 C 00050700 C SET NUMBER OF CDP (MUM) AND NUMBER OF LINES (LUM) FOR RNMO 00050800 C AVERAGING. THESE ARE NORMALLY SET TO PARAMETER VALUES NUMC 00050900 C AND NUML BUT CAN NOT BE GREATER THAN 1/4 NUMBER OF CDP OR 00051000 C LINES IN DATA. 00051100 C 00051200 MUM = NUMC 00051300 IF (MUM*4+3 .GT. LD) MUM = (LD-3) / 4 00051400 MUMP = MUM + 1 00051500 LDM = LD - MUM 00051600 LDMP = LDM + 1 00051700 LUM = NUML 00051800 IF (LUM*4+3 .GT. NLINE) LUM = (NLINE-3) / 4 00051900 LUMP = LUM + 1 00052000 NLNM = NLINE - LUM 00052100 NLNMP = NLNM + 1 00052200 LUMLD = LUM * LD 00052300 C 00052400 C CLEAR SCRATCH ARRAYS AND XCDP AND XRNMO ARRAYS. 00052500 C 00052600 CALL ARSET(ZA, LDD2, 0.) 00052700 CALL ARSET(ZB, LDD2, 0.) 00052800 CALL ARSET(ZC, LDD2, 0.) 00052900 CALL ARSET(XCDP, LDD2, 0.) 00053000 CALL ARSET(XRNMO, LDD2, 0.) 00053100 C 00053200 C PREPARE TO EVALUATE CDP SUMS DESCRIBED IN STEP 1. 00053300 C 00053400 II = 0 00053500 III = 0 00053600 KEY2 = 1 00053700 KEY = 0 00053800 L = 1 00053900 KK = 0 00054000 C 00054100 C INITIALIZE SUMS FOR NEW CDP. 00054200 C 00054300 10 SW = 0. 00054400 SYD = 0. 00054500 SY = 0. 00054600 SX = 0. 00054700 SX2 = 0. 00054800 KK = KK + 1 00054900 III = III + 1 00055000 IIIC = ICDP(III) 00055100 IF (KEY2 .EQ. 0) GO TO 20 00055200 C 00055300 C UPDATE POINTERS, ETC. FOR NEW LINES 00055400 C 00055500 II = II + 1 00055600 IILD = (II-1) * LD 00055700 NTT = LINE(II) 00055800 KEY2 = 0 00055900 KK = KCDP(II) 00056000 C 00056100 20 DO 30 I = L, IIIC 00056200 Q2 = Q(I) 00056300 C 00056400 C IF TRACE NOT REJECTED , UPDATE SUMS 00056500 C 00056600 IF(Q2 .GT. 1.0E6) GO TO 30 00056700 D2 = D(I) 00056800 SW = SW + 1.0 00056900 SYD = Q2 * D2 + SYD 00057000 SY = Q2 + SY 00057100 SX = D2 + SX 00057200 SX2 = D2 * D2 + SX2 00057300 C 00057400 30 CONTINUE 00057500 C 00057600 C END OF CDP. UPDATE POINTERS AND IF ALL TRACES ARE NOT REJECTED 00057700 C STORE KERNELS DESCRIBED IN STEP 1. 00057800 C 00057900 L = IIIC + 1 00058000 IF (IIIC .EQ. NTT) KEY2 = 1 00058100 IF (II .EQ. NLINE .AND. IIIC .EQ. NTT) KEY = 1 00058200 IF (SW .EQ. 0.) GO TO 40 00058300 INDX = KK + IILD 00058400 B = SX / SW 00058500 ZA(INDX) = SYD - SY * B 00058600 ZB(INDX) = SX2 - SX * B 00058700 XCDP(INDX) = SY / SW 00058800 XRNMO(INDX) = B 00058900 C 00059000 40 CONTINUE 00059100 IF (KEY .EQ. 0) GO TO 10 00059200 C 00059300 C STEP 1 COMPLETE. BEGIN PART OF STEP 2 WHICH FORMS Q IN APPENDIX A.00059400 C (NOT THE SAME AS Q HERE.) 00059500 C 00059600 K = 0 00059700 C 00059800 C LOOP THROUGH EACH TIME 00059900 C 00060000 DO 100 I = 1, NLINE 00060100 SD = 0. 00060200 KK = K + MUM 00060300 KL = K - MUM 00060400 C 00060500 C IF CDP AVERAGE RANGE NOT ZERO INITIALIZE SLIDING SUM FOR Q IN 00060600 C STEP 2 AT BEGINNING OF LINE WITH PARTIAL SUM. Q WILL BE STOPPED 00060700 C IN VARIABLE ZC. 00060800 C 00060900 IF (MUM .LE. 0) GO TO 70 00061000 C 00061100 DO 50 J = 1, MUM 00061200 50 SD = SD + ZB(J+K) 00061300 C 00061400 DO 60 J = 1, MUM 00061500 SD = SD + ZB(J+KK) 00061600 60 ZC(J+K) = SD 00061700 C 00061800 C FORM SLIDING SUM FOR 1 ON MIDDLE OF LINE IN STEP 2. 00061900 C 00062000 70 DO 80 J = MUMP, LDM 00062100 SD = SD + ZB(J+KK) 00062200 ZC(J+K) = SD 00062300 C 00062400 80 SD = SD - ZB(J+KL) 00062500 C 00062600 C IF CDP AVERAGE RANGE IS NOT ZERO FINISH OFF SLIDING SUM FOR Q 00062700 C AT THE END OF THE LINE. 00062800 C 00062900 IF (MUM .LE. 0) GO TO 100 00063000 C 00063100 DO 90 J = LDMP, LD 00063200 ZC(J+K) = SD 00063300 90 SD = SD - ZB(J+KL) 00063400 C 00063500 100 K = K + LD 00063600 C 00063700 C END STEP 2, Q CALCULATION. BEGIN STEP 3, H CALCULATION. CLEAR ZB 00063800 C WHICH WILL NOW BE USED TO STORE H. 00063900 C 00064000 CALL ARSET (ZB, LDD2, 0.) 00064100 C 00064200 C IF NUMBER OF LINES FOR RNMO AVERAGE IS ZERO FORM H BY SLIDING SUM 00064300 C FOR ENTIRE GRID. 00064400 C 00064500 IF (LUM .GT. 0) GO TO 130 00064600 C 00064700 DO 120 I = 1, LD 00064800 SD = 0. 00064900 KK = I + LUMLD 00065000 KL = I - LUMLD 00065100 JJ = LD * (LUMP-1) 00065200 C 00065300 DO 110 J = LUMP, NLNM 00065400 SD = SD + ZC(JJ+KK) 00065500 ZB(JJ+I) = SD 00065600 SD = SD - ZC(JJ+KL) 00065700 C 00065800 110 JJ = JJ + LD 00065900 C 00066000 120 CONTINUE 00066100 C 00066200 GO TO 190 00066300 C 00066400 C FOR CASES WHERE NUMBER OF LINES RNMO AVERAGE IS NOT ZERO LOOP 00066500 C THROUGH CONSTANT CDP CROSS LINES TO FORM H FOR STEP 3. 00066600 C 00066700 130 DO 180 I = 1, LD 00066800 SD = 0. 00066900 KK = I + LUMLD 00067000 KL = I - LUMLD 00067100 JJ = 0 00067200 C 00067300 C INITIALIZE SLIDING SUM FOR H AT BEGINNING OF CROSS LINE. 00067400 C 00067500 DO 140 J = 1, LUM 00067600 SD = SD + ZC(JJ+I) 00067700 140 JJ = JJ + LD 00067800 C 00067900 JJ = 0 00068000 C 00068100 DO 150 J = 1, LUM 00068200 SD = SD + ZC(JJ+KK) 00068300 ZB(JJ+I) = SD 00068400 150 JJ = JJ + LD 00068500 C 00068600 JJ = LD * (LUMP-1) 00068700 C 00068800 C FORM H FOR MIDDLE OF CROSSLINE. 00068900 C 00069000 DO 160 J = LUMP, NLNM 00069100 SD = SD + ZC(JJ+KK) 00069200 ZB(JJ+I) = SD 00069300 SD = SD - ZC(JJ+KL) 00069400 160 JJ = JJ + LD 00069500 C 00069600 JJ = LD * (NLNMP-1) 00069700 C 00069800 DO 170 J = NLNMP, NLINE 00069900 ZB(JJ+I) = SD 00070000 SD = SD - ZC(JJ+KL) 00070100 170 JJ = JJ + LD 00070200 C 00070300 180 CONTINUE 00070400 C 00070500 C END OF STEP #3 H CALCULATION. BEGIN STEP #2 P CALCULATION. 00070600 C CLEAR ZC WHICH WILL NOW BE REUSED TO STORE P. 00070700 C 00070800 190 CALL ARSET(ZC, LDD2, 0.) 00070900 K = 0 00071000 C 00071100 C IF CDP AVERAGE RANGE IS ZERO FORM P FOR ENTIRE GRID. 00071200 C 00071300 IF (MUM .GT. 0) GO TO 220 00071400 C 00071500 DO 210 I = 1, NLINE 00071600 SD = 0. 00071700 KK = K + MUM 00071800 KL = K - MUM 00071900 C 00072000 DO 200 J = MUMP, LDM 00072100 SD = SD + ZA(J+KK) 00072200 ZC(J+K) = SD 00072300 200 SD = SD - ZA(J+KL) 00072400 C 00072500 210 K = K + LD 00072600 C 00072700 GO TO 280 00072800 C 00072900 C FOR CASE WHERE CDP AVERAGE RANGE IS NOT ZERO BEGIN P CALCULATION 00073000 C BY LOOPING THROUGH EACH LINE. 00073100 C 00073200 220 DO 270 I = 1, NLINE 00073300 SD = 0. 00073400 KK = K + MUM 00073500 KL = K - MUM 00073600 C 00073700 C INITIALIZE SLIDING SUM FOR P AT BEGINNING OF LINE. 00073800 C 00073900 DO 230 J = 1, MUM 00074000 230 SD = SD + ZA(J+K) 00074100 C 00074200 DO 240 J = 1, MUM 00074300 SD = SD + ZA(J+KK) 00074400 240 ZC(J+K) = SD 00074500 C 00074600 C COMPUTE P FOR MIDDLE OF LINE. 00074700 C 00074800 DO 250 J = MUMP, LDM 00074900 SD = SD + ZA(J+KK) 00075000 ZC(J+K) = SD 00075100 250 SD = SD - ZA(J+KL) 00075200 C 00075300 C COMPUTE P FOR END OF LINE. 00075400 C 00075500 DO 260 J = LDMP, LD 00075600 ZC(J+K) = SD 00075700 260 SD = SD - ZA(J+KL) 00075800 C 00075900 270 K = K + LD 00076000 C 00076100 C END OF STEP #2 P CALCULATION. BEGIN STEP #3 G CALCULATION. 00076200 C NOTE G IS NOT STORED BUT IS USED TO IMMEDIATELY COMPLETE STEP 00076300 C 3 BY COMPUTING ALPHA AND BETA. 00076400 C 00076500 280 IF (LUM .GT. 0) GO TO 310 00076600 C 00076700 C SPECIAL CASE FOR ZERO CDP IN RNM0 AVERAGE CYCLE THRU CROSS LINES 00076800 C 00076900 DO 300 I = 1, LD 00077000 SD = 0. 00077100 KK = I + LUMLD 00077200 KL = I - LUMLD 00077300 JJ = LD * (LUMP-1) 00077400 C 00077500 C CYCLE THRU THE CDP 00077600 C 00077700 DO 290 J = LUMP, NLNM 00077800 JJPI = JJ + I 00077900 SD = SD + ZC(JJ+KK) 00078000 B = 0. 00078100 IF (DABS(ZB(JJPI)) .GT. 1.E-30) B = SD / ZB(JJPI) 00078200 C 00078300 C STORE ALPHA IN XCDP. STORE BETA IN XRNM0. 00078400 C 00078500 XCDP(JJPI) = XCDP(JJPI) - XRNMO(JJPI) * B 00078600 XRNMO(JJPI) = B 00078700 SD = SD - ZC(JJ+KL) 00078800 C 00078900 290 JJ = JJ + LD 00079000 C 00079100 300 CONTINUE 00079200 C 00079300 GO TO 370 00079400 C 00079500 C CASE WHERE THE NUMBER OF CDP IN REN0 AVERAGE IS NOT ZERO. CYCLE 00079600 C THRU THE CROSS LINES 00079700 C 00079800 310 DO 360 I = 1, LD 00079900 SD = 0. 00080000 KK = I + LUMLD 00080100 KL = I - LUMLD 00080200 JJ = 0 00080300 C 00080400 C INITIALIZE SLIDING SUM G ON CROSS LINE. 00080500 C 00080600 DO 320 J = 1, LUM 00080700 SD = SD + ZC(JJ+I) 00080800 320 JJ = JJ + LD 00080900 C 00081000 JJ = 0 00081100 C 00081200 DO 330 J = 1, LUM 00081300 JJPI = JJ + I 00081400 SD = SD + ZC(JJ+KK) 00081500 B = 0. 00081600 IF (DABS(ZB(JJPI)) .GT. 1.E-30) B = SD / ZB(JJPI) 00081700 C 00081800 C STORE ALPHA IN XCDP. STORE BETA IN XRNM0. 00081900 C 00082000 XCDP(JJPI) = XCDP(JJPI) - XRNMO(JJPI) * B 00082100 XRNMO(JJPI) = B 00082200 C 00082300 330 JJ = JJ + LD 00082400 C 00082500 JJ = LD * (LUMP-1) 00082600 C 00082700 DO 340 J = LUMP, NLNM 00082800 JJPI = JJ + I 00082900 SD = SD + ZC(JJ+KK) 00083000 B = 0. 00083100 IF (DABS(ZB(JJPI)) .GT. 1.E-30) B = SD / ZB(JJPI) 00083200 C 00083300 C STORE ALPHA IN XCDP. STORE BETA IN XRNM0. 00083400 C 00083500 XCDP(JJPI) = XCDP(JJPI) - XRNMO(JJPI) * B 00083600 XRNMO(JJPI) = B 00083700 SD = SD - ZC(JJ+KL) 00083800 C 00083900 340 JJ = JJ + LD 00084000 C 00084100 JJ = LD * (NLNMP-1) 00084200 C 00084300 C FINISH SLIDING SUM G ON END OF CROSS LINES. 00084400 C 00084500 DO 350 J = NLNMP, NLINE 00084600 JJPI = JJ + I 00084700 B = 0. 00084800 IF (DABS(ZB(JJPI)) .GT. 1.E-30) B = SD / ZB(JJPI) 00084900 C 00085000 C STORE ALPHA IN XCDP. STORE BETA IN XRNM0. 00085100 C 00085200 XCDP(JJPI) = XCDP(JJPI) - XRNMO(JJPI) * B 00085300 XRNMO(JJPI) = B 00085400 SD = SD - ZC(JJ+KL) 00085500 C 00085600 350 JJ = JJ + LD 00085700 C 00085800 360 CONTINUE 00085900 C 00086000 C BEGIN STEP #4. COMPUTE THE TRANSFORMED LAGS. 00086100 C 00086200 370 KT = 0 00086300 KK = 0 00086400 III = 0 00086500 C 00086600 C SCAN THRU THE LINES. 00086700 C 00086800 DO 400 I = 1, NLINE 00086900 KL = KT + 1 00087000 KT = LINE(I) 00087100 K = KCDP(I) 00087200 III = III + 1 00087300 IIIC = ICDP(III) 00087400 KKK = KK + K 00087500 C 00087600 C SCAN THRU THE CDP PRESENT IN EACH LINE. 00087700 C 00087800 DO 390 J = KL, KT 00087900 IF (J .LE. IIIC) GO TO 380 00088000 K = K + 1 00088100 KKK = KK + K 00088200 III = III + 1 00088300 IIIC = ICDP(III) 00088400 C 00088500 380 Q2 = Q(J) 00088600 C 00088700 C IF TRACE IS NOT REJECTED REMOVE CONSTANT AND QUADRATIC COMPONENTS 00088800 C FROM LAG VALUES. 00088900 C 00089000 IF (Q2 .LE. 1.0E6) Q(J) = Q2 - XCDP(KKK) - XRNMO(KKK)*D(J) 00089100 C 00089200 390 CONTINUE 00089300 C 00089400 400 KK = KK + LD 00089500 C 00089600 RETURN 00089700 END 00089800