C 
C Subroutine file 'apvar.for': Applications and processing of the
C results of complete ray tracing --- Part2: Travel-time variations
C
C Date: 2003, January 18
C Coded by Ludek Klimes
C
C This file consists of the following external procedures:
C     AP29... Subroutine designed to evaluate the variations of the
C             travel time with respect to the model coefficients.
C             It has to be called once at each point along the ray at
C             which the computed quantities are stored, i.e. after each
C             invocation of the subroutine AP00 which reads the
C             quantities into the common block /POINTC/.
C             The subroutine can also optionally calculate the
C             variations of the line integral of the density for the
C             gravimetric inversion.  This option is indicated by NY=8
C             in include file pointc.inc.
C             AP29
C     AP29A...Auxiliary subroutine to AP29.
C             AP29A
C
C=======================================================================
C
C     
C
      SUBROUTINE AP29(NSUM,SUM)
      INTEGER NSUM
      REAL SUM(NSUM)
C
C This subroutine evaluates variations of the travel time with respect
C to the model coefficients.  It has to be called once at each point
C along the ray at which the computed quantities are stored, i.e. after
C each invocation of the subroutine AP00 which reads the quantities into
C the common block /POINTC/.
C The subroutine can also optionally calculate the variations of the
C line integral of the density for the gravimetric inversion.  This
C option is indicated by NY=8 in common block /POINTC/ declared in
C include file pointc.inc.
C Subroutine PARM2 is called to evaluate the material parameters at the
C current point and, at a structural interface, also subroutine SRFC2 is
C called to evaluate the function describing the interface.  After the
C invocation of PARM2 or SRFC2, respectively, subroutine VAR6 is called
C to recall the variations of the model parameters or of the interface,
C with respect to the model coefficients.  If the user replaces the
C subroutine file 'parm.for' or 'srfc.for' by his own version, it is his
C own responsibility to call subroutines VAR1 to VAR5 (see the file
C 'var.for') in such a way that the required variations are stored when
C returning from his own subroutine PARM2 or SRFC2.
C
C Input:
C     NSUM... Total number of the coefficients describing the model.
C     SUM...  Array of dimension at least NSUM, in which the variations
C             of the travel time with respect to the model coefficients
C             are accumulated.  Its elements are set to zeros at the
C             initial point of the ray by this subroutine.
C
C Output:
C     SUM...  Variations of the travel time (from the initial point of
C             the ray to the current point along the ray) with respect
C             to the model coefficients.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C To calculate the variations of the line integral of the density for
C the gravimetric inversion, the folowing variables should be defined:
C     NYF=0
C     NY=8
C     IPT,ICB1I,ICB1,ISRF... The same meaning as for ray tracing.
C     YI(1),Y(1)... Arclength along a straight line.
C     YI(3:5),Y(3:5)... Coordinates, possibly curvilinear.
C     YI(6:8),Y(6:8)... Derivatives of the coordinates with respect to
C             the arclength.
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      INTEGER KOOR
      EXTERNAL KOOR,METRIC,SRFC2,VAR6,AP28
C     SMVPRD
C     PARM2,VELOC
C
C Date: 2003, January 18
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:
C
      INTEGER MFUN
      PARAMETER (MFUN=64)
      INTEGER NFUN1,IFUN1(MFUN),NFUN2,IFUN2(MFUN),ISRF1,II,IBI
      REAL B0I,B1I,B2I,B3I,FUN1(2*MFUN),FUN2(2*MFUN)
      REAL X1,PIN1,PIN2,PIN3,C(3),P1,P2,P3,AUX0
      SAVE NFUN1,IFUN1,FUN1,ISRF1,X1,PIN1,PIN2,PIN3
