C SUBROUTINE FILE 'APVAR.FOR': APPLICATIONS AND PROCESSING OF THE
C RESULTS OF COMPLETE RAY TRACING --- PART2: TRAVEL-TIME VARIATIONS
C
C BY LUDEK KLIMES
C
C THIS FILE CONSISTS OF THE FOLLOWING EXTERNAL PROCEDURES:
C     AP28... SUBROUTINE DESIGNED TO PERFORM THE NUMERICAL QUADRATURE OF
C             THE SET OF GIVEN FUNCTIONS ALONG A RAY.  IT HAS TO BE
C             CALLED ONCE AT EACH POINT ALONG THE RAY IN WHICH THE
C             COMPUTED QUANTITIES ARE STORED, I.E. AFTER EACH INVOCATION
C             OF THE SUBROUTINE AP00 WHICH READS THE QUANTITIES INTO THE
C             COMMON BLOCK /POINTC/.
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     AP29A...AUXILIARY SUBROUTINE TO AP29.
C
C DATE: 1994, JANUARY 15
C CODED BY LUDEK KLIMES
C
C=======================================================================
C
      SUBROUTINE AP28(NSUM,SUM,IX,NDER,STEP,
     *                                NFUN1,IFUN1,FUN1,NFUN2,IFUN2,FUN2)
      INTEGER NSUM,IX,NDER,NFUN1,IFUN1(*),NFUN2,IFUN2(NFUN2)
      REAL SUM(NSUM),STEP,FUN1(*),FUN2(NFUN2*NDER)
