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
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                                                                 ```