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