CTITLESADDSB -- SUBINTERVAL DECON 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. D. THOMPSON 00000020 CA DESIGNER D. D. THOMPSON 00000030 CA LANGUAGE FORTRAN 77 00000040 CA SYSTEM IBM & CRAY 00000041 CA WRITTEN 10/01/76 00000050 C REVISED MO-DA-YR 00000060 C REVISED 06-25-84 TRA TO RETURN NORM FACTOR. 00000061 C REVISED 08-06-84 TRA TO RETURN NUM ITERATIONS. 00000062 C REVISED 05-20-85 TWH REMOVED VS FORTRAN WARNING 00000063 C MESSAGES AND ADAPTED TO RUN ON 00000064 C IBM & CRAY SYSTEMS. 00000065 C 00000070 CA 00000080 CA 00000090 CA CALL SADDSB (X,R,W,KW,IL,IH,WS,ERRMAX,COEFF1,FRACT,COEFFR, 00000100 CA MXSTPS,IPRNT,S,KEY,PR,Z,XX,V,G,GT,TZ,IBUF, 00000110 CA EBUF,IPR,VAL,NSTP) 00000120 CA 00000130 CA INPUT X = TRACE ARRAY R4 00000140 CA OUTPUT R = REFLECTION COEFFICIENT ARRAY (BLENDED R4 00000150 CA WITH NONZERO INPUT) 00000160 CA OUTPUT W = COMPUTED WAVELET (NOT USED IF KEY=0) R4 00000170 CA OUTPUT NSTP = COUNT OF NUMBER OF ITERATIONS. R4 00000171 CA INPUT KW = LENGTH OF WAVELET I4 00000180 CA INPUT IL = FIRST TRACE SAMPLE TO USE I4 00000190 CA INPUT IH = LAST TRACE SAMPLE TO USE I4 00000200 CA INPUT WS = STARTING WAVELET ARRAY R4 00000210 CA INPUT ERRMAX = MAX % ERROR R4 00000220 CA INPUT COEFF1 = STARTING COEFFICIENT R4 00000230 CA INPUT FRACT = FRACTION TO SCALE COEFFICIENT R4 00000240 CA INPUT COEFFR = COEFFICIENT TO USE FOR RC ONLY R4 00000250 CA INPUT MXSTPS = MAX NUMBER ITERATIONS I4 00000260 CA INPUT IPRNT = PRINT CODE I4 00000270 CA INPUT S = SAMPLE PERIOD OF TRACE R4 00000280 CA INPUT KEY = FLAG FOR DECON TYPE I4 00000290 CA 0 COMPUTE R ONLY (NO ITERATIONS) 00000300 CA 3 COMPUTE R AND W BY ITERATION. 00000310 CA PR = WORK ARRAY LENGTH (KZ*(KZ+1))/2 R4 00000320 CA WHERE 00000330 CA KZ = KWMAX/2 + KX/2 00000340 CA KX = IH - IL + 1 00000350 CA Z = WORK ARRAY LENGTH KZ R4 00000360 CA XX = WORK ARRAY " KX R4 00000370 CA V = WORK ARRAY " KZ R4 00000380 CA G = WORK ARRAY " KZ R4 00000390 CA GT = WORK ARRAY " KZ R4 00000400 CA TZ = WORK ARRAY " KZ R4 00000410 CA IBUF = WORK ARRAY " 5*KZ I4 00000420 CA EBUF = WORK ARRAY R4 00000430 CA LENGTH 4 * (MAX # ITERATIONS + 1) 00000440 CA INPUT IPR = PRINT UNIT NUMBER I4 00000450 CA OUTPUT VAL = NORMALIZATION FACTOR. R4 00000451 CA 00000460 CA 00000470 CA THIS ROUTINE DOES DIRECT DECONVOLUTION ON A SINGLE SUBINTERVAL. 00000480 CA IT DOES A RAMP-ON-RAMP-OFF BLENDING WITH ANY PREVIOUSLY 00000490 CA COMPUTED RC'S. GIVEN A STARTING WAVELET IT CAN EITHER ITERATE 00000500 CA TO COMPUTE BOTH WAVELET AND RC'S OR USE A STARTING WAVELET TO 00000510 CA COMPUTE RC'S ONLY. 00000520 C 00000530 C SUBROUTINES CALLED: ARMVE (S1ATP) 00000540 C SAPITR 00000550 C SADDIA 00000560 C SADDIT 00000570 C MSTRT 00000580 C 00000590 SUBROUTINE SADDSB(X,R,W,KW,IL,IH,WS,ERRMAX,COEFF1,FRACT,COEFFR, 00000600 * MXSTPS,IPRNT,S,KEY,PR,Z,XX,V,G,GT,TZ,IBUF,EBUF,IPR,VAL,NSTP) 00000610 CAEND 00000630 C 00000640 REAL X(1),R(1),W(1),WS(1),PR(1),Z(1),XX(1),V(1),G(1), 00000650 *GT(1),TZ(1) 00000660 REAL EBUF(1) 00000670 INTEGER IPRNT(5),IBUF(1) 00000680 C 00000681 IF (1 .EQ. 2) CALL S1ATP 00000682 C 00000683 KX=IH-IL+1 00000690 KR=KW/2+KX/2 00000700 KA = KW/4 + 1 00000710 KB = KR-KA 00000720 CALL ARMVE(WS,Z(KR+1),KW) 00000730 TL = S * IL 00000740 TH = S * IH 00000750 C 00000751 IF(IPRNT(5).NE.0 .OR. (KEY.NE.0 .AND. IPRNT(3).NE.0).OR. 00000770 *IPRNT(4) .NE.0) WRITE(IPR, 9000 ) TL,TH 00000780 C 00000790 C CODE FROM OLD SUBROUTINE ADJUST 00000800 C 00000810 IT=MXSTPS+1 00000820 COEFF=COEFF1/FRACT 00000830 C 00000840 10 IF(KEY.EQ.0) GO TO 70 00000850 IF(IPRNT(4).EQ.0) GO TO 20 00000860 WRITE(IPR, 9010 )(Z(I+KR),I=1,KW) 00000870 C 00000880 20 DO 60 00000890 * MAIN=1,IT 00000900 NSTP=MAIN-1 00000910 IF(MAIN.NE.1) GO TO 30 00000920 CALL MSTRT(Z,X(IL),XX,KR,KW,KX,KEY,G,GT,TZ,V) 00000930 CALL SADDIA(Z,X(IL),XX,KR,KW,KX,PE,KEY,PR,V,G,TZ, 00000940 * IPR) 00000950 DW=0. 00000960 DR=0. 00000970 GO TO 40 00000980 C 00000990 30 CALL SADDIT(Z,X(IL),XX,PE,DW,DR,COEFF,FRACT) 00001000 C 00001010 40 CONTINUE 00001020 IF(PE.LT.0.) GO TO 80 00001030 IF(IPRNT(4).EQ.0) GO TO 50 00001040 WRITE(IPR, 9020 )NSTP,PE,DW,DR,COEFF,(Z(KR+I),I=1,KW) 00001050 WRITE(IPR, 9030 )(Z(I),I=KA,KB) 00001060 C 00001070 50 IF(PE.LE.ERRMAX .OR. MAIN.EQ.IT) GO TO 80 00001080 IF(IPRNT(3).NE.0 .OR. IPRNT(5).NE.0)CALL SAPITR(Z,KR 00001090 * ,KW,PE,DR,DW,COEFF,IBUF,EBUF,0,IPRNT,IPR) 00001100 C 00001110 60 CONTINUE 00001120 C 00001130 GO TO 80 00001140 C 00001150 70 NSTP= 0 00001160 COEFF = COEFFR 00001170 CALL SADDIA(Z,X(IL),XX,KR,KW,KX,PE,1,PR,V,G,TZ,IPR) 00001180 CALL SADDIT(Z,X(IL),XX,PE,DW,DR,COEFF,FRACT) 00001190 IF(IPRNT(4).EQ.0) GO TO 80 00001200 WRITE(IPR, 9020 )NSTP,PE,DW,DR,COEFF,(Z(KR+I),I=1,KW) 00001210 WRITE(IPR, 9030 )(Z(I),I=KA,KB) 00001220 80 IF(IPRNT(3).NE.0 .OR. IPRNT(5).NE.0)CALL SAPITR(Z,KR,KW, 00001230 *PE,DR,DW,COEFF,IBUF,EBUF,1,IPRNT,IPR) 00001240 C 00001250 90 CONTINUE 00001260 C 00001270 C THIS ENDS THE CODE TAKEN FROM 00001280 C THE OLD ADJUST ROUTINE 00001290 C 00001300 LA=KW/4+1 00001310 LB=KR-LA+1 00001320 INDX= (IL-KW+1)/2-1 00001330 IF(LA+INDX.LE.0) LA=1-INDX 00001340 LAX=LA+INDX 00001350 LBX=LB+INDX 00001360 C 00001370 C COMPUTE SUM OF SQUARES FOR UNNORMALIZED WAVELET. 00001380 C 00001390 VAL=0. 00001400 C 00001410 DO 100 00001420 * I=1,KW 00001430 C 00001440 100 VAL=VAL+Z(I+KR)**2 00001450 C 00001460 VAL = SQRT(VAL) 00001470 C 00001480 C COMPUTE OVERLAP PARAMETERS LAAX, LBBX. 00001490 C 00001500 LSUM=LAX+LBX 00001510 C 00001520 DO 110 00001530 * I=LAX,LBX 00001540 IF(R(I).EQ.0.) GO TO 120 00001550 C 00001560 110 CONTINUE 00001570 C 00001580 LAAX = LBX 00001590 GO TO 130 00001600 C 00001610 120 LAAX=I 00001620 C 00001630 130 DO 140 00001640 * I=LAX,LBX 00001650 J=LSUM-I 00001660 IF(R(J).EQ.0.) GO TO 150 00001670 C 00001680 140 CONTINUE 00001690 C 00001700 150 LBBX=J 00001710 IF(LBBX-LAAX) 230 , 230 , 160 00001720 C 00001730 C LOAD UNBLENDED MIDDLE INTERVAL. 00001740 C 00001750 160 DO 170 00001760 * I=LAAX,LBBX 00001770 C 00001780 170 R(I)=Z(I-INDX)*VAL 00001790 C 00001800 IF(LAAX.EQ.LAX) GO TO 190 00001810 C 00001820 C BLEND TOP OVERLAP. 00001830 C 00001840 P=1./(LAAX-LAX) 00001850 C 00001860 DO 180 00001870 * I=LAX,LAAX 00001880 Q=(LAAX-I)*P 00001890 C 00001900 180 R(I)=VAL*(Q*R(I)+(1.-Q)*Z(I-INDX)) 00001910 C 00001920 190 IF(LBX.EQ.LBBX) GO TO 210 00001930 C 00001940 C BLEND BOTTOM OVERLAP. 00001950 C 00001960 P=1./(LBX-LBBX) 00001970 C 00001980 DO 200 00001990 * I=LBBX,LBX 00002000 Q=(I-LBBX)*P 00002010 C 00002020 200 R(I)=VAL*(Q*R(I)+(1.-Q)*Z(I-INDX)) 00002030 C 00002040 210 IF(KEY.EQ.0) RETURN 00002050 C 00002060 C NORMALIZE WAVELET AND RETURN IN W IF KEY NOT 00002070 C ZERO. 00002080 C 00002090 VAL=1./VAL 00002100 C 00002110 DO 220 00002120 * I=1,KW 00002130 C 00002140 220 W(I)=VAL*Z(KR+I) 00002150 C 00002160 RETURN 00002170 C 00002180 C BLEND EQUALLY OVER WHOLE INTERVAL. 00002190 C 00002200 230 DO 240 00002210 * I=LAX,LBX 00002220 C 00002230 240 R(I)=.5*(R(I)+Z(I-INDX)) * VAL 00002240 C 00002250 GO TO 210 00002260 C 00002270 9000 FORMAT('0',5X,'RESULTS FOR WINDOW',F10.3,' TO',F10.3) 00002280 C 00002290 9010 FORMAT('0INITIAL WAVELET'/(' ',10F13.6)) 00002300 C 00002310 9020 FORMAT('0'/'0'/'0STEP',I3,', TRACE ERROR =',F7.2,'%, WAVELET ', 00002320 *'CHANGE =',F7.2,'%, REF COEFF CHANGE =',F7.2,'%, COEFF =', 00002330 *F8.5/'0'/'0WAVELET'/(' ',10F13.6)) 00002340 C 00002350 9030 FORMAT('0'/'0REF COEFF'/(' ',10F13.6)) 00002360 END 00002370