C
C     NFUN1...Number of functions having nonzero values or nonzero first
C             derivatives at the previous point along the ray.
C     IFUN1...IFUN(1:NFUN1)... Indices in the array SUM corresponding to
C             the functions having nonzero values or nonzero first
C             derivatives at the previous point along the ray.
C     FUN1... FUN(1:NFUN1)... Values of the functions having nonzero
C               values or nonzero first derivatives at the previous
C               point along the ray.
C             FUN(NFUN1+1:2*NFUN1)... First derivatives with respect to
C               the independent variable along the ray at the previous
C               point along the ray, of the functions having nonzero
C               values or nonzero first derivatives.
C             At the first point after the initial point of the ray,
C             the values of NFUN1, IFUN and FUN1 correspond to the
C             initial (zero) point of the ray.
C     NFUN2...Number of functions having nonzero values or nonzero
C             first derivatives at the current point along the ray.
C     IFUN2...Indices in the array SUM corresponding to the functions
C             having nonzero values or nonzero first derivatives at the
C             current point along the ray.
C     FUN2... FUN(1:NFUN2)... Values of the functions having nonzero
C               values or nonzero first derivatives at the current point
C               along the ray.
C             FUN(NFUN2+1:2*NFUN2)... First derivatives with respect to
C               the independent variable along the ray at the previous
C               point along the ray, of the functions having nonzero
C               values or nonzero first derivatives.
C     ISRF1...Index of the surface covering the interface, updated by
C             subroutine AP28.
C     II...   Loop variable (sequential number of the required
C             variation).
C     IBI...  Absolute index of the function coefficient.
C     B0I,B1I,B2I,B3I... Variation of the functional value and the three
C             first derivatives, with respect to the IBI-th coefficient
C             of the model.
C     X1...   Independent variable along the ray, updated by AP28.
C     PIN1,PIN2,PIN3... Contravariant components of the slowness vector
C             at the point of incidence.
C     C...    Coordinates.
C     P1,P2,P3... Contravariant components of the slowness vector.
C     AUX0... Temporary storage location.
C
C.......................................................................
C
C     Initial point of the ray:
      IF(IPT.LE.1) THEN
        PIN1=0.
        PIN2=0.
        PIN3=0.
        CALL AP29A(NY,ICB1I,YI,C,P1,P2,P3,MFUN,NFUN1,IFUN1,FUN1)
      END IF
C
C     Another point of the ray:
      IF(NYF.GT.0) THEN
        CALL AP29A
     *          (NY,ICB1F,YF,C,P1,P2,P3,MFUN,NFUN2,IFUN2,FUN2)
      ELSE
        CALL AP29A
     *          (NY,ICB1 ,Y ,C,P1,P2,P3,MFUN,NFUN2,IFUN2,FUN2)
      END IF
C
C     Numerical quadrature:
      CALL AP28(NSUM,SUM,1,2,0.,X1,ISRF1,
     *                                NFUN1,IFUN1,FUN1,NFUN2,IFUN2,FUN2)
C
C     Structural interface:
      IF(ISRF1.NE.0) THEN
        IF(PIN1.EQ.0..AND.PIN2.EQ.0..AND.PIN3.EQ.0.) THEN
C         incident ray:
          PIN1=P1
          PIN2=P2
          PIN3=P3
        ELSE
C         Reflected/transmitted ray:
C         Including the variation of the travel time with respect to the
C         structural interface
          CALL SRFC2(IABS(ISRF1),C,VD)
          IF(KOOR().NE.0) THEN
            CALL METRIC(C,GSQRD,G,GAMMA)
            AUX0=VD(2)*(G(7)*VD(2)+2.*(G(8)*VD(3)+G(10)*VD(4))) +
     *           VD(3)*(G(9)*VD(3)+2.*G(11)*VD(4)) + VD(4)*G(12)*VD(4)
          ELSE
            AUX0=VD(2)*VD(2)+VD(3)*VD(3)+VD(4)*VD(4)
          END IF
          AUX0=( VD(2)*(P1-PIN1)+VD(3)*(P2-PIN2)+VD(4)*(P3-PIN3) )/AUX0
          II=0
   30     CONTINUE
            II=II+1
            CALL VAR6(1,II,NFUN2,IBI,B0I,B1I,B2I,B3I)
            IF(II.LE.NFUN2) THEN
              IF(IBI.GT.NSUM) THEN
C               729
                CALL ERROR('729 in AP29: Too small input array SUM')
C               Dimension NSUM of input array SUM should be increased.
              END IF
              SUM(IBI)=SUM(IBI)+AUX0*B0I
            END IF
          IF(II.LT.NFUN2) GO TO 30
          PIN1=0.
          PIN2=0.
          PIN3=0.
        END IF
      END IF
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C
C     
C
      SUBROUTINE AP29A(NY,ICB1,Y,C,P1,P2,P3,MFUN,NFUN,IFUN,FUN)
      INTEGER NY,ICB1,MFUN,NFUN,IFUN(MFUN)
      REAL Y(8),C(3),P1,P2,P3,FUN(2*MFUN)
