C
C Subroutine WELL to extract a 1-D vertical profile from the given model C C C SUBROUTINE WELL(COOR,ZSTEP,M,N,XS,IS,XC,IC,D,A,B,R) INTEGER M,N,IS(M),IC(M) REAL COOR(3),ZSTEP,XS(3,M),XC(3,M),D(M),A(M),B(M),R(M) C C Subroutine calculating the material parameters along a given vertical C profile in the model. C C Input: C COOR... Array containing coordinates X1, X2, X3 of a point at C the vertical profile. C ZSTEP...Basic step in the vertical direction. The profile is C composed of line segments. The segments are constructed C from the bottom of the model box towards the Earth C surface. The length of a segment is ZSTEP if the segment C does not intersect an interface. Otherwise, a segment C terminates at the interface, and a new segment extents C from the interface upwards. On output, the segments are C ordered from the Earth surface downwards. The bottommost C segment is the last but one segment, and is followed by C the point (segment of zero length) situated at the bottom C of the model box. C M... Dimension of arrays XS,IS,XC,IC,D,A,B,R. C None of the input parameters are altered. C C Output: C N... Number of segments along the vertical profile. C The material parameters for the segments are calculated at C the central points of the segments. The first point is C situated at the centre of the first segment below the C Earth surface, the last point at the bottom of the model C box. C XS(1:3,I)... Coordinates of the top endpoint of the I-th segment. C IS(I)...Index of the interface if the top endpoints of the I-th C segment is the intersection with the interface, C supplemented by a sign determining whether the segment is C situated at the positive or negative side of the C interface. C Zero otherwise. C XC(1:3,I)... Coordinates of the centre of the I-th segment. C The material parameters for the segment are calculated at C this point. C IC(I)...Index of the complex block in which the I-th segment is C situated. C D(I)... Length of the I-th segment along the vertical profile. C The topmost segment inside each complex block terminates C at the structural interface, the lengths of other C segments are given by ZSTEP. Segment N-1 is situated C directly above the bottom of the model box, D(N)=0. C The points are situated at the centres of the segments. C A(I)... P-wave velocity at the centre of the I-th segment. C B(I)... S-wave velocity at the centre of the I-th segment. C R(I)... Density at the centre of the I-th segment. C Output parameters XS,IS,XC,IC are usually used during inversion, C output parameters D,A,B,R are usually used for forward calculations. C C Common block /MODELC/: INCLUDE 'model.inc' C model.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL NSRFC,SRFC2 INTEGER NSRFC C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for'). C C Subroutines and external functions required: EXTERNAL PARM2,VELOC C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C VELOC... File 'model.for' of the package 'MODEL'. C C Date: 2004, April 5 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C auxmod.inc C C....................................................................... C C Auxiliary storage locations for comparison with program VDISP: CHARACTER*80 FILTST * INTEGER LU1 * PARAMETER (LU1=1) C C Auxiliary storage locations: INTEGER KEND(1),KBOUND(6),IY(8),IX,I REAL X0,X1,X2,XA,XB,ERR,Z REAL Y1(3),D1(3),Y2(3),D2(3),YA(3),DA(3),YB(3),DB(3) DATA KBOUND/1,-1,2,-2,3,-3/ C C....................................................................... C C Test: CALL RSEP3T('SURFTEST',FILTST,' ') IF(FILTST.NE.' ') THEN * OPEN(LU1,FILE=FILTST,STATUS='OLD') * CLOSE(LU1) N=8 A(1)=5.64 B(1)=3.47 R(1)=2.70 D(1)= 6.0 A(2)=6.15 B(2)=3.64 R(2)=2.80 D(2)= 10.5 A(3)=6.60 B(3)=3.85 R(3)=2.85 D(3)= 18.7 A(4)=8.10 B(4)=4.72 R(4)=3.30 D(4)= 80.0 A(5)=8.20 B(5)=4.54 R(5)=3.44 D(5)=100.0 A(6)=8.30 B(6)=4.51 R(6)=3.53 D(6)=100.0 A(7)=8.70 B(7)=4.76 R(7)=3.60 D(7)= 80.0 A(8)=9.30 B(8)=5.12 R(8)=3.76 D(8)= 0.0 RETURN END IF C ERR=0.000100*ZSTEP C X0=0. X1=0. X2=ZSTEP XA=ZSTEP DO 11 I=1,3 D1(I)=0. D2(I)=0. DA(I)=0. Y1(I)=COOR(I) Y2(I)=COOR(I) YA(I)=COOR(I) 11 CONTINUE I=2*IABS(IVERT) IF(IVERT.GT.0) THEN I=I-1 END IF Z=BOUNDM(I) I=IABS(IVERT) D1(I)=FLOAT(ISIGN(1,IVERT)) D2(I)=D1(I) DA(I)=D2(I) Y1(I)=Z Y2(I)=Y1(I)+D1(I)*ZSTEP YA(I)=Y2(I) CALL BLOCK(Y1,0,0,IY(6),IY(4),IY(5)) IF(IY(6).NE.0) THEN C WELL-01 CALL ERROR('WELL-01: Nonzero interface at the model bottom') END IF IF(IY(5).EQ.0) THEN C WELL-02 CALL ERROR('WELL-02: Zero complex block at the model bottom') END IF C Coordinates: IX=M XS(1,IX)=Y1(1) XS(2,IX)=Y1(2) XS(3,IX)=Y1(3) IS(IX)=0 XC(1,IX)=Y1(1) XC(2,IX)=Y1(2) XC(3,IX)=Y1(3) IC(IX)=IY(5) C Material parameters: D(IX)=0. CALL PARM2(IABS(IY(5)),Y1,UP,US,R(IX),QP,QS) CALL VELOC(IY(5),UP,US,QP,QS,A(IX),B(IX),VD,QL) C C Loop over segments: 20 CONTINUE IX=IX-1 IF(IX.LE.0) THEN C WELL-03 CALL ERROR('WELL-03: Too many segments along the profile') END IF IY(6)=0 IY(7)=0 IY(8)=0 CALL CDE(0,0,KEND,6,KBOUND,BOUNDM, * 1,3,3,IY,ERR,X0,X1,Y1,D1,X2,Y2,D2,XA,YA,DA,XB,YB,DB) IF(IY(6).GT.100.AND.IX.EQ.M-1) THEN C WELL-04 CALL ERROR('WELL-04: Crossing the bottom model boundary') END IF IF(IY(6).GT.100) THEN C WELL-05 CALL ERROR('WELL-05: Free space not found') END IF C Centre and length of the segment: Z=0.5*(Y1(I)+YA(I)) D(IX)=XB-X1 C Initial values for the next segment: X1=XB X2=X1+ZSTEP XA=X2 Y1(I)=YB(I) Y2(I)=Y1(I)+D1(I)*ZSTEP YA(I)=Y2(I) C Coordinates: XS(1,IX)=YB(1) XS(2,IX)=YB(2) XS(3,IX)=YB(3) IS(IX)=IY(6) YB(I)=Z XC(1,IX)=YB(1) XC(2,IX)=YB(2) XC(3,IX)=YB(3) IC(IX)=IY(5) C Material parameters: CALL PARM2(IABS(IY(5)),YB,UP,US,R(IX),QP,QS) CALL VELOC(IY(5),UP,US,QP,QS,A(IX),B(IX),VD,QL) C Initial values for the next segment: IF(IY(6).NE.0) THEN IY(4)=IY(7) IY(5)=IY(8) END IF IF(IY(5).NE.0) GO TO 20 C IX=IX-1 N=M-IX DO 80 I=1,N XS(1,I)=XS(1,I+IX) XS(2,I)=XS(2,I+IX) XS(3,I)=XS(3,I+IX) IS(I)=IS(I+IX) XC(1,I)=XC(1,I+IX) XC(2,I)=XC(2,I+IX) XC(3,I)=XC(3,I+IX) IC(I)=IC(I+IX) D(I)=D(I+IX) A(I)=A(I+IX) B(I)=B(I+IX) R(I)=R(I+IX) 80 CONTINUE C RETURN END C C======================================================================= C