C
C THIS SUBROUTINE PERFORMS THE NUMERICAL QUADRATURE OF THE SET OF GIVEN
C FUNCTIONS ALONG A WHOLE RAY.  IT HAS TO BE CALLED ONCE AT EACH POINT
C ALONG THE RAY IN 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
C INPUT:
C     NSUM... TOTAL NUMBER OF THE FUNCTIONS TO BE NUMERICALLY INTEGRATED
C             ALONG THE RAY.
C     SUM...  ARRAY OF DIMENSION AT LEAST NSUM, IN WHICH THE INTEGRALS
C             OF THE GIVEN FUNCTIONS ARE ACCUMULATED.  ITS ELEMENTS ARE
C             SET TO ZEROS AT THE INITIAL POINT OF THE RAY BY THIS
C             SUBROUTINE.
C     IX...   SPECIFIES THE INDEPENDENT VARIABLE ALONG THE RAY:
C             IX=0 INDEPENDENT VARIABLE IS X, I.E. THE SAME AS FOR THE
C               RAY TRACING.
C             IX=1 INDEPENDENT VARIABLE IS THE TRAVEL TIME.
C     NDER... NDER=1 IF JUST THE FUNCTIONAL VALUES OF THE INTEGRATED
C               FUNCTIONS ARE SUBMITTED.  THEN THE RELATIVE ERROR OF THE
C               NUMERICAL QUADRATURE IS PROPORTIONAL TO THE THIRD POWER
C               OF THE STEP ALONG THE RAY (SEE THE PARAMETER STORE IN
C               THE INPUT DATA (3) FOR THE FILE 'RAY.FOR').
C               WHEN INTEGRATING A B-SPLINE IN A REGULAR GRID, THE ERROR
C               IS ABOUT 0.01 FOR THE STEP OF HALF THE SIZE OF THE GRID
C               INTERVAL.
C             NDER=2 IF BOTH THE FUNCTIONAL VALUES AND FIRST DERIVATIVES
C               OF THE INTEGRATED FUNCTIONS ARE SUBMITTED.  THEN THE
C               RELATIVE ERROR OF THE NUMERICAL QUADRATURE IS
C               PROPORTIONAL TO THE FOURTH POWER OF THE STEP ALONG THE
C               RAY (SEE THE PARAMETER STORE IN THE INPUT DATA (3) FOR
C               THE FILE 'RAY.FOR').
C               WHEN INTEGRATING A B-SPLINE IN A REGULAR GRID, THE ERROR
C               IS ABOUT 0.01 FOR THE STEP OF THE SIZE OF THE GRID
C               INTERVAL.
C     STEP... STEP IN THE INDEPENDENT VARIABLE ALONG THE RAY (SEE THE
C             PARAMETER STORE IN THE INPUT DATA (3) FOR THE FILE
C             'RAY.FOR').  REQUIRED JUST IF NDER=1.  IF NDER=1 AND STEP
C             HAS NOT THE CORRECT VALUE, THE RELATIVE ERROR OF THE
C             NUMERICAL QUADRATURE IS PROPORTIONAL TO THE SECOND POWER
C             OF THE ACTUAL STEP ALONG THE RAY.  WHEN INTEGRATING A
C             B-SPLINE WITH NDER=1 AND STEP=0, IN A REGULAR GRID, THE
C             ERROR IS ABOUT 0.01 FOR THE STEP OF THE SIZE OF 0.4 GRID
C             INTERVAL.
C     NFUN1...NUMBER OF FUNCTIONS HAVING NONZERO VALUES (OR NONZERO
C             FIRST DERIVATIVES IF NDER=2) AT THE PREVIOUS POINT ALONG
C             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 IF NDER=2) AT THE PREVIOUS POINT ALONG THE
C             RAY.
C     FUN1... FUN(1:NFUN1)... VALUES OF THE FUNCTIONS HAVING NONZERO
C               VALUES (OR NONZERO FIRST DERIVATIVES IF NDER=2) AT THE
C               PREVIOUS POINT ALONG THE RAY.
C             FUN(NFUN1+1:2*NFUN1)... FOR NDER=2, FIRST DERIVATIVES WITH
C               RESPECT TO THE INDEPENDENT VARIABLE ALONG THE RAY AT THE
C               PREVIOUS POINT ALONG THE RAY, OF THE FUNCTIONS HAVING
C               NONZERO VALUES OR NONZERO FIRST DERIVATIVES.
C             IF THIS SUBROUTINE IS INVOKED AT THE FIRST POINT AFTER THE
C             INITIAL POINT OF THE RAY, THE INPUT VALUES OF NFUN1, IFUN
C             AND FUN1 CORRESPOND TO THE INITIAL (ZERO) POINT OF THE
C             RAY.  AT THE SUBSEQUENT POINTS ALONG THE RAY, THE INPUT
C             VALUES OF NFUN1, IFUN AND FUN1 ARE THE OUTPUT FROM THE
C             PREVIOUS INVOCATION OF THIS SUBROUTINE.
C     NFUN2...NUMBER OF FUNCTIONS HAVING NONZERO VALUES (OR NONZERO
C             FIRST DERIVATIVES IF NDER=2) AT THE CURRENT POINT ALONG
C             THE RAY.
C     IFUN2...INDICES IN THE ARRAY SUM CORRESPONDING TO THE FUNCTIONS
C             HAVING NONZERO VALUES (OR NONZERO FIRST DERIVATIVES IF
C             NDER=2) AT THE CURRENT POINT ALONG THE RAY.
C     FUN2... FUN(1:NFUN2)... VALUES OF THE FUNCTIONS HAVING NONZERO
C               VALUES (OR NONZERO FIRST DERIVATIVES IF NDER=2) AT THE
C               CURRENT POINT ALONG THE RAY.
C             FUN(NFUN2+1:2*NFUN2)... FOR NDER=2, FIRST DERIVATIVES WITH
C               RESPECT TO THE INDEPENDENT VARIABLE ALONG THE RAY AT THE
C               PREVIOUS POINT ALONG THE RAY, OF THE FUNCTIONS HAVING
C               NONZERO VALUES OR NONZERO FIRST DERIVATIVES.
C
C OUTPUT:
C     SUM...  INTEGRALS OF THE GIVEN FUNCTIONS WITH RESPECT TO THE
C             INDEPENDENT VARIABLE ALONG THE RAY, FROM THE INITIAL POINT
C             OF THE RAY TO THE CURRENT POINT ALONG THE RAY (STORED IN
C             THE COMMON BLOCK /POINTC/).
C     NFUN1,IFUN1,FUN1... COPIES OF THE INPUT VALUES OF NFUN2, IFUN2 AND
C             FUN2.
C
C COMMON BLOCK /POINTC/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
C
      REAL X1,X2,W1,W2
      INTEGER ISRF1,ISRF2,ISUM,I
      SAVE X1,ISRF1