C
C Auxiliary subroutine to AP29.
C
C Input:
C     NY...   If NY=8, line integral of the density for the gravimetric
C             inversion is considered.
C             Otherwise, line integral of slowness for the travel-time
C             inversion is considered.
C     ICB1... Index of the complex block.
C     Y...    Quantities computed along a ray.
C     MFUN... Array dimension.
C
C Output:
C     C...    Coordinates.
C     P1,P2,P3... Contravariant components of the slowness vector.
C     NFUN... Number of variations.
C     IFUN... Indices of variations.
C     FUN...  FUN(1:NFUN)... Values of variations.
C             FUN(NFUN+1:2*NFUN)... First derivatives of variations
C               with respect to the independent variable along the ray.
C
C Subroutines and external functions required:
      INTEGER KOOR
      EXTERNAL KOOR,METRIC,SMVPRD,PARM2,VELOC,VAR6
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
      REAL AUX0,AUX1,AUX2,AUX3,AUX4
      INTEGER NEXPS,IVAL,II
      PARAMETER (NEXPS=0)
      REAL B0I,B1I,B2I,B3I
C
C     AUX0,AUX1,AUX2,AUX3,AUX4... Auxiliary storage locations
C             for local model parameters or temporary variables.
C     IVAL... Index of the function describing the model.
C             IVAL=1 for P-wave,
C             IVAL=2 for S-wave.
C     II...   Loop variable (sequential number of the required
C             variation).
C     B0I,B1I,B2I,B3I... Variation of the functional value and the three
C             first derivatives, with respect to the IBI-th coefficient
C             of the model.
C
C.......................................................................
C
C     Assignments:
      C(1)=Y(3)
      C(2)=Y(4)
      C(3)=Y(5)
C
      IF(NY.EQ.8) THEN
C       Calculating the variations for gravimetric inversion:
        P1=Y(6)
        P2=Y(7)
        P3=Y(8)
        IF(ICB1.EQ.0) THEN
          NFUN=0
          RETURN
        END IF
        CALL PARM2(IABS(ICB1),Y(3),UP,US,AUX0,AUX1,AUX2)
        AUX0=1.
        AUX1=P1
        AUX2=P2
        AUX3=P3
        AUX4=0.
        IVAL=3
      ELSE
C       Calculating the variations for travel-time inversion:
C
C       Contravariant components of the slowness vector:
        IF(KOOR().NE.0) THEN
          CALL METRIC(Y(3),GSQRD,G,GAMMA)
          CALL SMVPRD(G(7),Y(6),Y(7),Y(8),P1,P2,P3)
        ELSE
          P1=Y(6)
          P2=Y(7)
          P3=Y(8)
        END IF
C
C       Material parameters:
        CALL PARM2(IABS(ICB1),Y(3),UP,US,AUX0,AUX1,AUX2)
        CALL VELOC(ICB1,UP,US,AUX1,AUX2,AUX3,AUX4,VD,AUX0)
C       Material parameters and their variations are defined.
C
        AUX0=-VD(1)**(-NEXPS-1)
        AUX4=-VD(1)**(-NEXPS-NEXPS+1)
        AUX1=AUX4*P1
        AUX2=AUX4*P2
        AUX3=AUX4*P3
        AUX4=-FLOAT(NEXPS+1)*(VD(2)*AUX1+VD(3)*AUX2+VD(4)*AUX3)/VD(1)
        IF(ICB1.GT.0) THEN
C         P-wave:
          IVAL=1
        ELSE
C         S-wave:
          IVAL=2
        END IF
      END IF
C
C     Recalling the variations:
      II=0
   20 CONTINUE
        II=II+1
        CALL VAR6(IVAL,II,NFUN,IFUN(II),B0I,B1I,B2I,B3I)
        IF(II.LE.NFUN) THEN
          IF(NFUN.GT.MFUN) THEN
C           730
            CALL ERROR('730 in AP29: Array index out of range')
C           Dimension MFUN of arrays IFUN1, FUN1, IFUN2, FUN2 should
C           be increased.
          END IF
          FUN(II)=AUX0*B0I
          FUN(NFUN+II)=AUX1*B1I+AUX2*B2I+AUX3*B3I+AUX4*B0I
        END IF
      IF(II.LT.NFUN) GO TO 20
      RETURN
      END
C
C=======================================================================
C