CTITLE SARINVT -- MULTI-LAYER REFRACTION (STATICS) MODEL PREDICTION C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA CA AUTHOR J.V.S. HARVEY CA LANGUAGE VS FORTRAN (77) CA SYSTEM IBM ONLY CA REWRITTEN 18 NOV 1988 CA REVISED 21 JAN 1992 RDK MAKE ACCUMULATOR VARIABLE DTC CA REAL*8 FOR BETTER ACCURACY. CA CA CA THIS SUBROUTINE PERFORMS A MULTI-LAYER REFRACTION INVERSION CA WITH RESIDUALS FOR STATIC SHIFT ESTIMATION. CA CA CA CALL SARINVT( IPR, KPBUGF, LCTPSP, FBNLAB, CA NUMHRZ, INVPRM, EXTOPT, DXGRID, NRESID, CA DATUMF, SURFAV, CDPBEG, CDPEND, CA NVPICK, VLOCAT, VELOC1, CA NDATUM, DLOCAT, UDATUM, CA NCDPN, NCDPT, NSHOT, NRECV, NGRID, CA CDPNID, CDPNUM, CDPREF, CA SHOTID, SHOTNO, RECVID, CA ORTNSR, IOCREF, RTIMES, CA ZMODEL, PMODEL, STATIC, CA CMODEL, SMODEL, RMODEL, CA MINELV, MAXELV, MINSLW, MAXSLW, CA MTIMES, SRESID, RRESID, CA V1AVG, REPVEL, PSEUDO, CA NWORK, MODTMP, LSFOLD, WORKAR, CA CPUTIM, CPUDEL, ERR1 ) CA CA CA IN/OUT ARGUMENT TYPE DESCRIPTION CA CA IN IPR I4 DEBUG OUTPUT UNIT NUMBER CA IN KPBUGF I4 DEBUG FLAG CA AT LEAST 0 FOR DUMP OF SCALARS CA AT LEAST 2 FOR DUMP OF COORDINATES CA IN LCTPSP I4 NUMBER OF TRACES PER SHOT CA IN FBNLAB CH8 FIRST-BREAK NODE LABELS CA CA CA IN NUMHRZ I4 NUMBER OF HORIZONS CA CA UPDATE INVPRM I2 HORIZON INVERSION PARAMETERS CA (HRZ,01) = NUMBER OF PASSES CA (HRZ,02) = DIAGNOSTICS LEVEL CA (HRZ,03) = CDP WINDOW WIDTH CA (HRZ,04) = OFFSET WINDOW WIDTH CA (HRZ,05) = VELOCITY SMOOTHING LENGTH CA (HRZ,06) = CRIT. OFFSET SMOOTHING LENGTH CA (HRZ,07) = DIP ANGLE SMOOTHING LENGTH CA (HRZ,08) = ELEVATION SMOOTHING LENGTH CA (HRZ,09) = LEFT END REGRESSION LENGTH CA (HRZ,10) = RIGHT END REGRESSION LENGTH CA ( 2-D ARRAY DIMENSIONED: 8 BY 10 ) CA CA IN EXTOPT CH4 EXTRAPOLATION OPTIONS CA (HRZ,1) = LEFT END OPTION CA (HRZ,2) = RIGHT END OPTION CA ( 2-D ARRAY DIMENSIONED: 8 BY 2 ) CA CA IN DXGRID R4 ANALYSIS GRID SPACING CA CA IN NRESID I4 NUMBER OF RESIDUAL ITERATIONS CA ( 0 DISABLES RESIDUALS ) CA CA CA IN DATUMF CH4 PROCESS DATUM OPTION FLAG CA 'FCDP' = FLOATING DATUM CA 'MINS' = AVG. FLOATING DATUM CA 'USER' = USER-DEFINED DATUM CA CA 'MIGD' = MIGRATION DATUM CA 'REFL' = REFERENCE REFLECTOR FOCUSING CA 'FLAT' = REFERENCE REFLECTOR FLATTENING CA CA IN SURFAV R4 AVERAGE SURFACE ELEVATION CA IN CDPBEG I4 STARTING CDP FOR REPLACEMENT VELOCITY CA 0 = DEFAULT FOR NO LIMITS CA IN CDPEND I4 ENDING CDP FOR REPLACEMENT VELOCITY CA -99999 = DEFAULT FOR NO LIMITS CA CA CA IN NVPICK I4 NUMBER OF CURRENT VELOCITY PICKS CA IN VLOCAT I2 LOCATION INDICES OF VELOCITY PICKS CA ( 1-D ARRAY DIMENSIONED: 24 ) CA IN VELOC1 R4 SURFACE VELOCITY (V1) FUNCTION CA ( 1-D ARRAY DIMENSIONED: 24 ) CA CA IN NDATUM I4 NUMBER OF CURRENT DATUM PICKS CA IN DLOCAT I2 LOCATION INDICES OF DATUM PICKS CA ( 1-D ARRAY DIMENSIONED: 24 ) CA IN UDATUM R4 USER-DEFINED DATUM FUNCTION CA ( 1-D ARRAY DIMENSIONED: 24 ) CA CA CA IN NCDPN I4 NUMBER OF CDPN LOCATIONS CA IN NCDPT I4 NUMBER OF TRACES PER CDPN CA IN NSHOT I4 NUMBER OF SHOTS PICKED CA IN NRECV I4 NUMBER OF RECEIVER POSITIONS FOUND CA IN NGRID I4 NUMBER OF GRID X LOCATIONS CA CA CA IN CDPNID R4 CDP IDENTIFICATION ARRAY CA (1,I) = X COORDINATE (IN-LINE) CA IN UNITS OF GRID-POINT POSITION CA (2,I) = SURFACE ELEVATION CA IN ORIGINAL UNITS CA ( 2-D ARRAY DIMENSIONED: 2 BY NCDPN ) CA IN CDPNUM I2 CDP NUMBERS (LABELS) CA ( 1-D ARRAY DIMENSIONED: NCDPN ) CA CA IN CDPREF I2 CDP-TO-SHOT CROSS-REFERENCE INDICES CA (1,SIDE,CDPT,CDPN) = SHOT INDEX CA (2,SIDE,CDPT,CDPN) = RECEIVER INDEX CA (3,SIDE,CDPT,CDPN) = ORTN (TRACE) CA ( 4-D ARRAY DIMENSIONED: 3 BY 2 BY NCDPT BY NCDPN ) CA CA CA IN SHOTID R4 SHOTPOINT IDENTIFICATION ARRAY CA (1,I) = X COORDINATE (IN-LINE) CA IN UNITS OF GRID-POINT POSITION CA (2,I) = Y COORDINATE (CROSS-LINE) CA IN UNITS OF GRID-POINT POSITION CA (3,I) = SURFACE ELEVATION CA IN ORIGINAL UNITS CA (4,I) = RECIPROCAL RECEIVER NUMBER CA AS INDEX PLUS FRAC. DISTANCE CA (5,I) = EQUIVALENT CDPN LOCATION CA AS INDEX PLUS FRAC. DISTANCE CA ( 2-D ARRAY DIMENSIONED: 5 BY NSHOT ) CA CA IN SHOTNO I2 PICKED SHOTPOINT NUMBER LIST CA ( 1-D ARRAY DIMENSIONED: NSHOT ) CA CA CA IN RECVID R4 RECEIVER IDENTIFICATION ARRAY CA (1,I) = X COORDINATE (IN-LINE) CA IN UNITS OF GRID-POINT POSITION CA (2,I) = Y COORDINATE (CROSS-LINE) CA IN UNITS OF GRID-POINT POSITION CA (3,I) = SURFACE ELEVATION CA IN ORIGINAL UNITS CA (4,I) = RECIPROCAL SHOT LOCATION NUMBER CA AS INDEX PLUS FRAC. DISTANCE CA (5,I) = EQUIVALENT CDPN LOCATION CA AS INDEX PLUS FRAC. DISTANCE CA ( 2-D ARRAY DIMENSIONED: 5 BY NRECV ) CA CA CA IN ORTNSR I2 TRACE CROSS-REFERENCE FOR SHOT AND RECEIVER CA ORTNSR(I,J) = TRACE NUMBER (ORTN) CA OF SHOT I AT REC J CA = 0 IF OUT OF RANGE CA ( 2-D ARRAY DIMENSIONED: NSHOT BY NRECV ) CA CA IN IOCREF I2 REFERENCE CDPT LOCATIONS FOR INWARD AND CA OUTWARD CANCELLATION CA (1,SIDE,HORIZ,CDPN) = INWARD REFERENCE CA FIRST BREAKS ONLY CA (DEFAULTS TO -999) CA (2,SIDE,HORIZ,CDPN) = OUTWARD REFERENCE CA FIRST BREAKS ONLY CA (DEFAULTS TO -999) CA (3,SIDE,HORIZ,CDPN) = OUTWARD REFERENCE CA INCLUDING PHANTOMS CA (DEFAULTS TO -999) CA ( 4-D ARRAY DIMENSIONED: 3 BY 2 BY NUMHRZ BY NCDPN ) CA CA IN RTIMES R4 REMAPPED REFRACTION TIMES IN SECONDS CA DEFAULTED TO ZERO FOR NO PICKED TIME CA ( 3-D ARRAY DIMENSIONED: LCTPSP BY NSHOT BY NUMHRZ ) CA CA CA OUT ZMODEL R4 THE ELEVATION MODEL FOR ALL HORIZONS CA (GRID,CURHRZ,1) = ELEVATION ESTIMATE CA (GRID,CURHRZ,2) = TANGENT OF DIP CA ( 3-D ARRAY DIMENSIONED: NGRID BY 0:NUMHRZ BY 2 ) CA CA OUT PMODEL R4 THE SLOWNESS MODEL FOR THE LAYER BELOW CA EACH HORIZON CA (GRID,CURHRZ,1) = SLOWNESS ESTIMATE CA (GRID,CURHRZ,2) = VELOCITY RMS ERROR CA ( 3-D ARRAY DIMENSIONED: NGRID BY 0:NUMHRZ BY 2 ) CA CA OUT STATIC R4 MODEL-PREDICTED RAYPATH STATIC SHIFT IN SEC CA *** DEFAULTS TO 1.0E+10 FOR NO ESTIMATE *** CA ( 2-D ARRAY DIMENSIONED: LCTPSP BY NSHOT ) CA CA CA OUT CMODEL R8 MODEL VALUES ASSOCIATED WITH CDP LOCATIONS CA (CDPN,1) = VERTICAL RAY TIME (SEC) CA (CDPN,2) = DATUM ELEVATION CA (CDPN,3) = STATIC REFERENCE TIME CA (CDPN,4) = "FLAT" VERTICAL RAY TIME (SEC) CA ( 2-D ARRAY DIMENSIONED: NCDPN BY 9 ) CA CA OUT SMODEL R8 SHOT-CONSISTENT OVERBURDEN VALUES CA (SHOT,1) = VERTICAL RAY TIMES CA ( 2-D ARRAY DIMENSIONED: NSHOT BY 5 ) CA CA OUT RMODEL R8 RECV-CONSISTENT OVERBURDEN VALUES CA (RECV,1) = VERTICAL RAY TIMES CA ( 2-D ARRAY DIMENSIONED: NRECV BY 5 ) CA CA CA OUT MINELV R4 ELEVATION MODEL MINIMA (PER HORIZON) CA ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) CA OUT MAXELV R4 ELEVATION MODEL MAXIMA (PER HORIZON) CA ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) CA CA OUT MINSLW R4 SLOWNESS MODEL MINIMA (PER HORIZON) CA ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) CA OUT MAXSLW R4 SLOWNESS MODEL MAXIMA (PER HORIZON) CA ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) CA CA CA OUT MTIMES R4 MODEL-PREDICTED TIMES IN SECONDS CA SIGN INDICATES SOURCE: CA POSITIVE FOR MODELLED PICK CA ZERO FOR MISSING PICK CA NEGATIVE FOR PHANTOMED PICK CA ABSOLUTE VALUE GIVES PROPER TIME CA ( 3-D ARRAY DIMENSIONED: LCTPSP BY NSHOT BY NUMHRZ ) CA CA OUT SRESID R4 SHOT-CONSISTENT LEAST-SQUARES RESIDUALS CA (SHOT,1) = POSITIVE RAY (L-R) RESIDUAL CA (SHOT,2) = NEGATIVE RAY (R-L) RESIDUAL CA ( 2-D ARRAY DIMENSIONED: NSHOT BY 2 ) CA CA OUT RRESID R4 RECV-CONSISTENT OVERBURDEN VALUES CA (RECV,1) = POSITIVE RAY (L-R) RESIDUAL CA (RECV,2) = NEGATIVE RAY (R-L) RESIDUAL CA ( 2-D ARRAY DIMENSIONED: NRECV BY 2 ) CA CA CA OUT V1AVG R4 AVERAGE SURFACE VELOCITY CA CA OUT REPVEL R4 REPLACEMENT VELOCITY FOR STATICS CALC. CA *** NEAR-ZERO VALUE INDICATES ERROR *** CA OUT PSEUDO R4 PSEUDO-DATUM ELEVATION CA CA CA IN NWORK I4 LENGTH OF WORK VECTORS CA ( -- MAX OF: NCDPN, NSHOT, NRECV -- CA CA WORK MODTMP R8 TEMPORARY MODEL VALUES CA ( 2-D ARRAY DIMENSIONED: NWORK BY 8 ) CA WORK LSFOLD I4 FOLD/DIAGNOSTIC FLAGS CA ( 1-D ARRAY DIMENSIONED: NWORK ) CA CA WORK WORKAR R8 WORK AREA FOR SMOOTHING CA ( 1-D ARRAY DIMENSIONED: NWORK BY 5 ) CA CA CA UPDATE CPUTIM I4 CUMMULATIVE CPU TIME CA UPDATE CPUDEL I4 RELATIVE CPU TIME CA CA OUT ERR1 I4 ERROR FLAG CA -1 WARNINGS ISSUED CA 0 FOR NO ERROR CA 1 INVERSION ERROR CA 2 INVALID DATUM POSITION CA 3 INVALID V1 POSITION CA CA 4 V1 EXTRAPOLATION BUG CA 5 Z1 EXTRAPOLATION BUG CA CA********************************************************************** CA** **** CA** NOTE: THE FOLLOWING EVALUATIONS GIVE STATIC SHIFTS **** CA** **** CA** SHOT SHIFT = CMODEL(CDPN,3) - SMODEL(SHOT,1) **** CA** RECV SHIFT = CMODEL(CDPN,3) - RMODEL(RECV,1) **** CA** **** CA********************************************************************** C C SUBROUTINES CALLED: SARGRMR -- GRM INVERSION WITH RESIDUALS C SARFILL -- INTERPOLATION/SMOOTHING/EXTRAPOLATION C C ARSET -- ARRAY BLOCK FILL (W/ SCALAR) C CLOCK -- CPU TIME C C*********************************************************************** C SUBROUTINE SARINVT( IPR, KPBUGF, LCTPSP, FBNLAB, * NUMHRZ, INVPRM, EXTOPT, DXGRID, NRESID, * DATUMF, SURFAV, CDPBEG, CDPEND, * NVPICK, VLOCAT, VELOC1, * NDATUM, DLOCAT, UDATUM, * NCDPN, NCDPT, NSHOT, NRECV, NGRID, * CDPNID, CDPNUM, CDPREF, * SHOTID, SHOTNO, RECVID, * ORTNSR, IOCREF, RTIMES, * ZMODEL, PMODEL, STATIC, * CMODEL, SMODEL, RMODEL, * MINELV, MAXELV, MINSLW, MAXSLW, * MTIMES, SRESID, RRESID, * V1AVG, REPVEL, PSEUDO, * NWORK, MODTMP, LSFOLD, WORKAR, * CPUTIM, CPUDEL, ERR1 ) IMPLICIT INTEGER (A-Z) C CHARACTER*8 FBNLAB C INTEGER*2 INVPRM(8,10) CHARACTER*4 EXTOPT(8,2) REAL DXGRID C CHARACTER*4 DATUMF REAL SURFAV C INTEGER*2 VLOCAT(24) REAL VELOC1(24) INTEGER*2 DLOCAT(24) REAL UDATUM(24) C REAL CDPNID(2,NCDPN) INTEGER*2 CDPNUM(NCDPN) INTEGER*2 CDPREF(3,2,NCDPT,NCDPN) C REAL SHOTID(5,NSHOT) INTEGER*2 SHOTNO(NSHOT) REAL RECVID(5,NRECV) C INTEGER*2 ORTNSR(NSHOT,NRECV) INTEGER*2 IOCREF(3,2,NUMHRZ,NCDPN) REAL RTIMES(LCTPSP,NSHOT,NUMHRZ) C REAL ZMODEL(NGRID,0:NUMHRZ,2) REAL PMODEL(NGRID,0:NUMHRZ,2) C REAL STATIC(LCTPSP,NSHOT) C REAL*8 CMODEL(NCDPN,9) REAL*8 SMODEL(NSHOT,5) REAL*8 RMODEL(NRECV,5) C REAL MAXELV(0:NUMHRZ) REAL MINELV(0:NUMHRZ) C REAL MAXSLW(0:NUMHRZ) REAL MINSLW(0:NUMHRZ) C REAL MTIMES(LCTPSP,NSHOT,NUMHRZ) REAL SRESID(NSHOT,2) REAL RRESID(NRECV,2) C REAL V1AVG REAL REPVEL REAL PSEUDO C REAL*8 MODTMP(NWORK,8) INTEGER LSFOLD(NWORK) REAL*8 WORKAR(NWORK,5) C---------------------------------------------------------------------- C C INTERNAL CONSTANTS C REAL UNDEFR REAL*8 UNDEFV PARAMETER ( UNDEFR = +1.0E+10 ) PARAMETER ( UNDEFV = +1.0D+10 ) C---------------------------------------------------------------------- C C LOCAL SCALAR VARIABLES C REAL COUNT REAL DP REAL DT REAL*8 DTC REAL DTL REAL DTR REAL DTSUM REAL DX REAL DXL REAL DXR REAL DX1 REAL DX2 REAL DZ REAL MAXVEL REAL MAXVRT REAL MINVEL REAL MINVRT REAL PA REAL PB REAL PC REAL PNL REAL PNR REAL PX REAL P1 REAL P1AVG REAL SLOPE REAL SLOW REAL T0AVG REAL V REAL VX REAL X REAL XA REAL XB REAL XL REAL XR REAL X1 REAL X2 REAL X3 REAL ZA REAL ZB REAL ZX REAL Z0 REAL Z1 REAL Z1L REAL Z1R REAL Z2 REAL Z2L REAL Z2R REAL Z3 C REAL*8 P0 C ------------------------------------------- C C LOCAL INTEGER CONSTANTS C INTEGER*2 ITEMP C---------------------------------------------------------------------- C C CHARACTER ARRAYS C CHARACTER*8 HRZNUM /'12345678'/ C CHARACTER*132 OUTPUT C---------------------------------------------------------------------- C C FUNCTIONS C REAL ABS REAL FLOAT REAL SQRT C C*********************************************************************** C*** **** C*** INITIALZATION **** C*** **** C*********************************************************************** C C DISPLAY INVERSION PARAMETERS C CALL CLOCK( CPUDEL, CPUTIM ) C WRITE( IPR, 8000 ) NUMHRZ, CPUTIM, ( HRZNUM(I:I), I=1,NUMHRZ ) C WRITE( IPR, 8001 ) ( INVPRM(I,1), I=1,NUMHRZ ) WRITE( IPR, 8002 ) ( INVPRM(I,2), I=1,NUMHRZ ) WRITE( IPR, 8003 ) ( INVPRM(I,3), I=1,NUMHRZ ) WRITE( IPR, 8004 ) ( INVPRM(I,4), I=1,NUMHRZ ) WRITE( IPR, 8005 ) ( INVPRM(I,5), I=1,NUMHRZ ) WRITE( IPR, 8006 ) ( INVPRM(I,6), I=1,NUMHRZ ) WRITE( IPR, 8007 ) ( INVPRM(I,7), I=1,NUMHRZ ) WRITE( IPR, 8008 ) ( INVPRM(I,8), I=1,NUMHRZ ) C WRITE( IPR, 8011 ) ( EXTOPT(I,1), I=1,NUMHRZ ) WRITE( IPR, 8012 ) ( INVPRM(I,9), I=1,NUMHRZ ) C WRITE( IPR, 8021 ) ( EXTOPT(I,2), I=1,NUMHRZ ) WRITE( IPR, 8022 ) ( INVPRM(I,10), I=1,NUMHRZ ) C WRITE( IPR, 8030 ) DXGRID, SURFAV C====================================================================== C C INITIALIZE FOLD C CALL ARSET( LSFOLD, NWORK, 1 ) C C*********************************************************************** C*** **** C*** SURFACE ATTRIBUTES **** C*** **** C*********************************************************************** C C INITIALIZE SLOWNESS SCATTER FOR SURFACE C DO 5 GRID = 1, NGRID PMODEL(GRID,0,1) = UNDEFR PMODEL(GRID,0,2) = 0.0 5 CONTINUE C====================================================================== C C CONSTANT VELOCITY C IF( NVPICK .EQ. 1 ) THEN V1AVG = VELOC1(1) P1AVG = 1.0/V1AVG C DO 10 SHOT = 1, NSHOT SMODEL(SHOT,1) = P1AVG 10 CONTINUE C DO 20 RECV = 1, NRECV RMODEL(RECV,1) = P1AVG 20 CONTINUE C DO 30 CDPN = 1, NCDPN CMODEL(CDPN,1) = P1AVG 30 CONTINUE C DO 40 GRID = 1, NGRID PMODEL(GRID,0,1) = P1AVG 40 CONTINUE C====================================================================== C C SORT THE VELOCITY PICKS C ELSE IF( NVPICK .GT. 0 ) THEN DO 50 RPICK = 2, NVPICK PICKA = 1 C LA = VLOCAT(PICKA) IF( LA .GT. 0 ) THEN XA = SHOTID(1,LA) ELSE IF( LA .LT. 0 ) THEN LA = 0 - LA XA = RECVID(1,LA) ELSE ERR1 = 3 RETURN ENDIF C DO 45 PICKB = 2, NVPICK LB = VLOCAT(PICKB) IF( LB .GT. 0 ) THEN XB = SHOTID(1,LB) ELSE IF( LB .LT. 0 ) THEN LB = 0 - LB XB = RECVID(1,LB) ELSE ERR1 = 3 RETURN ENDIF C IF( XB .LT. XA ) THEN ITEMP = VLOCAT(PICKA) VLOCAT(PICKA) = VLOCAT(PICKB) VLOCAT(PICKB) = ITEMP C VX = VELOC1(PICKA) VELOC1(PICKA) = VELOC1(PICKB) VELOC1(PICKB) = VX ELSE XA = XB ENDIF C PICKA = PICKB 45 CONTINUE 50 CONTINUE C---------------------------------------------------------------------- C C INTERPOLATE THE VELOCITY PICKS C LA = VLOCAT(1) IF( LA .GT. 0 ) THEN XA = SHOTID(1,LA) ELSE LA = 0 - LA XA = RECVID(1,LA) ENDIF PA = 1.0/VELOC1(1) PICKA = 1 V1AVG = VELOC1(1) C DO 95 PICKB = 2, NVPICK V1AVG = V1AVG + VELOC1(PICKB) C LB = VLOCAT(PICKB) IF( LB .GT. 0 ) THEN XB = SHOTID(1,LB) ELSE LB = 0 - LB XB = RECVID(1,LB) ENDIF PB = 1.0/VELOC1(PICKB) C DP = ( PB - PA ) DX = ( XB - XA ) IF( ABS( DX ) .GT. 1.0E-10 ) THEN SLOPE = DP/DX PC = PA - SLOPE*XA C ------------------------------------------- C C SURFACE SLOWNESS AT SHOTPOINTS C DO 60 SHOT = 1, NSHOT X = SHOTID(1,SHOT) P0 = UNDEFV C IF( X .EQ. XA ) THEN P0 = PA C ELSE IF( X .EQ. XB ) THEN P0 = PB C ELSE IF( XA .LT. X .AND. X .LT. XB ) THEN P0 = PC + SLOPE*X ENDIF C SMODEL(SHOT,1) = P0 60 CONTINUE C ------------------------------------------- C C SURFACE SLOWNESS AT RECEIVERS C DO 70 RECV = 1, NRECV X = RECVID(1,RECV) P0 = UNDEFV C IF( X .EQ. XA ) THEN P0 = PA C ELSE IF( X .EQ. XB ) THEN P0 = PB C ELSE IF( XA .LT. X .AND. X .LT. XB ) THEN P0 = PC + SLOPE*X ENDIF C RMODEL(RECV,1) = P0 70 CONTINUE C ------------------------------------------- C C SURFACE SLOWNESS AT DEPTHPOINTS C DO 80 CDPN = 1, NCDPN X = CDPNID(1,CDPN) P0 = UNDEFV C IF( X .EQ. XA ) THEN P0 = PA C ELSE IF( X .EQ. XB ) THEN P0 = PB C ELSE IF( XA .LT. X .AND. X .LT. XB ) THEN P0 = PC + SLOPE*X ENDIF C CMODEL(CDPN,1) = P0 80 CONTINUE C ------------------------------------------- C C SURFACE SLOWNESS FOR MODEL C DO 90 GRID = 1, NGRID X = GRID*1.0 P1 = UNDEFR C IF( X .EQ. XA ) THEN P1 = PA C ELSE IF( X .EQ. XB ) THEN P1 = PB C ELSE IF( XA .LT. X .AND. X .LT. XB ) THEN P1 = PC + SLOPE*X ENDIF C PMODEL(GRID,0,1) = P1 90 CONTINUE ENDIF C PICKA = PICKB PA = PB XA = XB 95 CONTINUE C V1AVG = V1AVG/NVPICK ELSE ERR1 = 3 RETURN ENDIF C---------------------------------------------------------------------- C C FILL HOLES IN SLOWNESS C C IN/OUT ARGUMENT TYPE DESCRIPTION C C IN INTOPT CH4 INTERNAL HOLE FILLING OPTION C 'SKIP' FOR SKIP HOLES C 'FILL' FILL HOLES (LINEARLY) C IN UNDEFV R8 VALUE OF HOLE ("NOT COMPUTED") C C C IN SMFLAG CH4 SMOOTHING OPTION (OVER WINDOW) C 'SKIP' FOR NO SMOOTHING C C 'WSUM' WEIGHTED SUM C 'WLSF' WEIGHTED LEAST-SQUARES SMOOTHING C 'WAVG' WEIGHTED AVERAGING SMOOTHING C C 'ULSF' UNWEIGHTED LEAST-SQUARES SMOOTHING C 'GRAD' UNWEIGHTED LEAST-SQUARES SLOPE C C IN SMTHL I4 SMOOTHING WINDOW LENGTH (TOTAL) C *** MUST BE AN ODD NUMBER < NGRID *** C C C IN LEXOPT CH4 LEFT EXTENSION OPTION C 'NONE' FOR NO EXTRAPOLATION C 'CONS' FOR CONSTANT VALUE C 'ELEV' FOR CONSTANT VALUE W/ TAPER C 'DIP ' FOR CONSTANT SLOPE C IN LREGRL I4 LEFT REGRESSION LENGTH C C IN REXOPT CH4 RIGHT EXTENSION OPTION C 'NONE' FOR NO EXTRAPOLATION C 'CONS' FOR CONSTANT VALUE C 'ELEV' FOR CONSTANT VALUE W/ TAPER C 'DIP ' FOR CONSTANT SLOPE C IN RREGRL I4 RIGHT REGRESSION LENGTH C C C IN NGRID I4 NUMBER OF GRID POINTS IN MODEL C C IN RMFOLD I4 MODEL FOLD (WEIGHTING PER GRID POINT) C ( 1-D ARRAY DIMENSIONED: NGRID ) C C UPDATE RMODEL R8 MODEL ARRAY (OR ITS SLOPE) C ( 1-D ARRAY DIMENSIONED: NGRID ) C C WORK WORKAR R8 WORK AREA FOR SMOOTHING C ( 2-D ARRAY DIMENSIONED: NGRID BY 5 ) C C CALL SARFILL( 'FILL', UNDEFV, 'SKIP', 0, * 'ELEV', 0, 'ELEV', 0, * NSHOT, LSFOLD, SMODEL, WORKAR ) C CALL SARFILL( 'FILL', UNDEFV, 'SKIP', 0, * 'ELEV', 0, 'ELEV', 0, * NRECV, LSFOLD, RMODEL, WORKAR ) C CALL SARFILL( 'FILL', UNDEFV, 'SKIP', 0, * 'ELEV', 0, 'ELEV', 0, * NCDPN, LSFOLD, CMODEL, WORKAR ) C ------------------------------------------- C C EXTEND SLOWNESS ON REFERENCE GRID C FIRST = 0 LAST = NGRID C DO 125 GRID = 1, NGRID IF( PMODEL(GRID,0,1) .NE. UNDEFR ) THEN LAST = GRID IF( FIRST .EQ. 0 ) FIRST = GRID ENDIF 125 CONTINUE C IF( FIRST .GT. 1 ) THEN PX = PMODEL(FIRST,0,1) C JX = FIRST - 1 C DO 150 GRID = 1, JX PMODEL(GRID,0,1) = PX 150 CONTINUE ENDIF C IF( LAST .LT. NGRID ) THEN PX = PMODEL(LAST,0,1) C JX = LAST + 1 C DO 175 GRID = JX, NGRID PMODEL(GRID,0,1) = PX 175 CONTINUE ENDIF C ------------------------------------------- C C COMPLETE CDP MODEL DEFINITION AT SURFACE C DO 190 CDPN = 1, NCDPN CMODEL(CDPN,2) = 0.0 CMODEL(CDPN,3) = CDPNID(2,CDPN) CMODEL(CDPN,4) = 0.0 190 CONTINUE C====================================================================== C C SET SURFACE DIP TO ZERO C DO 210 GRID = 1, NGRID ZMODEL(GRID,0,1) = UNDEFR ZMODEL(GRID,0,2) = 0.0 210 CONTINUE C ------------------------------------------- C C INTIALIZE THE MODEL TO THE SURFACE ELEVATION C FIRST = 0 LAST = NGRID C X = 0.0 DO 280 GRID = 1, NGRID X = X + 1.0 C ------------------------------------------- C C FIT CDP SURFACE VALUES FIRST C XA = CDPNID(1,1) ZA = CDPNID(2,1) C DO 225 CDPN = 2, NCDPN XB = CDPNID(1,CDPN) ZB = CDPNID(2,CDPN) C DZ = ( ZB - ZA ) DX = ( XB - XA ) C IF( ABS( DX ) .GT. 1.0E-10 ) THEN SLOPE = DZ/DX C IF( XA .LE. X .AND. X .LE. XB ) THEN ZMODEL(GRID,0,1) = ZA + SLOPE*( X - XA ) ENDIF C ZA = ZB XA = XB ENDIF 225 CONTINUE C ------------------------------------------- C C FIT RECV SURFACE VALUES SECOND C IF( ZMODEL(GRID,0,1) .EQ. UNDEFR ) THEN XA = RECVID(1,1) ZA = RECVID(3,1) C DO 250 RECV = 2, NRECV XB = RECVID(1,RECV) ZB = RECVID(3,RECV) C DZ = ( ZB - ZA ) DX = ( XB - XA ) C IF( ABS( DX ) .GT. 1.0E-10 ) THEN SLOPE = DZ/DX C IF( XA .LE. X .AND. X .LE. XB ) THEN ZMODEL(GRID,0,1) = ZA + SLOPE*( X - XA ) ENDIF C ZA = ZB XA = XB ENDIF 250 CONTINUE C ------------------------------------------- C C FIT SHOT SURFACE VALUES LAST C IF( ZMODEL(GRID,0,1) .EQ. UNDEFR ) THEN IF( X .LE. SHOTID(1,1) ) THEN ZMODEL(GRID,0,1) = SHOTID(3,1) C ELSE IF( X .GE. SHOTID(1,NSHOT) ) THEN ZMODEL(GRID,0,1) = SHOTID(3,NSHOT) C ELSE XA = SHOTID(1,1) ZA = SHOTID(3,1) C DO 275 SHOT = 2, NSHOT XB = SHOTID(1,SHOT) ZB = SHOTID(3,SHOT) C DZ = ( ZB - ZA ) DX = ( XB - XA ) C IF( ABS( DX ) .GT. 1.0E-10 ) THEN SLOPE = DZ/DX C IF( XA .LE. X .AND. X .LE. XB ) THEN ZMODEL(GRID,0,1) = ZA + SLOPE*( X - XA ) ENDIF C ZA = ZB XA = XB ENDIF 275 CONTINUE ENDIF ENDIF ENDIF C ------------------------------------------- C C NOTE FIRST AND LAST DEFINED LOCATIONS C IF( ZMODEL(GRID,0,1) .NE. UNDEFR ) THEN LAST = GRID IF( FIRST .EQ. 0 ) FIRST = GRID ENDIF 280 CONTINUE C---------------------------------------------------------------------- C C EXTEND ELEVATION C IF( FIRST .GT. 1 ) THEN ZX = ZMODEL(FIRST,0,1) C JX = FIRST - 1 C DO 285 GRID = 1, JX ZMODEL(GRID,0,1) = ZX 285 CONTINUE ENDIF C IF( LAST .LT. NGRID ) THEN ZX = ZMODEL(LAST,0,1) C JX = LAST + 1 C DO 290 GRID = JX, NGRID ZMODEL(GRID,0,1) = ZX 290 CONTINUE ENDIF C====================================================================== C C DEFINE LIMITS C MINSLW(0) = 1.0E+10 MAXSLW(0) = 0.0 C MINELV(0) = 1.0E+10 MAXELV(0) = 0.0 C DO 295 GRID = 1, NGRID P1 = PMODEL(GRID,0,1) Z1 = ZMODEL(GRID,0,1) C IF( P1 .EQ. UNDEFR ) THEN ERR1 = 4 ELSE IF( MINSLW(0) .GT. P1 ) MINSLW(0) = P1 IF( MAXSLW(0) .LT. P1 ) MAXSLW(0) = P1 ENDIF C IF( Z1 .EQ. UNDEFR ) THEN ERR1 = 5 ELSE IF( MINELV(0) .GT. Z1 ) MINELV(0) = Z1 IF( MAXELV(0) .LT. Z1 ) MAXELV(0) = Z1 ENDIF 295 CONTINUE C C*********************************************************************** C*** **** C*** INVERSION **** C*** **** C*********************************************************************** C C IN/OUT ARGUMENT TYPE DESCRIPTION C C IN IPR I4 PRINT UNIT NUMBER C IN KPBUGF I4 DEBUG PRINT FLAG C C IN DATUMF CH4 STATICS APPLICATION (DATUM) TYPE C 'USER' FOR USER-DEFINED DATUM C 'FCDP' FLOATING CDP DATUM C 'MINS' AVG. FLOATING CDP DATUM C 'MIGD' MIGRATION DATUM C 'REFL' REFERENCE REFLECTOR FOCUSING C 'FLAT' REFERENCE REFLECTOR FLATTENING C C C UPDATE INVPRM I2 HORIZON INVERSION PARAMETERS C (HRZ,01) = NUMBER OF PASSES C (HRZ,02) = DIAGNOSTICS LEVEL C (HRZ,03) = CDP WINDOW WIDTH C (HRZ,04) = OFFSET WINDOW WIDTH C (HRZ,05) = VELOCITY SMOOTHING LENGTH C (HRZ,06) = CRIT. OFFSET SMOOTHING LENGTH C (HRZ,07) = DIP ANGLE SMOOTHING LENGTH C (HRZ,08) = ELEVATION SMOOTHING LENGTH C (HRZ,09) = LEFT END REGRESSION LENGTH C (HRZ,10) = RIGHT END REGRESSION LENGTH C ( 2-D ARRAY DIMENSIONED: 8 BY 10 ) C C IN EXTOPT CH4 EXTRAPOLATION OPTIONS C (HRZ,1) = LEFT END OPTION C (HRZ,2) = RIGHT END OPTION C ( 2-D ARRAY DIMENSIONED: 8 BY 2 ) C C IN DXGRID R4 GRID SPACING C C IN NRESID I4 NUMBER OF RESIDUAL ITERATIONS C ( 0 DISABLES RESIDUALS ) C C C IN NUMHRZ I4 NUMBER OF HORIZONS IN ANALYSIS C IN NSHOT I4 NUMBER OF SHOTPOINTS PICKED C IN NRECV I4 NUMBER OF UNIQUE RECEIVER POSITIONS C IN NCDPN I4 NUMBER OF CDP LOCATIONS C IN NGRID I4 NUMBER OF GRID X LOCATIONS C IN LCTPSP I4 MAX. NUMBER OF TRACES PER SHOT (LCTPSP) C IN NCDPT I4 MAX. NUMBER OF TRACES PER CDPN SIDE C C IN SHOTNO I2 PICKED SHOTPOINT NUMBER LIST C ( 1-D ARRAY DIMENSIONED: NSHOT ) C C IN CDPNUM I2 CDP NUMBERS (LABELS) C ( 1-D ARRAY DIMENSIONED: NCDPN ) C C IN ORTNSR I2 TRACE CROSS-REFERENCE FOR SHOT AND RECEIVER C ORTNSR(I,J) = TRACE NUMBER (ORTN) C OF SHOT I AT REC J C = 0 IF OUT OF RANGE C ( 2-D ARRAY DIMENSIONED: NSHOT BY NRECV ) C C C IN SHOTID R4 SHOTPOINT IDENTIFICATION ARRAY C (1,I) = X COORDINATE (IN-LINE) C IN UNITS OF GRID-POINT POSITION C (2,I) = Y COORDINATE (CROSS-LINE) C IN UNITS OF GRID-POINT POSITION C (3,I) = SURFACE ELEVATION C IN ORIGINAL UNITS C (4,I) = RECIPROCAL RECEIVER NUMBER C AS INDEX PLUS FRAC. DISTANCE C (5,I) = EQUIVALENT CDPN LOCATION C AS INDEX PLUS FRAC. DISTANCE C ( 2-D ARRAY DIMENSIONED: 5 BY NSHOT ) C C IN RECVID R4 RECEIVER IDENTIFICATION ARRAY C (1,I) = X COORDINATE (IN-LINE) C IN UNITS OF GRID-POINT POSITION C (2,I) = Y COORDINATE (CROSS-LINE) C IN UNITS OF GRID-POINT POSITION C (3,I) = SURFACE ELEVATION C IN ORIGINAL UNITS C (4,I) = RECIPROCAL SHOT LOCATION NUMBER C AS INDEX PLUS FRAC. DISTANCE C (5,I) = EQUIVALENT CDPN LOCATION C AS INDEX PLUS FRAC. DISTANCE C ( 2-D ARRAY DIMENSIONED: 5 BY NRECV ) C C IN CDPNID R4 CDP IDENTIFICATION ARRAY C (1,I) = X COORDINATE (IN-LINE) C IN UNITS OF GRID-POINT POSITION C (2,I) = SURFACE ELEVATION C IN ORIGINAL UNITS C ( 2-D ARRAY DIMENSIONED: 2 BY NCDPN ) C C IN CDPREF I2 CDP-TO-SHOT CROSS-REFERENCE INDICES C (1,SIDE,CDPT,CDPN) = SHOT INDEX C (2,SIDE,CDPT,CDPN) = RECEIVER INDEX C (3,SIDE,CDPT,CDPN) = ORTN (TRACE) C ( 4-D ARRAY DIMENSIONED: 3 BY 2 BY NCDPT BY NCDPN ) C C IN RTIMES R4 REFRACTION TIMES IN SECONDS FOR CURRENT C HORIZON (ORIGINAL PICKS PLUS PHANTOMING) C SIGN INDICATES SOURCE: C POSITIVE FOR ORIGINAL PICK C ZERO FOR MISSING PICK C NEGATIVE FOR PHANTOMED PICK C ABSOLUTE VALUE GIVES PROPER TIME C ( 3-D ARRAY DIMENSIONED: LCTPSP BY NSHOT BY NUMHRZ ) C C IN IOCREF I2 REFERENCE CDPT LOCATIONS FOR INWARD AND C OUTWARD CANCELLATION C (1,SIDE,HORIZ,CDPN) = INWARD REFERENCE C FIRST BREAKS ONLY C (DEFAULTS TO -999) C (2,SIDE,HORIZ,CDPN) = OUTWARD REFERENCE C FIRST BREAKS ONLY C (DEFAULTS TO -999) C (3,SIDE,HORIZ,CDPN) = OUTWARD REFERENCE C INCLUDING PHANTOMS C (DEFAULTS TO -999) C ( 4-D ARRAY DIMENSIONED: 3 BY 2 BY NUMHRZ BY NCDPN ) C C C UPDATE ZMODEL R4 THE ELEVATION MODEL FOR ALL HORIZONS C (GRID,CURHRZ,1) = ELEVATION ESTIMATE C (GRID,CURHRZ,2) = TANGENT OF DIP C ( 3-D ARRAY DIMENSIONED: NGRID BY 0:NUMHRZ BY 2 ) C C UPDATE MINELV R4 ELEVATION MODEL MINIMA (PER HORIZON) C ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) C UPDATE MAXELV R4 ELEVATION MODEL MAXIMA (PER HORIZON) C ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) C C C UPDATE PMODEL R4 THE SLOWNESS MODEL FOR THE LAYER BELOW C EACH HORIZON C (GRID,CURHRZ,1) = SLOWNESS ESTIMATE C (GRID,CURHRZ,2) = VELOCITY RMS ERROR C ( 3-D ARRAY DIMENSIONED: NGRID BY 0:NUMHRZ BY 2 ) C C UPDATE MINSLW R4 SLOWNESS MODEL MINIMA (PER HORIZON) C ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) C UPDATE MAXSLW R4 SLOWNESS MODEL MAXIMA (PER HORIZON) C ( 1-D ARRAY DIMENSIONED: 0:NUMHRZ ) C C C UPDATE SMODEL R8 SHOT-CONSISTENT OVERBURDEN VALUES C (SHOT,1) = SURFACE SLOWNESS -- INPUT C (SHOT,2) = LEFT RAY TIME C (SHOT,3) = LEFT INLINE CORRECTION C (SHOT,4) = RIGHT RAY TIME C (SHOT,5) = RIGHT INLINE CORRECTION C ( 2-D ARRAY DIMENSIONED: NSHOT BY 5 ) C C UPDATE RMODEL R8 RECV-CONSISTENT OVERBURDEN VALUES C (RECV,1) = SURFACE SLOWNESS -- INPUT C (RECV,2) = LEFT RAY TIME C (RECV,3) = LEFT INLINE CORRECTION C (RECV,4) = RIGHT RAY TIME C (RECV,5) = RIGHT INLINE CORRECTION C ( 2-D ARRAY DIMENSIONED: NRECV BY 5 ) C C WORK CMODEL R8 TEMPORARY MODEL VALUES C (CDPN,1) = SLOWNESS C (CDPN,2) = VELOCITY SCATTER C (CDPN,3) = ELEVATION C (CDPN,4) = TANGENT OF DIP C (CDPN,5) = CRITICAL OFFSET C (CDPN,6) = OVERLYING SLOWNESS C (CDPN,7) = OVERLYING ELEVATION C (CDPN,8) = << SCRATCH ENTRY >> C (CDPN,9) = << SCRATCH ENTRY >> C ( 2-D ARRAY DIMENSIONED: NCDPN BY 9 ) C C C OUT MTIMES R4 MODEL-PREDICTED TIMES IN SECONDS C SIGN INDICATES SOURCE: C POSITIVE FOR MODELLED PICK C ZERO FOR MISSING PICK C NEGATIVE FOR PHANTOMED PICK C ABSOLUTE VALUE GIVES PROPER TIME C ( 3-D ARRAY DIMENSIONED: LCTPSP BY NSHOT BY NUMHRZ ) C C OUT SRESID R4 SHOT-CONSISTENT LEAST-SQUARES RESIDUALS C (SHOT,1) = POSITIVE RAY (L-R) RESIDUAL C (SHOT,2) = NEGATIVE RAY (R-L) RESIDUAL C ( 2-D ARRAY DIMENSIONED: NSHOT BY 2 ) C C OUT RRESID R4 RECV-CONSISTENT OVERBURDEN VALUES C (RECV,1) = POSITIVE RAY (L-R) RESIDUAL C (RECV,2) = NEGATIVE RAY (R-L) RESIDUAL C ( 2-D ARRAY DIMENSIONED: NRECV BY 2 ) C C C IN NWORK I4 LENGTH OF WORK VECTORS C ( -- MAX OF: NCDPN, NSHOT, NRECV -- C C WORK MODTMP R8 TEMPORARY MODEL VALUES C ( 2-D ARRAY DIMENSIONED: NWORK BY 8 ) C WORK LSFOLD I4 FOLD/DIAGNOSTIC FLAGS C ( 1-D ARRAY DIMENSIONED: NWORK ) C C WORK WORKAR R8 WORK AREA FOR SMOOTHING C ( 1-D ARRAY DIMENSIONED: NWORK BY 5 ) C C C UPDATE CPUTIM I4 CUMMULATIVE CPU TIME C UPDATE CPUDEL I4 RELATIVE CPU TIME C C OUT ERR1 I4 INVERSION ERROR COUNT FLAG C 0 FOR NO ERROR C 1 NO HORIZONS INVERTED C C IF( ERR1 .EQ. 0 ) THEN CALL SARGRMR( IPR, KPBUGF, DATUMF, * INVPRM, EXTOPT, DXGRID, NRESID, * NUMHRZ, NSHOT, NRECV, NCDPN, * NGRID, LCTPSP, NCDPT, * SHOTNO, CDPNUM, ORTNSR, * SHOTID, RECVID, CDPNID, * CDPREF, RTIMES, IOCREF, * ZMODEL, MINELV, MAXELV, * PMODEL, MINSLW, MAXSLW, * SMODEL, RMODEL, CMODEL, * MTIMES, SRESID, RRESID, * NWORK, MODTMP, LSFOLD, WORKAR, * CPUTIM, CPUDEL, ERR1 ) ENDIF C C*********************************************************************** C*** **** C*** VERTICAL RAY TIMES **** C*** **** C*********************************************************************** C C REDEFINE PSEUDO-DATUM IF DEFAULTED C IF( ERR1 .EQ. 0 ) THEN IF( PSEUDO .EQ. UNDEFR ) PSEUDO = MINELV(NUMHRZ) C====================================================================== C C INITIALIZE STATIC REFERENCE ARRAYS C DO 610 SHOT = 1, NSHOT DO 610 ORTN = 1, LCTPSP STATIC(ORTN,SHOT) = 1.0E+10 610 CONTINUE C DO 615 SHOT = 1, NSHOT SMODEL(SHOT,1) = 0.0 615 CONTINUE C DO 620 RECV = 1, NRECV RMODEL(RECV,1) = 0.0 620 CONTINUE C DO 625 CDPN = 1, NCDPN CMODEL(CDPN,1) = 0.0 CMODEL(CDPN,2) = 0.0 CMODEL(CDPN,3) = 0.0 CMODEL(CDPN,4) = 0.0 625 CONTINUE C---------------------------------------------------------------------- C C CALCULATE VERTICAL RAY TIMES C MINVRT = 1.0E+10 MAXVRT = 0.0 NUMLAY = NUMHRZ + 1 C GRIDL = 1 XL = 1.0 C DO 675 GRIDR = 2, NGRID XR = XL + 1.0 C Z2L = PSEUDO Z2R = PSEUDO CURHRZ = NUMHRZ C DO 650 REFHRZ = 1, NUMLAY Z1L = ZMODEL(GRIDL,CURHRZ,1) Z1R = ZMODEL(GRIDR,CURHRZ,1) C PNL = PMODEL(GRIDL,CURHRZ,1) PNR = PMODEL(GRIDR,CURHRZ,1) C DTL = ( Z1L - Z2L )*PNL DTR = ( Z1R - Z2R )*PNR C ------------------------------------------- C C CDP VERTICAL RAY TIME C DO 630 CDPN = 1, NCDPN DX1 = CDPNID(1,CDPN) - XL DX2 = XR - CDPNID(1,CDPN) C IF( 0.0 .LE. DX1 .AND. DX1 .LT. 1.0 ) THEN DT = DTL*DX2 + DTR*DX1 C CMODEL(CDPN,1) = CMODEL(CDPN,1) + DT C IF( CMODEL(CDPN,1) .LT. MINVRT ) * MINVRT = CMODEL(CDPN,1) C IF( CMODEL(CDPN,1) .GT. MAXVRT ) * MAXVRT = CMODEL(CDPN,1) ENDIF 630 CONTINUE C ------------------------------------------- C C SHOTPOINT VERTICAL RAY TIME C DO 635 SHOT = 1, NSHOT DX1 = SHOTID(1,SHOT) - XL DX2 = XR - SHOTID(1,SHOT) C IF( 0.0 .LE. DX1 .AND. DX1 .LT. 1.0 ) THEN DT = DTL*DX2 + DTR*DX1 C SMODEL(SHOT,1) = SMODEL(SHOT,1) + DT ENDIF 635 CONTINUE C ------------------------------------------- C C RECEIVER VERTICAL RAY TIME C DO 640 RECV = 1, NRECV DX1 = RECVID(1,RECV) - XL DX2 = XR - RECVID(1,RECV) C IF( 0.0 .LE. DX1 .AND. DX1 .LT. 1.0 ) THEN DT = DTL*DX2 + DTR*DX1 C RMODEL(RECV,1) = RMODEL(RECV,1) + DT ENDIF 640 CONTINUE C ------------------------------------------- C C UPDATE TERMS FOR NEXT LAYER C CURHRZ = CURHRZ - 1 Z2L = Z1R Z2R = Z1R 650 CONTINUE C ------------------------------------------- C C "FLAT" (CDP) VERTICAL RAY TIME C DTL = ( SURFAV - Z1L )*PNL DTR = ( SURFAV - Z1R )*PNR C DO 660 CDPN = 1, NCDPN DX1 = CDPNID(1,CDPN) - XL DX2 = XR - CDPNID(1,CDPN) C IF( 0.0 .LE. DX1 .AND. DX1 .LT. 1.0 ) THEN CMODEL(CDPN,4) = CMODEL(CDPN,1) * + ( DTL*DX2 + DTR*DX1 ) ENDIF 660 CONTINUE C GRIDL = GRIDR XL = XR 675 CONTINUE C C*********************************************************************** C*** **** C*** STATIC SHIFT ESTIMATION **** C*** **** C*********************************************************************** C C CDP LIMITS FOR REPLACEMENT VELOCITY C BEGCDP = 1 ENDCDP = NCDPN C IF( CDPBEG .LE. CDPEND ) THEN DO 705 CDPI = 1, NCDPN CDPN = CDPNUM(CDPN) C IF( CDPBEG .EQ. CDPN ) BEGCDP = CDPI IF( CDPEND .EQ. CDPN ) ENDCDP = CDPI 705 CONTINUE ENDIF C====================================================================== C C FCDP OR MINS DATUM C IF( DATUMF .EQ. 'FCDP' .OR. DATUMF .EQ. 'MINS' .OR. * DATUMF .EQ. 'REFL' .OR. DATUMF .EQ. 'FLAT' ) THEN DTC = 0.0 COUNT = 0.0 C ------------------------------------------- C C REPLACEMENT VELOCITY C IF( REPVEL .EQ. 0.0 ) THEN NSUM = 0 C DO 725 CDPN = BEGCDP, ENDCDP DO 725 CSIDE = 1, 2 DO 720 CDPT = 1, NCDPT SHOT = CDPREF(1,CSIDE,CDPT,CDPN) RECV = CDPREF(2,CSIDE,CDPT,CDPN) C IF( SHOT .GT. 0 .AND. RECV .GT. 0 ) THEN DTC = DTC + SMODEL(SHOT,1) * + RMODEL(RECV,1) NSUM = NSUM + 1 ENDIF 720 CONTINUE 725 CONTINUE C IF( DTC .NE. 0.0 ) THEN DZ = 2.0*NSUM*( SURFAV - PSEUDO ) REPVEL = DZ/DTC DTC = 0.0 ENDIF ENDIF C ------------------------------------------- C C AVERAGE VERTICAL RAY TIME C DO 775 CDPN = 1, NCDPN DO 750 CSIDE = 1, 2 DO 745 CDPT = 1, NCDPT SHOT = CDPREF(1,CSIDE,CDPT,CDPN) RECV = CDPREF(2,CSIDE,CDPT,CDPN) C IF( SHOT .GT. 0 .AND. RECV .GT. 0 ) THEN DTC = DTC + SMODEL(SHOT,1) * + RMODEL(RECV,1) COUNT = COUNT + 2.0 ENDIF 745 CONTINUE 750 CONTINUE C ------------------------------------------- C C FLOATING CDP DATUM C IF( DATUMF .EQ. 'FCDP' .AND. COUNT .GT. 0.0 ) THEN DT = DTC/COUNT DZ = ( PSEUDO + REPVEL*DT ) COUNT = 0.0 DTC = 0.0 C CMODEL(CDPN,2) = DZ CMODEL(CDPN,3) = DT ENDIF 775 CONTINUE C ------------------------------------------- C C MINIMUM LINE STATIC (MINS) DATUM C CCC IF( DATUMF .EQ. 'MINS' .AND. COUNT .GT. 0.0 ) THEN IF( DATUMF .NE. 'FCDP' .AND. COUNT .GT. 0.0 ) THEN DT = DTC/COUNT DZ = ( PSEUDO + REPVEL*DT ) C DO 790 CDPN = 1, NCDPN CMODEL(CDPN,2) = DZ CMODEL(CDPN,3) = DT 790 CONTINUE ENDIF C====================================================================== C C REFERENCE REFLECTOR FOCUSING C IF( DATUMF .EQ. 'REFL' .OR. DATUMF .EQ. 'FLAT' ) THEN T0AVG = 0.0 C ------------------------------------------- C C REDEFINE STATIC TO FOCUS FIRST HORIZON C DO 825 CDPN = 1, NCDPN X2 = CDPNID(1,CDPN) C GRIDL = X2 GRIDR = GRIDL + 1 C PRVHRZ = NUMHRZ DT = 0.0 DZ = 0.0 C IF( GRIDL .EQ. NGRID ) THEN ZB = ZMODEL(NGRID,NUMHRZ,1) C DO 810 HRZ = 1, NUMHRZ CURHRZ = PRVHRZ PRVHRZ = CURHRZ - 1 C Z2 = ZMODEL(NGRID,PRVHRZ,1) SLOW = PMODEL(NGRID,PRVHRZ,1) Z1 = ZMODEL(NGRID,CURHRZ,1) C DZ = DZ + ( Z2 - Z1 ) DT = DT + ( Z2 - Z1 )*SLOW 810 CONTINUE C ELSE DXL = X2 - GRIDL DXR = GRIDR - X2 C ZB = ZMODEL(GRIDL,NUMHRZ,1)*DXR * + ZMODEL(GRIDR,NUMHRZ,1)*DXL C DO 815 HRZ = 1, NUMHRZ CURHRZ = PRVHRZ PRVHRZ = CURHRZ - 1 C Z2 = ZMODEL(GRIDL,PRVHRZ,1)*DXR * + ZMODEL(GRIDR,PRVHRZ,1)*DXL C SLOW = PMODEL(GRIDL,PRVHRZ,1)*DXR * + PMODEL(GRIDR,PRVHRZ,1)*DXL C Z1 = ZMODEL(GRIDL,CURHRZ,1)*DXR * + ZMODEL(GRIDR,CURHRZ,1)*DXL C DZ = DZ + ( Z2 - Z1 ) DT = DT + ( Z2 - Z1 )*SLOW 815 CONTINUE ENDIF C T0AVG = T0AVG + 2.0*DT/FLOAT(NCDPN) DZ = 2.0*DZ C DO 825 CSIDE = 1, 2 DO 825 CDPT = 1, NCDPT SHOT = CDPREF(1,CSIDE,CDPT,CDPN) RECV = CDPREF(2,CSIDE,CDPT,CDPN) ORTN = CDPREF(3,CSIDE,CDPT,CDPN) C IF( SHOT .GT. 0 .AND. RECV .GT. 0 * .AND. ORTN .GT. 0 ) THEN X1 = SHOTID(1,SHOT) Z1 = SHOTID(3,SHOT) C X3 = RECVID(1,RECV) Z3 = RECVID(3,RECV) C ------------------------------------------- C C CALCULATE MOVEOUT TERM C DX = SQRT( ( X3 - X1 )**2 + DZ**2 ) DTSUM = DX*SLOW C ------------------------------------------- C C CALCULATE STATIC TERM C DX = SQRT( ( X2 - X1 )**2 + ( Z1 - ZB )**2 ) * + SQRT( ( X3 - X2 )**2 + ( Z3 - ZB )**2 ) DT = DX*SLOW C ------------------------------------------- C C MODIFY STATIC FOR FOCUSING C STATIC(ORTN,SHOT) = DTSUM - DT ENDIF 825 CONTINUE C---------------------------------------------------------------------- C C ADD IN RESIDUALS C CCC DX = X3 - X1 C CCC IF( DX .GT. 0.0 ) THEN CCC STATIC(ORTN,SHOT) = STATIC(ORTN,SHOT) CCC * - SL(SHOT) - RR(RECV) CCC ELSE CCC STATIC(ORTN,SHOT) = STATIC(ORTN,SHOT) CCC * - SR(SHOT) - RL(RECV) CCC ENDIF C---------------------------------------------------------------------- C C REDEFINE STATIC TO FLATTEN FIRST HORIZON C IF( DATUMF .EQ. 'FLAT' ) THEN DO 840 CDPN = 1, NCDPN X2 = CDPNID(1,CDPN) C GRIDL = X2 GRIDR = GRIDL + 1 C PRVHRZ = 0 DT = 0.0 - T0AVG C IF( GRIDL .EQ. NGRID ) THEN DO 830 CURHRZ = 1, NUMHRZ Z2 = ZMODEL(NGRID,PRVHRZ,1) SLOW = PMODEL(NGRID,PRVHRZ,1) Z1 = ZMODEL(NGRID,CURHRZ,1) C DT = DT + 2.0*( Z2 - Z1 )*SLOW C PRVHRZ = CURHRZ 830 CONTINUE C ELSE DXL = X2 - GRIDL DXR = GRIDR - X2 C DO 835 CURHRZ = 1, NUMHRZ Z2 = ZMODEL(GRIDL,PRVHRZ,1)*DXR * + ZMODEL(GRIDR,PRVHRZ,1)*DXL C SLOW = PMODEL(GRIDL,PRVHRZ,1)*DXR * + PMODEL(GRIDR,PRVHRZ,1)*DXL C Z1 = ZMODEL(GRIDL,CURHRZ,1)*DXR * + ZMODEL(GRIDR,CURHRZ,1)*DXL C DT = DT + 2.0*( Z2 - Z1 )*SLOW C PRVHRZ = CURHRZ 835 CONTINUE ENDIF C DO 840 CSIDE = 1, 2 DO 840 CDPT = 1, NCDPT SHOT = CDPREF(1,CSIDE,CDPT,CDPN) RECV = CDPREF(2,CSIDE,CDPT,CDPN) ORTN = CDPREF(3,CSIDE,CDPT,CDPN) C IF( SHOT .GT. 0 .AND. RECV .GT. 0 * .AND. ORTN .GT. 0 ) THEN STATIC(ORTN,SHOT) = STATIC(ORTN,SHOT) - DT ENDIF 840 CONTINUE ENDIF ENDIF C====================================================================== C C WAVE-EQUATION (MIGRATION) DATUM C ELSE IF( DATUMF .EQ. 'MIGD' ) THEN C ------------------------------------------- C C REPLACEMENT VELOCITY C IF( REPVEL .EQ. 0.0 ) THEN DTC = 0.0 DO 845 CDPN = BEGCDP, ENDCDP DTC = DTC + CMODEL(CDPN,1) 845 CONTINUE C IF( DTC .NE. 0.0 ) THEN DZ = NCDPN*( SURFAV - PSEUDO ) REPVEL = DZ/DTC ENDIF ENDIF C ------------------------------------------- C C DATUM AND REFERENCE TIME C DO 850 CDPN = 1, NCDPN DT = CMODEL(CDPN,1) DZ = PSEUDO + REPVEL*DT C CMODEL(CDPN,2) = DZ CMODEL(CDPN,3) = DT 850 CONTINUE C====================================================================== C C USER-DEFINED DATUM C ELSE IF( DATUMF .EQ. 'USER' ) THEN IF( NDATUM .GT. 0 ) THEN PICKA = 1 C I = DLOCAT(PICKA) IF( I .GT. 0 ) THEN XA = SHOTID(1,I) ELSE J = 0 - I XA = RECVID(1,J) ENDIF C DO 870 PICKB = 2, NDATUM I = DLOCAT(PICKB) IF( I .GT. 0 ) THEN XB = SHOTID(1,I) ELSE J = 0 - I XB = RECVID(1,J) ENDIF C DX = ( XB - XA ) IF( ABS( DX ) .GT. 1.0E-10 ) THEN DZ = UDATUM(PICKB) - UDATUM(PICKA) SLOPE = DZ/DX C Z1 = UDATUM(PICKA) Z0 = Z1 - SLOPE*XA C DO 860 CDPN = 1, NCDPN X = CDPNID(1,CDPN) C IF( X .LE. XA ) THEN IF( PICKA .EQ. 1 ) THEN CMODEL(CDPN,2) = Z0 + SLOPE*X ENDIF C ELSE IF( X .EQ. XA ) THEN CMODEL(CDPN,2) = UDATUM(PICKA) C ELSE IF( X .GE. XB ) THEN IF( X .EQ. XB ) THEN CMODEL(CDPN,2) = UDATUM(PICKB) C ELSE IF( PICKB .EQ. NDATUM ) THEN CMODEL(CDPN,2) = Z0 + SLOPE*X ENDIF C ELSE CMODEL(CDPN,2) = Z0 + SLOPE*X ENDIF 860 CONTINUE ENDIF C PICKA = PICKB XA = XB 870 CONTINUE C ------------------------------------------- C C REFERENCE TIME FOR STATIC SHIFT ESTIMATES C DO 875 CDPN = 1, NCDPN DZ = 2.0*( CMODEL(CDPN,2) - PSEUDO ) DT = DZ/REPVEL C CMODEL(CDPN,3) = DT 875 CONTINUE ELSE ERR1 = 2 RETURN ENDIF ENDIF C ERR1 = 0 C C*********************************************************************** C*** **** C*** PRINTOUT **** C*** **** C*********************************************************************** C C OUTPUT CPU TIME C CALL CLOCK( CPUDEL, CPUTIM ) WRITE( IPR, 8400 ) CPUTIM C====================================================================== C C LIST MODEL VALUES C MINVEL = +1.0E+10 MAXVEL = -1.0E+10 C DO 950 L = 1, 2 IF( L .EQ. 1 ) WRITE( IPR, 8410 ) IF( L .EQ. 2 ) WRITE( IPR, 8420 ) C DO 905 I = 1,132 OUTPUT(I:I) = ' ' 905 CONTINUE C ------------------------------------------- C OUTPUT(4:7) = 'GRID' C I = 18 J = I + 6 OUTPUT(I:J) = 'SURFACE' I = J + 3 C NXTHRZ = 2 DO 920 CURHRZ = 1, NUMHRZ J = I + 8 OUTPUT(I:J) = 'REFRACTOR' I = J + 2 OUTPUT(I:I) = HRZNUM(CURHRZ:CURHRZ) I = I + 1 OUTPUT(I:I) = '(' I = I + 1 OUTPUT(I:I) = FBNLAB(CURHRZ:CURHRZ) I = I + 1 OUTPUT(I:I) = FBNLAB(NXTHRZ:NXTHRZ) I = I + 1 OUTPUT(I:I) = ')' I = I + 3 C NXTHRZ = NXTHRZ + 1 920 CONTINUE C WRITE( IPR, 8450 ) OUTPUT C ------------------------------------------- C OUTPUT(4:7) = '----' C I = 18 J = I + 6 OUTPUT(I:J) = '-------' I = J + 3 C DO 925 CURHRZ = 1, NUMHRZ J = I + 15 OUTPUT(I:J) = '---------------' I = I + 17 925 CONTINUE C WRITE( IPR, 8450 ) OUTPUT C ------------------------------------------- C C CLEAR LINE BUFFER C DO 940 GRID = 1, NGRID DO 930 I = 1,132 OUTPUT(I:I) = ' ' 930 CONTINUE C ------------------------------------------- C C OUTPUT VALUES FOR CURRENT GRID C I = 3 J = I + 4 WRITE( OUTPUT(I:J), 8425 ) GRID C I = 12 DO 935 CURHRZ = 0, NUMHRZ J = I + 12 C IF( L .EQ. 1 ) THEN WRITE( OUTPUT(I:J), 8430 ) ZMODEL(GRID,CURHRZ,1) ELSE V = 1.0/PMODEL(GRID,CURHRZ,1) WRITE( OUTPUT(I:J), 8430 ) V C IF( MINVEL .GT. V ) MINVEL = V IF( MAXVEL .LT. V ) MAXVEL = V ENDIF C I = I + 17 935 CONTINUE C WRITE( IPR, 8450 ) OUTPUT 940 CONTINUE C IF( L .EQ. 1 ) WRITE( IPR, 8490 ) PSEUDO IF( L .EQ. 2 ) WRITE( IPR, 8495 ) MINVEL, MAXVEL 950 CONTINUE C====================================================================== C C STATIC APPLICATION PARAMETERS C WRITE( IPR, 8600 ) REPVEL IF( DATUMF .EQ. 'FCDP' ) WRITE( IPR, 8610 ) IF( DATUMF .EQ. 'MINS' ) WRITE( IPR, 8620 ) IF( DATUMF .EQ. 'USER' ) WRITE( IPR, 8630 ) IF( DATUMF .EQ. 'MIGD' ) WRITE( IPR, 8640 ) IF( DATUMF .EQ. 'REFL' ) WRITE( IPR, 8641 ) IF( DATUMF .EQ. 'FLAT' ) WRITE( IPR, 8642 ) T0AVG C IF( KPBUGF .GE. 0 ) THEN WRITE( IPR, 8650 ) C DO 960 CDPN = 1, NCDPN WRITE( IPR, 8675 ) CDPN, CMODEL(CDPN,2), CMODEL(CDPN,1), * CMODEL(CDPN,3), CMODEL(CDPN,4) 960 CONTINUE ENDIF C WRITE( IPR, 8695 ) MINVRT, MAXVRT ENDIF RETURN C C*********************************************************************** C*** **** C*** FORMAT STATEMENTS **** C*** **** C*********************************************************************** C 8000 FORMAT(' ',130('*')//56X,I1,'-LAYER INVERSION', * T100,'CPU TIME = ',I12// * ' REFRACTOR: ',8(9X,A1)) 8001 FORMAT(' NUMBER OF PASSES: ',8I10) 8002 FORMAT(' DIAGNOSTICS LEVEL: ',8I10) 8003 FORMAT(' CDP WINDOW WIDTH: ',8I10) 8004 FORMAT(' OFFSET WINDOW WIDTH: ',8I10) 8005 FORMAT(' VELOCITY SMOOTHING LENGTH: ',8I10) 8006 FORMAT(' CRIT. OFFSET SMOOTHING LENGTH: ',8I10) 8007 FORMAT(' DIP ANGLE SMOOTHING LENGTH: ',8I10) 8008 FORMAT(' ELEVATION SMOOTHING LENGTH: ',8I10) C 8011 FORMAT('0 LEFT END EXTRAPOLATION OPTION: ',8(6X,A4)) 8012 FORMAT(' LEFT END REGRESSION LENGTH: ',8I10) C 8021 FORMAT('0 RIGHT END EXTRAPOLATION OPTION: ',8(6X,A4)) 8022 FORMAT(' RIGHT END REGRESSION LENGTH: ',8I10) C 8030 FORMAT('0 ANALYSIS GRID SPACING: ',F12.3/ * ' AVERAGE SURFACE ELEVATION: ',F12.3) C---------------------------------------------------------------------- C C MODEL VALUES C 8400 FORMAT(1X,130('=')//55X,' FINAL MODEL'// * T100,'CPU TIME = ',I12) 8410 FORMAT(55X,'ELEVATION MODEL VALUES:'/) 8420 FORMAT(1X,130('-')//55X,'LAYER VELOCITY MODEL VALUES:'/) C 8425 FORMAT(I5) 8430 FORMAT(F12.3) 8450 FORMAT(A132) C 8490 FORMAT('0PSEUDO-DATUM ELEVATION: ',F12.3) 8495 FORMAT('0TOTAL LAYER VELOCITY RANGE: ',F12.3,' TO ',F12.3) C---------------------------------------------------------------------- C C STATIC APPLICATION PARAMETERS C 8600 FORMAT(1X,130('=')//25X,'STATIC APPLICATION PARAMETERS'// * ' REPLACEMENT VELOCITY = ',F12.3) 8610 FORMAT(' FLOATING CDPN DATUM (FCDP)'/) 8620 FORMAT(' AVERAGED FLOATING CDPN DATUM (MINS)'/) 8630 FORMAT(' USER-DEFINED DATUM'/) 8640 FORMAT(' WAVE-EQUATION (MIGRATION) DATUM'/) 8641 FORMAT(' REFERENCE REFLECTOR FOCUSING'/) 8642 FORMAT(' REFERENCE REFLECTOR FLATTENING AT TIME',F7.3/) C C CDP PARAMETERS C 8650 FORMAT(' DATUM VERTICAL STATIC "FLAT" '/ * ' CDPN ELEVATION RAY TIME REF TIME RAY TIME'/ * ' ---- --------- -------- -------- --------'/) 8675 FORMAT(I5,4F12.3) C C END OF STATICS PARAMETERS C 8695 FORMAT('0TOTAL VERTICAL TIME RANGE: ',F12.3,' TO ',F12.3/ * ' ',130('=')) END