C
C     X1...   INDEPENDENT VARIABLE ALONG THE RAY AT THE PREVIOUS POINT.
C     X2...   INDEPENDENT VARIABLE ALONG THE RAY AT THE CURRENT POINT.
C     ISRF1...INDEX OF THE SURFACE AT WHICH THE PREVIOUS POINT IS
C             SITUATED, ZERO INSIDE A COMPLEX BLOCK.
C     ISRF2...INDEX OF THE SURFACE AT WHICH THE CURRENT POINT IS
C             SITUATED, ZERO INSIDE A COMPLEX BLOCK.
C     W1,W2...WEIGTING COEFFICIENTS OF THE NUMERICAL QUADRATURE.
C     ISUM... INDEX IN THE ARRAY SUM CORRESPONDING TO THE FUNCTION UNDER
C             CONSIDERATION.
C     I...    LOOP VARIABLE.
C
C.......................................................................
C
C     INTEGRALS ARE SET TO ZEROS AT THE INITIAL POINT OF THE RAY
      IF(IPT.LE.1) THEN
        DO 10 ISUM=1,NSUM
          SUM(ISUM)=0.
   10   CONTINUE
        ISRF1=1
        IF(IX.LE.0) THEN
          X1=0.
        ELSE
          X1=YI(IX)
        END IF
      END IF
      IF(NYF.GT.0) THEN
        ISRF2=ISRFF
        IF(IX.LE.0) THEN
          X2=XF
        ELSE
          X2=YF(IX)
        END IF
      ELSE
        ISRF2=ISRF
        IF(IX.LE.0) THEN
          X2=X
        ELSE
          X2=Y(IX)
        END IF
      END IF
C
C     NUMERICAL QUADRATURE
      IF(X2.NE.X1) THEN
        W1=(X2-X1)/2.
        W2=W1
        IF(NDER.EQ.1) THEN
          IF(ISRF1.NE.0) THEN
            IF(ISRF2.EQ.0) THEN
C             FIRST INTERVAL OF THE RAY ELEMENT
              W1=W1-(STEP*STEP)/(12.*(X2-X1))
              W2=W2+(STEP*STEP)/(12.*(X2-X1))
            END IF
          ELSE
            IF(ISRF2.NE.0) THEN
C             LAST INTERVAL OF THE RAY ELEMENT
              W1=W1+(STEP*STEP)/(12.*(X2-X1))
              W2=W2-(STEP*STEP)/(12.*(X2-X1))
            END IF
          END IF
        END IF
        DO 21 I=1,NFUN1
          ISUM=IFUN1(I)
          SUM(ISUM)=SUM(ISUM)+W1*FUN1(I)
   21   CONTINUE
        DO 22 I=1,NFUN2
          ISUM=IFUN2(I)
          SUM(ISUM)=SUM(ISUM)+W2*FUN2(I)
   22   CONTINUE
        IF(NDER.EQ.2) THEN
          W1=((X2-X1)**2)/12.
          W2=-W1
          DO 31 I=1,NFUN1
            ISUM=IFUN1(I)
            SUM(ISUM)=SUM(ISUM)+W1*FUN1(NFUN1+I)
   31     CONTINUE
          DO 32 I=1,NFUN2
            ISUM=IFUN2(I)
            SUM(ISUM)=SUM(ISUM)+W2*FUN2(NFUN2+I)
   32     CONTINUE
        END IF
      END IF
C
C     COPYING NFUN2,IFUN2,FUN2 INTO NFUN1,IFUN1,FUN1
      NFUN1=NFUN2
      DO 91 I=1,NFUN2
        IFUN1(I)=IFUN2(I)
        FUN1(I)=FUN2(I)
   91 CONTINUE
      DO 92 I=NFUN2+1,NDER*NFUN2
        FUN1(I)=FUN2(I)
   92 CONTINUE
