C
C Subroutine file 'apvar.for': Applications and processing of the C results of complete ray tracing -- Part2: Travel-time variations C C Date: 2005, April 24 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: 2005, April 24 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.AND.ISRF1.LE.100) 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