C
      X1=X2
      ISRF1=ISRF2
      RETURN
      END
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 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/:
      INTEGER IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,NY,ICB1,ISRF
      INTEGER ICB1I,IEND,ISHEET
      REAL XF,YLF(6),YF(39),X,YL(6),Y(39),YLI(6),YI(25)
      COMMON/POINTC/IWAVE,IRAY,IPT,NYF,ICB1F,ISRFF,XF,YLF,YF,
     *                 NY ,ICB1 ,ISRF ,X ,YL ,Y,ICB1I,IEND,ISHEET,YLI,YI
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 ERROR MESSAGES:
C     729...  ARRAY INDEX OUT OF RANGE:
C             DIMENSION MFUN OF ARRAYS IFUN1, FUN1, IFUN2, FUN2 SHOULD
C             BE INCREASED.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR THE METRIC TENSOR ETC.:
C
      REAL G(12),GAMMA(18),GSQRD
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /APCAUX/ BUT TO SHARE THE MEMORY.
      COMMON/APCAUX/ G,GAMMA,GSQRD
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C
C.......................................................................
C
C     AUXILIARY STORAGE LOCATIONS:
C
      INTEGER ISR,II,IBI,MFUN
      PARAMETER (MFUN=64)
      INTEGER NFUN1,IFUN1(MFUN),NFUN2,IFUN2(MFUN)
      REAL B0I,B1I,B2I,B3I,FUN1(2*MFUN),FUN2(2*MFUN)
      REAL PIN1,PIN2,PIN3,C(3),P1,P2,P3,VD(10),AUX0
      SAVE NFUN1,IFUN1,FUN1,PIN1,PIN2,PIN3
C
C     ISR...  INDEX OF THE SURFACE COVERING THE INTERFACE.
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     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     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     VD...   AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS.
C     AUX0... TEMPORARY STORAGE LOCATION.
C
C.......................................................................

C     INITIAL POINT OF THE RAY:
      IF(IPT.LE.1) THEN
        PIN1=0.
        PIN2=0.
        PIN3=0.
        CALL AP29A(ICB1I,0,YI,ISR,C,P1,P2,P3,MFUN,NFUN1,IFUN1,FUN1,VD)
      END IF
C
C     ANOTHER POINT OF THE RAY:
      IF(NYF.GT.0) THEN
        CALL AP29A
     *          (ICB1F,ISRFF,YF,ISR,C,P1,P2,P3,MFUN,NFUN2,IFUN2,FUN2,VD)
      ELSE
        CALL AP29A(ICB1,ISRF,Y,ISR,C,P1,P2,P3,MFUN,NFUN2,IFUN2,FUN2,VD)
      END IF
C
C     NUMERICAL QUADRATURE:
      CALL AP28(NSUM,SUM,1,2,0.,NFUN1,IFUN1,FUN1,NFUN2,IFUN2,FUN2)
C
C     STRUCTURAL INTERFACE:
      IF(ISR.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(ISR),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
              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
      SUBROUTINE AP29A(ICB1,ISRF,Y,ISR,C,P1,P2,P3,MFUN,NFUN,IFUN,FUN,VD)
      INTEGER ICB1,ISRF,ISR,MFUN,NFUN,IFUN(MFUN)
      REAL Y(8),C(3),P1,P2,P3,FUN(2*MFUN),VD(10)
C
C AUXILIARY SUBROUTINE TO AP29.
C
C INPUT:
C     ICB1... INDEX OF THE COMPLEX BLOCK.
C     ISRF... INDEX OF THE SURFACE COVERING THE INTERFACE.
C     Y...    QUANTITIES COMPUTED ALONG A RAY.
C     MFUN... ARRAY DIMENSION.
C     VD...   AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS.
C
C OUTPUT:
C     ISR...  INDEX OF THE SURFACE COVERING THE INTERFACE.
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 THE METRIC TENSOR ETC.:
C
      REAL G(12),GAMMA(18),GSQRD
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /APCAUX/ BUT TO SHARE THE MEMORY.
      COMMON/APCAUX/ G,GAMMA,GSQRD
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C
C.......................................................................
C
      REAL UP(10),US(10),AUX0,AUX1,AUX2,AUX3,AUX4
      INTEGER NEXPS,IVAL,II
      PARAMETER (NEXPS=0)
      REAL B0I,B1I,B2I,B3I
C
C     UP,US,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:
      ISR=ISRF
      C(1)=Y(3)
      C(2)=Y(4)
      C(3)=Y(5)
      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
      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
            PAUSE 'ERROR 729 IN AP29: ARRAY INDEX OUT OF RANGE'
            STOP
          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
