C
C Program CRTPFA calculating two prevailing-frequency
C coupling-ray-theory S-wave Green tensors along each input reference
C ray. The Green tensors are continuous along the reference ray.
C The input reference ray is then doubled. The first of the doubled
C rays is written with the first Green tensor, the second of the doubled
C rays with the second Green tensor.
C
C Version: 7.0
C Date: 2013, June 13
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University in Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.htm
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'SEP'/
C     'SEP'...Name of the file with input parameters.
C             Description of file SEP.
C             If blank, default values of the corresponding data are
C             considered.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file 'sep.for'.
C Names of input and output files related to ray tracing:
C     CRTOUT='string'...File with the names of the output files of
C             program CRT. Only names of the files 'CRT-R', 'CRT-S'
C             and 'CRT-I' are read. Description of file CRTOUT.
C             If CRTOUT is blank, default names are considered.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C               'CRT-S'=' ' and 'CRT-I'='r01i.out'
C               (usually sufficient for the first elementary wave)
C     CRTNEW='string'... File with the names of the output files
C             to contain the new system of doubled rays.
C             The files 'CRT-R' and 'CRT-I' and optionally 'CRT-S' are
C             generated.
C             The form of file CRTNEW is the same as of file CRTOUT.
C             Default: CRTNEW=' ' which means 'CRT-R'='r01-pf.out',
C               'CRT-S'=' ' and 'CRT-I'='r01i-pf.out'
C     MODEL='string'... Name of the input formatted file with the input
C             data for the anisotropic model to be used for the coupling
C             ray theory.
C             Description of file MODEL.
C             No default, MODEL must be specified.
C
C Prevailing-frequency approximation of the coupling ray theory:
C     SINGLF=positive real... Prevailing frequency for the
C             prevailing-frequency approximation of the coupling ray
C             theory.
C             No default, SINGLF must be specified.
C Auxiliary parameters for anisotropic models, need not be specified:
C     ERRWAN=real... Maximum error in propagator matrix along whole ray.
C             Default: ERRWAN=0.001
C     Parameters for optional second-order perturbations of anisotropic
C     travel times of S waves. The perturbations may be calculated only
C     along the isotropic or anisotropic rays which do not cross any
C     interface.
C     The perturbations can not be calculated for P waves.
C     Furthermore, in more complex models in the vicinity of caustics,
C     the second-order perturbations may be infinitely large and should
C     be avoided.
C             QIRAY=integer ... Optional second-order perturbations
C               of anisotropic travel times of S waves.
C               QIRAY=0: No second-order perturbations (default).
C               QIRAY=1: Second-order perturbations of anisotropic
C                 travel times of S waves are calculated and are taken
C                 into account when calculating the travel-time
C                 corrections.
C               Default: QIRAY=0
C             Parameters used for ray tracing. The values of the
C             parameters must be the same as used for ray tracing:
C               CRTANI=real... Switch to anisotropic ray tracing.
C                 Default: CRTANI=0.
C Optional parameters specifying the form of the real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C
C                                                  
C Input formatted file CRTOUT:
C     For general description of file CRTOUT refer to the file writ.for.
C Description specific to this program:
C     'CRT-R'.. Must be specified.
C     'CRT-S'.. May be specified.
C     'CRT-I'.. Must be specified.
C     'CRT-T'.. Not used.
C     Only the first data line is read in the current version.
C     If 'CRT-R' and 'CRT-I' only are specified, all the rays in file
C     'CRT-R' are doubled.  This option is suitable for example for
C     interpolation within ray tubes formed by the coupling rays.
C     If all the three 'CRT-R', 'CRT-S' and 'CRT-I' are specified,
C     only the rays common to files 'CRT-R' and 'CRT-S' are doubled.
C     This option is suitable for example for calculation of Green
C     function along the coupling rays.
C     Defaults for the first data line:
C             'CRT-R'='r01.out', 'CRT-S'=' ', 'CRT-I'='r01i.out'
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL RSEP1,RSEP3I,RSEP3R,RSEP3T,ERROR,WARN,AP00,AP03,AP05,
     *         AP06,AP21,AP28,APYYF,AN03,APW1,APW2,APW3,APW4,
     *         PARM2,PARM3,MODEL1,VELOC,NSRFC,HIVD2,HDER,
     *         KOOR,WAMAT,WAIJKL,WASUM,WASUM3,WASUM4,WASUM5,WACHRI,
     *         WAPROJ,WAVECR,WAANGL,PFWAN,PFREA1,PFREA2,PFINTP,PFPIGE,
     *         PFSHIF,PFPERT,LENGTH,SRFC2,EIGEN,VAL1,VAL2,CURVB1,
     *         SURFB1,VAL3B1,CURVBD,SURFBD,VAL3BD
      INTEGER  KOOR,NSRFC,LENGTH
      REAL WASUM,WASUM3,WASUM4,WASUM5,WAANGL
C     RSEP1,RSEP3I,RSEP3R,RSEP3T ... sep.for.
C     ERROR,WARN ... File error.for.
C     AP00,AP03,AP05,AP06,AP21,AP28,APYYF ... File ap.for.
C     AN03... File an.for.
C     APW1,APW2,APW3,APW4 ... File apw.for.
C     PARM2,PARM3 ... File parm.for.
C     MODEL1,VELOC,NSRFC ... File model.for.
C     HDER... File hder.for.
C     HIVD2.. File means.for.
C     KOOR... File metric.for.
C     WAMAT,WAIJKL,WASUM,WASUM3,WASUM4,WASUM5,WACHRI,WAPROJ,WAVECR,
C     WAANGL ... File wanpfa.for.
C     PFWAN,PFREA1,PFREA2,PFINTP,PFPIGE,PFSHIF,PFPERT ... This file.
C     LENGTH..File length.for.
C     EIGEN...File eigen.for.
C     SRFC2...File srfc.for.
C     VAL1,VAL2 ... File val.for.
C     CURVB1,SURFB1,VAL3B1,CURVBD,SURFBD,VAL3BD ... File fit.for.
C
C-----------------------------------------------------------------------
C
C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points
C of rays and to store the input SEP parameters:
      INCLUDE 'crtpfa.inc'
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER LU0,LUN
      INTEGER LU1,LU2,LU3,LU1N,LU2N,LU3N
      PARAMETER (LU0=1,LUN=2,LU1=3,LU3=4,LU1N=5,LU3N=6)
C                LU2 can be 0 or 7; LU2N can be 0 or 8
      CHARACTER*80 FILSEP,FILMOD,FILCRT,FILNEW
      CHARACTER*80 FILER,FILES,FILEI,FILOR,FILOS,FILOI
      CHARACTER*80 FILERN,FILESN,FILEIN,FILORN,FILOSN,FILOIN
      INTEGER I,IL
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+CRTPFA: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       CRTPFA-01
        CALL ERROR('CRTPFA-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
      WRITE(*,'(A)') '+CRTPFA: Working ...           '
C
C     Name of the model file
      CALL RSEP3T('MODEL',FILMOD,' ')
C
C     Reading the data for a possibly anisotropic model:
      IF(FILMOD.EQ.' ') THEN
C       CRTPFA-02
        CALL ERROR ('CRTPFA-02: Model not specified')
C       File MODEL must be specified.
      ELSE
        OPEN(LU0,FILE=FILMOD,STATUS='OLD')
        CALL MODEL1(LU0)
        CLOSE(LU0)
C       Checking whether the model is anisotropic
        CALL PARM4(I)
        IF(I.GT.0) THEN
C         Isotropic model
C         CRTPFA-03
          CALL ERROR ('CRTPFA-03: Isotropic model')
C         The model specified by parameter MODEL must be anisotropic.
        END IF
      END IF
      IF(KOOR().NE.0) THEN
C       CRTPFA-04
        CALL ERROR
     *      ('CRTPFA-04: Anisotropic model in curvilinear coordinates')
C       Current version of subroutine file 'wan.for' for weakly
C       anisotropic models is coded only for Cartesian coordinates.
      END IF
      IF(NSRFC().NE.0) THEN
C       CRTPFA-05
        CALL WARN('CRTPFA-05: Anisotropic model with interfaces')
C       This version of the program is not debugged for anisotropic
C       models with structural interfaces.
      END IF
C
C     Prevailing frequency and other parameters for coupling ray theory:
      CALL RSEP3R('SINGLF',SINGLF,0.)
      IF (SINGLF.EQ.0.) THEN
C       CRTPFA-06
        CALL ERROR ('CRTPFA-06: Zero prevailing frequency')
C       Prevailing frequency SINGLF must be nonzero.
      ENDIF
      CALL RSEP3R('ERRWAN',ERRWAN,0.001)
      CALL RSEP3I('QIRAY',KQIRAY,0)
      CALL RSEP3R('CRTANI',CRTANI, 0.)
C
C.......................................................................
C     Loop for the CRT output files (i.e. for the lines of file CRTOUT):
      FILOR=' '
      FILOS=' '
      FILOI=' '
      FILORN=' '
      FILOSN=' '
      FILOIN=' '
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      CALL RSEP3T('CRTNEW',FILNEW,' ')
      DO 70 IL=1,1
C     Recently only first line of CRTOUT is read.
C.......................................................................
C
C     Reading output filenames of the CRT program:
      IF(IL.EQ.1) THEN
        FILER='r01.out'
        FILES=' '
        FILEI='r01i.out'
        FILERN='r01-pf.out'
        FILESN=' '
        FILEIN='r01i-pf.out'
      ELSE
        FILER=' '
        FILES=' '
        FILEI=' '
        FILERN=' '
        FILESN=' '
        FILEIN=' '
      END IF
      LU2=0
      LU2N=0
      IF(FILCRT.NE.' ') THEN
        IF(IL.EQ.1) THEN
          OPEN(LU0,FILE=FILCRT,STATUS='OLD')
        END IF
        READ(LU0,*,END=90) FILER,FILES,FILEI
      END IF
      IF(FILER.EQ.' ') THEN
C       No new elementary wave
        GO TO 90
      END IF
      LU2N=0
      IF(FILNEW.NE.' ') THEN
        IF(IL.EQ.1) THEN
          OPEN(LUN,FILE=FILNEW,STATUS='OLD')
        END IF
        READ(LUN,*,END=4) FILERN,FILESN,FILEIN
      END IF
      IF((FILERN.EQ.' '.OR.FILEIN.EQ.' ').OR.
     *   (FILES.NE.' '.AND.FILESN.EQ.' ')) THEN
C       No filenames given
        GO TO 4
      END IF
C
C     Opening output files of the CRT program for input:
      IF(FILER.NE.FILOR) THEN
        IF(IL.GT.1) THEN
          CLOSE(LU1)
        END IF
        OPEN(LU1,FILE=FILER,FORM='UNFORMATTED',STATUS='OLD',ERR=2)
      END IF
      IF((FILES.NE.' ').AND.(FILES.NE.FILOS)) THEN
        LU2=7
        IF(IL.GT.1) THEN
          CLOSE(LU2)
        END IF
        OPEN(LU2,FILE=FILES,FORM='UNFORMATTED',STATUS='OLD',ERR=5)
      END IF
      IF(FILEI.NE.FILOI) THEN
        IF(IL.GT.1) THEN
          CLOSE(LU3)
        END IF
        OPEN(LU3,FILE=FILEI,FORM='UNFORMATTED',STATUS='OLD',ERR=3)
      END IF
      GO TO 9
C       Error messages:
    2   CONTINUE
C       CRTPFA-07
        CALL ERROR('CRTPFA-07: Cannot open the input file with rays')
C       The default name of the file is 'r01.out'.
    5   CONTINUE
C       CRTPFA-08
        CALL ERROR
     *  ('CRTPFA-08: Cannot open input file with points at surfaces')
C       The default name of the file is 's01.out'.
    3   CONTINUE
C       CRTPFA-09
        CALL ERROR
     *  ('CRTPFA-09: Cannot open the input file with initial points')
C       The default name of the file is 'r01i.out'.
    4   CONTINUE
C       CRTPFA-10
        CALL ERROR
     *  ('CRTPFA-10: No filenames of output files given')
C       File CRTOUT contains the names of input files CRT-R and CRT-I,
C       but file CRTNEW does not contain the names of the corresponding
C       input files.
    9 CONTINUE
C     Opening output files:
      IF(FILERN.NE.FILORN) THEN
        IF(IL.GT.1) THEN
          CLOSE(LU1N)
        END IF
        OPEN(LU1N,FILE=FILERN,FORM='UNFORMATTED')
      END IF
      IF((FILESN.NE.' ').AND.(FILESN.NE.FILOSN)) THEN
        LU2N=8
        IF(IL.GT.1) THEN
          CLOSE(LU2N)
        END IF
        OPEN(LU2N,FILE=FILESN,FORM='UNFORMATTED')
      END IF
      IF(FILEIN.NE.FILOIN) THEN
        IF(IL.GT.1) THEN
          CLOSE(LU3N)
        END IF
        OPEN(LU3N,FILE=FILEIN,FORM='UNFORMATTED')
      END IF
C
C.......................................................................
C
C     For all the rays of the files LU1 and LU2,
C     reading the results of the complete ray tracing,
C     calculating new values and writing files LU1N and LU2N:
      CALL PFWAN(LU1,LU2,LU3,LU1N,LU2N,LU3N)
C
C.......................................................................
C
      FILOR=FILER
      FILOS=FILES
      FILOI=FILEI
      FILORN=FILERN
      FILOSN=FILESN
      FILOIN=FILEIN
   70 CONTINUE
C     End of loop for input/output files.
   90 CONTINUE
C
C.......................................................................
C
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU1N)
      CLOSE(LU2N)
      IF(FILCRT.NE.' ') CLOSE(LU0)
      IF(FILNEW.NE.' ') CLOSE(LUN)
      WRITE(*,'(A)') '+CRTPFA: Done.                 '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE PFWAN(LU1,LU2,LU3,LU1N,LU2N,LU3N)
C
C----------------------------------------------------------------------
C
C Subroutine PFWAN to convert a system of common S-wave rays in files
C LU1,LU2 and LU3 into two systems of coupling-ray theory rays using the
C prevailing frequency approximation and to write them into files
C LU1N,LU2N and LU3N
      INTEGER LU1,LU2,LU3,LU1N,LU2N,LU3N
C     ------------------------------------------------------------------
C Input:
C     LU1 ... Number of the logical unit connected to the CRT output
C             file `CRT-R` with quantities along rays.
C     LU2 ... Number of the logical unit connected to the CRT output
C             file `CRT-S` with quantities at storing surface.
C             LU2=0 if the file is not to be read.
C     LU3 ... Number of the logical unit connected to the CRT output
C             file `CRT-I` with quantities at the initial points of the
C             rays.
C     LU1N... Number of the logical unit connected to the new file
C             `CRT-R` where the quantities along new coupling-ray-theory
C             rays will be written.
C     LU2N... Number of the logical unit connected to the new file
C             `CRT-S` where the quantities of new coupling-ray-theory
C             rays at storing surface will be written.
C             Not used if LU2=0.
C     LU3N... Number of the logical unit connected to the new file
C             `CRT-I` where the quantities at the initial points of the
C             new coupling-ray-theory rays will be written.
C No output
C
C.......................................................................
C Subroutines and external functions directly called:
      EXTERNAL AP00,APYYF,APW1,APW2,APW3,APW4,PFREA1,PFREA2,PFINTP,
     *         PFPIGE,PFSHIF
C     AP00,APYYF ... File ap.for.
C     APW1,APW2,APW3,APW4 ... File apw.for.
C     PFREA1,PFREA2,PFINTP,PFPIGE,PFSHIF ... This file.
C
C.......................................................................
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C ...........................
C
C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points
C of rays and to store the input SEP parameters:
      INCLUDE 'crtpfa.inc'
C ...........................
C
C Common block /RAM/ to store the points of the second coupling ray:
      INCLUDE 'ram.inc'
      INTEGER MISTOK,MSTOK
      PARAMETER (MISTOK=MRAM/2,MSTOK=MRAM/2)
      INTEGER ISTOK(MISTOK)
      REAL STOK(MSTOK)
      EQUIVALENCE (STOK,RAM),(ISTOK,RAM(MSTOK+1))
C
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY0
      REAL PIGET1,PIGET2,PIGR1(2,2),PIGI1(2,2),PIGR2(2,2),PIGI2(2,2)
      REAL PI
      PARAMETER (PI=3.1415926)
      LOGICAL LAPYYF,LWLU2,LSTOR2
C
C-----------------------------------------------------------------------
C
C
      NPTS=0
      ISTOK(1)=0
      LAPYYF=.FALSE.
      LWLU2=.FALSE.
      LSTOR2=.FALSE.
C     Loop for the points of rays
   50 CONTINUE
C       Reading the results of the complete ray tracing
        IF (LU2.EQ.0) THEN
          CALL AP00(0,LU1,LU3)
        ELSE
          CALL AP00(LU1,LU2,LU3)
        ENDIF
        IF((IPT.LE.1).AND.(ISTOK(1).NE.0)) THEN
C         New ray - recording the duplication of the last ray
          CALL APW4(LU1N,ISTOK,STOK)
          LAPYYF=.FALSE.
          LWLU2=.FALSE.
          LSTOR2=.FALSE.
        END IF
        IF(IWAVE.LT.1) THEN
C         End of rays
          GO TO 100
        END IF
        IF (LU2.NE.0) THEN
C         Reading and writing points O/S at storing surfaces:
          IF ((YF(1).LE.Y(1)).OR.LWLU2) THEN
C           The point O/F is before the point O/S, or the point O/S
C           was already written. Using point O/F:
            CALL APYYF
            LAPYYF=.TRUE.
            IF (YF(1).EQ.Y(1)) LSTOR2=.TRUE.
          ELSE
C           Point O/F is after point O/S, point O/S is to be used:
            LSTOR2=.TRUE.
          ENDIF
        ENDIF
        IF(IPT.LE.1)THEN
C         New ray - recording the initial point:
          CALL PFREA1
          IF (Y(1).EQ.YI(1)) THEN
C           The first point on the ray has the same travel time
C           as the point in the source. This may happen for rays
C           with zero length (rays terminating in the source).
C           This point will be skipped:
            GOTO 50
          ENDIF
        ENDIF
  60    CONTINUE
C         Recording the new point:
          CALL PFREA2
C         Calculating the Christoffel matrix and other values in the new
C         point(s), adding points by interpolation (if necessary):
          CALL PFINTP
C         Calculating the propagator matrix and travel-time corrections:
          CALL PFPIGE(PIGET1,PIGET2,PIGR1,PIGI1,PIGR2,PIGI2,DPIGD)
          IF(IPT.LE.1)THEN
C           Calculating the initial points of the first and
C           second coupling ray, recording them to the file CRT-I:
            IRAY0=IRAY
            IRAY=2*IRAY0-1
            PVI(1)=RFEVEC(1,1)
            PVI(2)=RFEVEC(2,1)
            PVI(3)=RFEVEC(3,1)
            PVI(4)=RFEVEC(4,1)
            PVI(5)=RFEVEC(5,1)
            PVI(6)=RFEVEC(6,1)
            CALL APW1(LU3N)
            IRAY=2*IRAY0
            CALL APW1(LU3N)
            IRAY=IRAY0
          ENDIF
          IF ((LU2.EQ.0).OR.LAPYYF) THEN
C           Calculating points O/F on the two coupling rays,
C           recording them to the file CRT-R:
            IRAY0=IRAY
C           Calculating and recording the point of the first coupling ray:
            IRAY=2*IRAY0-1
            PV(1)=RFEVEC(1,NPTS)
            PV(2)=RFEVEC(2,NPTS)
            PV(3)=RFEVEC(3,NPTS)
            PV(4)=RFEVEC(4,NPTS)
            PV(5)=RFEVEC(5,NPTS)
            PV(6)=RFEVEC(6,NPTS)
            Y(1)=PIGET1
            Y(6)=Y(6)*(1.-DPIGD)
            Y(7)=Y(7)*(1.-DPIGD)
            Y(8)=Y(8)*(1.-DPIGD)
            Y(28)=PIGR1(1,1)
            Y(29)=PIGI1(1,1)
            Y(30)=PIGR1(2,1)
            Y(31)=PIGI1(2,1)
            Y(32)=PIGR1(1,2)
            Y(33)=PIGI1(1,2)
            Y(34)=PIGR1(2,2)
            Y(35)=PIGI1(2,2)
            CALL APW2(LU1N)
C           Calculating and recording the point of the second coupling ray:
            IRAY=2*IRAY0
            Y(1)=PIGET2
            Y(6)=Y(6)*(1.+DPIGD)
            Y(7)=Y(7)*(1.+DPIGD)
            Y(8)=Y(8)*(1.+DPIGD)
            Y(28)=PIGR2(1,1)
            Y(29)=PIGI2(1,1)
            Y(30)=PIGR2(2,1)
            Y(31)=PIGI2(2,1)
            Y(32)=PIGR2(1,2)
            Y(33)=PIGI2(1,2)
            Y(34)=PIGR2(2,2)
            Y(35)=PIGI2(2,2)
            CALL APW3(ISTOK,MISTOK,STOK,MSTOK)
            IRAY=IRAY0
            IF (LAPYYF) THEN
              CALL APYYF
              LAPYYF=.FALSE.
            ENDIF
          ENDIF
          IF (LSTOR2) THEN
C           Calculating points O/S on the two coupling rays,
C           recording them to the file CRT-S:
            IF (LAPYYF) THEN
C             CRTPFA-11
              CALL ERROR('CRTPFA-11: Missorder in points along ray')
C             Point O/F should be at or behind point O/S, so APYYF
C             should not have been called and LAPYYF should be .FALSE..
            ENDIF
            IRAY0=IRAY
C           Calculating and recording the point of the first coupling ray:
            IRAY=2*IRAY0-1
            PV(1)=RFEVEC(1,NPTS)
            PV(2)=RFEVEC(2,NPTS)
            PV(3)=RFEVEC(3,NPTS)
            PV(4)=RFEVEC(4,NPTS)
            PV(5)=RFEVEC(5,NPTS)
            PV(6)=RFEVEC(6,NPTS)
            Y(1)=PIGET1
            Y(6)=Y(6)*(1.-DPIGD)
            Y(7)=Y(7)*(1.-DPIGD)
            Y(8)=Y(8)*(1.-DPIGD)
            Y(28)=PIGR1(1,1)
            Y(29)=PIGI1(1,1)
            Y(30)=PIGR1(2,1)
            Y(31)=PIGI1(2,1)
            Y(32)=PIGR1(1,2)
            Y(33)=PIGI1(1,2)
            Y(34)=PIGR1(2,2)
            Y(35)=PIGI1(2,2)
            CALL APW2(LU2N)
C           Calculating and recording the point of the second coupling ray:
            IRAY=2*IRAY0
            Y(1)=PIGET2
            Y(6)=Y(6)*(1.+DPIGD)
            Y(7)=Y(7)*(1.+DPIGD)
            Y(8)=Y(8)*(1.+DPIGD)
            Y(28)=PIGR2(1,1)
            Y(29)=PIGI2(1,1)
            Y(30)=PIGR2(2,1)
            Y(31)=PIGI2(2,1)
            Y(32)=PIGR2(1,2)
            Y(33)=PIGI2(1,2)
            Y(34)=PIGR2(2,2)
            Y(35)=PIGI2(2,2)
            CALL APW2(LU2N)
            IRAY=IRAY0
            LWLU2=.TRUE.
            LSTOR2=.FALSE.
            IF (YF(1).GT.Y(1)) THEN
C             Point O/F is behind point O/S, point O/S was recently
C             processed, now continuing by processing of point O/F:
              CALL APYYF
              LAPYYF=.TRUE.
              CALL PFSHIF(NPTS,1)
              GOTO 60
            ENDIF
          ENDIF
          IF (LAPYYF) THEN
            CALL APYYF
            LAPYYF=.FALSE.
          ENDIF
C         Shifting the values from "new" point to "old" point:
          CALL PFSHIF(NPTS,1)
      GO TO 50
C
  100 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE PFREA1
C
C----------------------------------------------------------------------
C Subroutine to store the unformatted output of program CRT.
C Subroutine assumes that subroutine AP00 of file 'ap.for' was called
C recently and the values are stored in arrays of common block POINTC,
C this subroutine calculates the necessary quantities and stores them
C into arrays of common block PFPTS.
C
C No input, the values are read from common block POINTC.
C No output, calculated quantities are written to common block PFPTS.
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C ...........................
C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points
C of rays and to store the input SEP parameters:
      INCLUDE 'crtpfa.inc'
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER NPTSO,NPTSN,NPTSA,ISTART,ISTOP,ISTEP,ISHIFT,I1,I2,I3
      REAL UP(10),US(10),RHO,QPA,QSA,VPA,VSA,VD(10),QL
      REAL HI(18),H(18),HUI(9)
C Subroutines and external functions directly called:
      EXTERNAL ERROR,AP03,AN03,AP05,AP06,PARM2,VELOC
C     ERROR.. File error.for.
C     AP03,AP05,AP06 ... File ap.for.
C     AN03... File an.for.
C     PARM2.. File parm.for.
C     VELOC.. File model.for.
C
C-----------------------------------------------------------------------
C
C     Recording initial point of a new ray:
      NPTS=1
C     Index of a point:
      IFIND(NPTS)=NPTS
C     Travel time:
      RFTT(NPTS)=YI(1)
C     Storing the travel time in the initial point:
      RFTT(0)=YI(1)
C     Coordinates:
      RFCOOR(1,NPTS)=YI(3)
      RFCOOR(2,NPTS)=YI(4)
      RFCOOR(3,NPTS)=YI(5)
C     Slowness vector:
      RFSLO1(NPTS)=YI(6)
      RFSLO2(NPTS)=YI(7)
      RFSLO3(NPTS)=YI(8)
C     Polarization vector:
      RFPOL1(NPTS)=YI(9)
      RFPOL2(NPTS)=YI(10)
      RFPOL3(NPTS)=YI(11)
C     Index of complex block:
      IFICB(NPTS)=ICB1I
C     S-wave velocity and their derivatives:
      IF (CRTANI.EQ.0) THEN
        RFVS(NPTS)=YLI(2)
        RFVS1(NPTS)=YLI(4)
        RFVS2(NPTS)=YLI(5)
        RFVS3(NPTS)=YLI(6)
      ELSE
        CALL PARM2(IFICB(NPTS),RFCOOR(1,NPTS),UP,US,RHO,QPA,QSA)
        CALL VELOC(IFICB(NPTS),UP,US,QPA,QSA,VPA,VSA,VD,QL)
        RFVS(NPTS)=VD(1)
        RFVS1(NPTS)=VD(2)
        RFVS2(NPTS)=VD(3)
        RFVS3(NPTS)=VD(4)
      ENDIF
C     Matrix of geometrical spreading:
      RFQMAT(1,NPTS)=YI(12)
      RFQMAT(2,NPTS)=YI(13)
      RFQMAT(3,NPTS)=YI(16)
      RFQMAT(4,NPTS)=YI(17)
C     Transformation matrix P:
      RFPMAT(1,NPTS)=YI(14)
      RFPMAT(2,NPTS)=YI(15)
      RFPMAT(3,NPTS)=YI(18)
      RFPMAT(4,NPTS)=YI(19)
C     Index of the surface:
      IFISRF(NPTS)=0
C     Covariant and contravariant components of the basis vectors
C     of the RCCS:
      CALL AP03(0,RFH(1,NPTS),H,HUI)
      IF (CRTANI.NE.0) THEN
        CALL AN03(ICB1I,YI,RFH(1,NPTS))
      ENDIF
      RETURN
C=======================================================================
C
      ENTRY PFREA2
C
C.......................................................................
C
C Entry designed to read the output of the CRT at other than initial
C point of a ray.
C
      NPTS=NPTS+1
      IF (NPTS.GT.MPTS) THEN
C       CRTPFA-12
        CALL ERROR('CRTPFA-12: NPTS greater than MPTS')
C       This error should not happen, NPTS should equal 2 here.
      ENDIF
C     Index of a point:
      IFIND(NPTS)=IFIND(NPTS-1)+1
C     Travel time:
      RFTT(NPTS)=Y(1)
C     Coordinates:
      RFCOOR(1,NPTS)=Y(3)
      RFCOOR(2,NPTS)=Y(4)
      RFCOOR(3,NPTS)=Y(5)
C     Slowness vector:
      RFSLO1(NPTS)=Y(6)
      RFSLO2(NPTS)=Y(7)
      RFSLO3(NPTS)=Y(8)
C     Polarization vector
      RFPOL1(NPTS)=Y(9)
      RFPOL2(NPTS)=Y(10)
      RFPOL3(NPTS)=Y(11)
C     Index of complex block:
      IFICB(NPTS)=ICB1
C     S-wave velocity and their derivatives:
      IF (CRTANI.EQ.0) THEN
        RFVS(NPTS)=YL(2)
        RFVS1(NPTS)=YL(4)
        RFVS2(NPTS)=YL(5)
        RFVS3(NPTS)=YL(6)
      ELSE
        CALL PARM2(IFICB(NPTS),RFCOOR(1,NPTS),UP,US,RHO,QPA,QSA)
        CALL VELOC(IFICB(NPTS),UP,US,QPA,QSA,VPA,VSA,VD,QL)
        RFVS(NPTS)=VD(1)
        RFVS1(NPTS)=VD(2)
        RFVS2(NPTS)=VD(3)
        RFVS3(NPTS)=VD(4)
      ENDIF
C     Matrix of geometrical spreading:
      CALL AP05(0,HUI,
     *  RFQMAT(1,NPTS),RFQMAT(2,NPTS),RFQMAT(3,NPTS),RFQMAT(4,NPTS))
C     Transformation matrix P:
      CALL AP06(0,HUI,
     *  RFPMAT(1,NPTS),RFPMAT(2,NPTS),RFPMAT(3,NPTS),RFPMAT(4,NPTS))
C     Index of the surface:
      IFISRF(NPTS)=ISRF
C     Covariant and contravariant components of the basis vectors
C     of the RCCS:
      CALL AP03(0,HI,RFH(1,NPTS),HUI)
      IF (CRTANI.NE.0) THEN
        CALL AN03(ICB1,Y,RFH(1,NPTS))
      ENDIF
      RETURN
C=======================================================================
C
      ENTRY PFSHIF(NPTSO,NPTSN)
C
C.......................................................................
C
C Entry designed to shift the values in arrays of common block PFPTS.
C The entry shifts all the values stored from NPTSO to NPTS so that
C the values are stored from NPTSN.
C Input:
C     NPTSO...Position in arrays of common block PFPTS from which the
C             values are moved.
C     NPTSN...Position in arrays of common block PFPTS to which the
C             values are moved. All the quantities stored in positions
C             NPTSO to NPTS are moved, and NPTS is then set to
C             NPTS+NPTSN-NPTSO.
C Example 1: NPTSO=3, NPTSN=4, NPTS=10:
C     Values in positions 3 to 10 will be shifted to positions 4 to 11,
C     position 3 will be free, NPTS will be 11 on output.
C Example 2: NPTSO=10, NPTSN=1, NPTS=10:
C     Value in position 10 will be shifted to position 1,
C     NPTS will be 1 on output, i.e. values originaly in positions 1 to
C     9 will be erased on output.
C
      IF (NPTSN.EQ.NPTSO) RETURN
      NPTSA=NPTS+NPTSN-NPTSO
      IF (NPTSA.GT.MPTS) THEN
C       CRTPFA-13
        CALL ERROR('CRTPFA-13: Arrays in common block PFPTS small')
C       The dimension MPTS of the arrays in common block PFPTS
C       is given in the include file crtpfa.inc.
      ENDIF
      IF (NPTSN.GT.NPTSO) THEN
        ISTART=NPTSA
        ISTOP=NPTSO+1
        ISTEP=-1
        ISHIFT=-1
      ELSE
        ISTART=NPTSN
        ISTOP=NPTSA
        ISTEP=1
        ISHIFT=NPTSO-NPTSN
      ENDIF
      DO 30, I1=ISTART,ISTOP,ISTEP
        I2=I1+ISHIFT
        IFIND(I1)=   IFIND(I2)
        IFICB(I1)=   IFICB(I2)
        IFISRF(I1)=  IFISRF(I2)
        RFTT(I1)=    RFTT(I2)
        RFSLO1(I1)=  RFSLO1(I2)
        RFSLO2(I1)=  RFSLO2(I2)
        RFSLO3(I1)=  RFSLO3(I2)
        RFPOL1(I1)=  RFPOL1(I2)
        RFPOL2(I1)=  RFPOL2(I2)
        RFPOL3(I1)=  RFPOL3(I2)
        RFVS(I1)=    RFVS(I2)
        RFVS1(I1)=   RFVS1(I2)
        RFVS2(I1)=   RFVS2(I2)
        RFVS3(I1)=   RFVS3(I2)
        RFGS1(I1)=   RFGS1(I2)
        RFGS2(I1)=   RFGS2(I2)
        RFGP(I1)=    RFGP(I2)
        RFANGL(I1)=  RFANGL(I2)
        RFTINT(I1)=  RFTINT(I2)
        DO 32, I3=1,4
          RFQMAT(I3,I1)=RFQMAT(I3,I2)
          RFPMAT(I3,I1)=RFPMAT(I3,I2)
  32    CONTINUE
        DO 34, I3=1,18
          RFH(I3,I1)=   RFH(I3,I2)
  34    CONTINUE
        DO 35, I3=1,9
          RFEVEC(I3,I1)= RFEVEC(I3,I2)
  35    CONTINUE
        DO 36, I3=1,3
          RFCOOR(I3,I1)=RFCOOR(I3,I2)
  36    CONTINUE
        DO 38, I3=1,2
          RFES1P(I3,I1)=RFES1P(I3,I2)
          RFES2P(I3,I1)=RFES2P(I3,I2)
  38    CONTINUE
  30  CONTINUE
      NPTS=NPTSA
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE PFINTP
C
C-----------------------------------------------------------------------
C Subroutine to check the angular change between two points
C on the ray, to compute all the auxiliary quantities in the points,
C and, if necessary, to add new points on the rays by interpolation.
C Computation of elastic parameters is completed by invocation
C of subroutine PARM3 of file 'parm.for'.
C Then Christoffel matrix is evaluated and its eigenvalues and
C eigenvectors are computed by subroutine EIGEN of file 'eigen.for'.
C In the next step the angular difference DELTFI is computed for
C each subinterval along the ray.
C If the value of DELTFI is greater than prescribed limit, new points
C are added on the ray.  Coordinates, slowness vector and polarization
C vector are interpolated by subroutine HIVD2 of the file 'means.for'.
C Material parameters, Christoffel matrix and their eigenvalues and
C eigenvectors are calculated in the new point. Matrix of geometrical
C spreading, transformation matrix P and covariant and contravariant
C components of the basis vectors of the RCCS are interpolated linearly.
C
C ...........................
C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points
C of rays and to store the input SEP parameters:
      INCLUDE 'crtpfa.inc'
C ...........................
C External functions directly called:
      EXTERNAL ERROR,HIVD2,PARM2,PARM3,VELOC,
     *         PFSHIF,WACHRI,WAPROJ,WAVECR,WAANGL
      REAL WAANGL
C     ERROR.. File error.for.
C     HIVD2.. File means.for.
C     PARM2,PARM3.. File parm.for.
C     VELOC.. File model.for.
C     PFSHIF..This file.
C     WACHRI,WAPROJ,WAVECR,WAANGL ... File wanpfa.for.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER NPTSE,LPTISO
      REAL UP(10),US(10),RHO,QP,QS,VP,VS,VD(10),QL,AA(10,21),QQ(21)
      REAL EE(9),DER(9)
      REAL DEFI,DELTFI,PI
      PARAMETER (PI=3.1415926)
      INTEGER I1
      INTEGER MQANT,MNEWP,NNEWP
      PARAMETER (MNEWP=50)
      PARAMETER (MQANT=56)
      REAL ROLDP(MQANT),RNEWP(MQANT,MNEWP)
      INTEGER KNEWP(MNEWP)
      REAL TTOLD,XOLD(3),DXOLD(3),POLD(3),DPOLD(3),EOLD(3),DEOLD(3)
      REAL TTNEW,XNEW(3),DXNEW(3),PNEW(3),DPNEW(3),ENEW(3),DENEW(3),AUX
      SAVE ROLDP,TTOLD,XOLD,DXOLD,POLD,DPOLD,EOLD,DEOLD
C
C     ROLDP(I),RNEWP(I,J)
C           I=1    ...   Travel time.
C             2-4  ...   Coordinates.
C             5-7  ...   Slowness vector.
C             8-10 ...   Polarization vector E1.
C             11   ...   qP eigenvalue of Christoffel matrix.
C             12   ...   qS1 eigenvalue of Christoffel matrix.
C             13   ...   qS2 eigenvalue of Christoffel matrix.
C             14   ...   First component of the qS1 eigenvector
C                  projected to the plane perpendicular to the ray,
C                  defined by two basis vectors of ray-centered
C                  coordinate system.
C             15   ...   Second component of the qS1 eigenvector.
C             16   ...   First component of the qS2 eigenvector.
C             17   ...   Second component of the qS2 eigenvector.
C             18-26...   Eigenvectors of qS1, qS2 and P wave.
C             27   ...   S-wave velocity.
C             28-30...   Derivatives of S-wave velocity.
C             31-34...   Matrix of geometrical spreading.
C             34-38...   Transformation matrix P.
C             39-56...   Covariant and contravariant components of
C                        the basis vectors of the RCCS.
C     KNEWP...Division index of points in RNEWP.
C     AA..    Values, first and second partial derivatives of real
C             parts of 21 reduced (divided by the density) elastic
C             parameters.  The order of the value, first and second
C             partial derivatives of each parameter Aij is:
C             Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33.
C             The order of parameters (second array index) is:
C             A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45,
C             A55,A16,A26,A36,A46,A56,A66.
C     RHO...  Density at the given point.
C     QQ ...  Imaginary parts of 21 reduced elastic parameters at the
C             given point, ordered as
C             Q11,Q12,Q22,Q13,Q23,Q33,Q14,Q24,Q34,Q44,Q15,Q25,Q35,Q45,
C             Q55,Q16,Q26,Q36,Q46,Q56,Q66.
C     EE ...  Eigenvectors of the Christoffel matrix.
C     DER ... Derivatives dx/dt (anisotropic).
C     DEFI .. Maximum allowed angular change for eigenvectors of
C             Christoffel matrix between neighboring points on the ray.
C     TTOLD...Travel time in the old point stored for interpolation.
C     XOLD ...Coordinates of the old point stored for interpolation.
C     DXOLD ..Derivatives in the old point stored for interpolation.
C     POLD ...Slowness vector in the old point stored for interpolation.
C     DPOLD ..Derivatives of the slowness vector in the old point.
C     EOLD ...Polarization vector in the old point.
C     DEOLD ..Derivatives of the polarization vector in the old point.
C     TTNEW,XNEW,DXNEW,PNEW,DPNEW,ENEW,DENEW ... The same quantities
C             in the new point. Here old point and new point mean two
C             consecutive points in which the output of CRT was
C             recorded.
C     NPTSE...Number of points along the ray, where all the quantities
C             have been checked.
C     LPTISO..LPTISO=1 in case that all already checked points were
C             isotropic.
C
C-----------------------------------------------------------------------
C
      IF (NPTS.NE.2) THEN
C       CRTPFA-14
        CALL ERROR('CRTPFA-14: Strange number of points in PFINTP')
C       This error should not appear, PFINTP should be called
C       for a subinterval formed by two points along the ray.
      ENDIF
C
      IF (IFIND(1).EQ.1) THEN
C       Initial point of the ray, calculating eigenvalues and
C       eigenvectors.
C       Reading the material parameters in the initial point:
        CALL PARM3(IFICB(1),RFCOOR(1,1),AA,RHO,QQ)
C       Computing eigenvectors and eigenvalues of
C       the Christoffel matrix in the initial point:
        CALL WACHRI(RFSLO1(1),RFSLO2(1),RFSLO3(1),
     *   AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7),
     *   AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14),
     *   AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21),
     *   RFGP(1),RFGS1(1),RFGS2(1),EE,DER)
C       Storing the eigenvectors of the Christoffel matrix
        RFEVEC(1,0)=EE(4)
        RFEVEC(2,0)=EE(5)
        RFEVEC(3,0)=EE(6)
        RFEVEC(4,0)=EE(7)
        RFEVEC(5,0)=EE(8)
        RFEVEC(6,0)=EE(9)
        RFEVEC(7,0)=EE(1)
        RFEVEC(8,0)=EE(2)
        RFEVEC(9,0)=EE(3)
        RFEVEC(1,1)=EE(4)
        RFEVEC(2,1)=EE(5)
        RFEVEC(3,1)=EE(6)
        RFEVEC(4,1)=EE(7)
        RFEVEC(5,1)=EE(8)
        RFEVEC(6,1)=EE(9)
        RFEVEC(7,1)=EE(1)
        RFEVEC(8,1)=EE(2)
        RFEVEC(9,1)=EE(3)
        IF (ABS(RFGS2(1)-RFGS1(1)).LT.0.00001) THEN
C         Isotropic case, no projection of eigenvectors:
          RFES1P(1,1)=1.
          RFES1P(2,1)=0.
          RFES2P(1,1)=0.
          RFES2P(2,1)=1.
          LPTISO=1
        ELSE
C         Projecting the qS eigenvectors to the plane
C         perpendicular to the ray:
          CALL WAPROJ(RFSLO1(1),RFSLO2(1),RFSLO3(1),
     *                RFPOL1(1),RFPOL2(1),RFPOL3(1),
     *                EE(4),EE(5),EE(6),EE(7),EE(8),EE(9),
     *                RFES1P(1,1),RFES1P(2,1),RFES2P(1,1),RFES2P(2,1))
          LPTISO=0
        ENDIF
        RFANGL(1)=0.
        RFTINT(1)=0.
      ENDIF
C
C     Quantities for future possible interpolation:
      TTOLD=RFTT(1)
C     Coordinates:
      XOLD(1)=RFCOOR(1,1)
      XOLD(2)=RFCOOR(2,1)
      XOLD(3)=RFCOOR(3,1)
C     dx/dt=p*v**2=p/sum(pi pi):
      AUX=(RFSLO1(1)**2+RFSLO2(1)**2+RFSLO3(1)**2)
      DXOLD(1)=RFSLO1(1)/AUX
      DXOLD(2)=RFSLO2(1)/AUX
      DXOLD(3)=RFSLO3(1)/AUX
C     Slowness:
      POLD(1)=RFSLO1(1)
      POLD(2)=RFSLO2(1)
      POLD(3)=RFSLO3(1)
C     dpi/dt=-1/v*dv/dxi:
      AUX=SQRT(AUX)
      DPOLD(1)=-RFVS1(1)/AUX
      DPOLD(2)=-RFVS2(1)/AUX
      DPOLD(3)=-RFVS3(1)/AUX
C     Polarization:
      EOLD(1)=RFPOL1(1)
      EOLD(2)=RFPOL2(1)
      EOLD(3)=RFPOL3(1)
C     dei/dt=v*sum(dv/dxk ek)*pi
      AUX=AUX*
     *(RFVS1(1)*RFPOL1(1)+RFVS2(1)*RFPOL2(1)+RFVS3(1)*RFPOL3(1))
      DEOLD(1)=AUX*RFSLO1(1)
      DEOLD(2)=AUX*RFSLO2(1)
      DEOLD(3)=AUX*RFSLO3(1)
C
      ROLDP( 1)=RFTT(1)
      ROLDP( 2)=RFCOOR(1,1)
      ROLDP( 3)=RFCOOR(2,1)
      ROLDP( 4)=RFCOOR(3,1)
      ROLDP( 5)=RFSLO1(1)
      ROLDP( 6)=RFSLO2(1)
      ROLDP( 7)=RFSLO3(1)
      ROLDP( 8)=RFPOL1(1)
      ROLDP( 9)=RFPOL2(1)
      ROLDP(10)=RFPOL3(1)
      ROLDP(11)=RFGP(1)
      ROLDP(12)=RFGS1(1)
      ROLDP(13)=RFGS2(1)
      ROLDP(14)=RFES1P(1,1)
      ROLDP(15)=RFES1P(2,1)
      ROLDP(16)=RFES2P(1,1)
      ROLDP(17)=RFES2P(2,1)
      ROLDP(18)=RFEVEC(1,1)
      ROLDP(19)=RFEVEC(2,1)
      ROLDP(20)=RFEVEC(3,1)
      ROLDP(21)=RFEVEC(4,1)
      ROLDP(22)=RFEVEC(5,1)
      ROLDP(23)=RFEVEC(6,1)
      ROLDP(24)=RFEVEC(7,1)
      ROLDP(25)=RFEVEC(8,1)
      ROLDP(26)=RFEVEC(9,1)
      ROLDP(27)=RFVS(1)
      ROLDP(28)=RFVS1(1)
      ROLDP(29)=RFVS2(1)
      ROLDP(30)=RFVS3(1)
      ROLDP(31)=RFQMAT(1,1)
      ROLDP(32)=RFQMAT(2,1)
      ROLDP(33)=RFQMAT(3,1)
      ROLDP(34)=RFQMAT(4,1)
      ROLDP(35)=RFPMAT(1,1)
      ROLDP(36)=RFPMAT(2,1)
      ROLDP(37)=RFPMAT(3,1)
      ROLDP(38)=RFPMAT(4,1)
      ROLDP(39)=RFH( 1,1)
      ROLDP(40)=RFH( 2,1)
      ROLDP(41)=RFH( 3,1)
      ROLDP(42)=RFH( 4,1)
      ROLDP(43)=RFH( 5,1)
      ROLDP(44)=RFH( 6,1)
      ROLDP(45)=RFH( 7,1)
      ROLDP(46)=RFH( 8,1)
      ROLDP(47)=RFH( 9,1)
      ROLDP(48)=RFH(10,1)
      ROLDP(49)=RFH(11,1)
      ROLDP(50)=RFH(12,1)
      ROLDP(51)=RFH(13,1)
      ROLDP(52)=RFH(14,1)
      ROLDP(53)=RFH(15,1)
      ROLDP(54)=RFH(16,1)
      ROLDP(55)=RFH(17,1)
      ROLDP(56)=RFH(18,1)
      NNEWP=0
      KNEWP(1)=0
      NPTSE=1
C
C
C     Loop along the ray:
  5   CONTINUE
      IF (NNEWP.LE.0) THEN
C       Reading the material parameters in a new point on the ray:
        IF (NNEWP.LT.0) THEN
C         CRTPFA-15
          CALL ERROR('CRTPFA-15: Negative number of points')
C         This error should not appear.
C         The number of new points should be zero or positive integer.
        ENDIF
        IF (NPTSE+1.GT.NPTS) THEN
C         CRTPFA-16
          CALL ERROR('CRTPFA-16: Point along ray missing')
C         Here we assume that we are checking the interval of points
C         NPTSE to NPTSE+1, and thus NPTS should be at least NPTSE+1
C         or more.  This error should not appear.
        ENDIF
        I1=NPTSE+1
        RNEWP( 1,1)=RFTT(I1)
        RNEWP( 2,1)=RFCOOR(1,I1)
        RNEWP( 3,1)=RFCOOR(2,I1)
        RNEWP( 4,1)=RFCOOR(3,I1)
        RNEWP( 5,1)=RFSLO1(I1)
        RNEWP( 6,1)=RFSLO2(I1)
        RNEWP( 7,1)=RFSLO3(I1)
        RNEWP( 8,1)=RFPOL1(I1)
        RNEWP( 9,1)=RFPOL2(I1)
        RNEWP(10,1)=RFPOL3(I1)
        RNEWP(27,1)=RFVS(I1)
        RNEWP(28,1)=RFVS1(I1)
        RNEWP(29,1)=RFVS2(I1)
        RNEWP(30,1)=RFVS3(I1)
        RNEWP(31,1)=RFQMAT(1,I1)
        RNEWP(32,1)=RFQMAT(2,I1)
        RNEWP(33,1)=RFQMAT(3,I1)
        RNEWP(34,1)=RFQMAT(4,I1)
        RNEWP(35,1)=RFPMAT(1,I1)
        RNEWP(36,1)=RFPMAT(2,I1)
        RNEWP(37,1)=RFPMAT(3,I1)
        RNEWP(38,1)=RFPMAT(4,I1)
        RNEWP(39,1)=RFH( 1,I1)
        RNEWP(40,1)=RFH( 2,I1)
        RNEWP(41,1)=RFH( 3,I1)
        RNEWP(42,1)=RFH( 4,I1)
        RNEWP(43,1)=RFH( 5,I1)
        RNEWP(44,1)=RFH( 6,I1)
        RNEWP(45,1)=RFH( 7,I1)
        RNEWP(46,1)=RFH( 8,I1)
        RNEWP(47,1)=RFH( 9,I1)
        RNEWP(48,1)=RFH(10,I1)
        RNEWP(49,1)=RFH(11,I1)
        RNEWP(50,1)=RFH(12,I1)
        RNEWP(51,1)=RFH(13,I1)
        RNEWP(52,1)=RFH(14,I1)
        RNEWP(53,1)=RFH(15,I1)
        RNEWP(54,1)=RFH(16,I1)
        RNEWP(55,1)=RFH(17,I1)
        RNEWP(56,1)=RFH(18,I1)
        NNEWP=1
        CALL PARM3(IFICB(I1),RNEWP(2,1),AA,RHO,QQ)
C       Computing the Christoffel matrix in the new point:
        CALL WACHRI(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1),
     *  AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7),
     *  AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14),
     *  AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),
     *  AA(1,21),RNEWP(11,1),RNEWP(12,1),RNEWP(13,1),EE,DER)
C       Storing the eigenvectors of the Christoffel matrix
        RNEWP(18,1)=EE(4)
        RNEWP(19,1)=EE(5)
        RNEWP(20,1)=EE(6)
        RNEWP(21,1)=EE(7)
        RNEWP(22,1)=EE(8)
        RNEWP(23,1)=EE(9)
        RNEWP(24,1)=EE(1)
        RNEWP(25,1)=EE(2)
        RNEWP(26,1)=EE(3)
C       Projecting the qS eigenvectors to the plane
C       perpendicular to the ray:
        CALL WAPROJ(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1),RNEWP(8,1),
     *  RNEWP(9,1),RNEWP(10,1),EE(4),EE(5),EE(6),EE(7),EE(8),EE(9),
     *  RNEWP(14,1),RNEWP(15,1),RNEWP(16,1),RNEWP(17,1))
C
C       Quantities for future possible interpolation:
        TTNEW=RFTT(I1)
C       Coordinates:
        XNEW(1)=RFCOOR(1,I1)
        XNEW(2)=RFCOOR(2,I1)
        XNEW(3)=RFCOOR(3,I1)
C       dx/dt=p*v**2=p/sum(pi pi):
        AUX=(RFSLO1(I1)**2+RFSLO2(I1)**2+RFSLO3(I1)**2)
        DXNEW(1)=RFSLO1(I1)/AUX
        DXNEW(2)=RFSLO2(I1)/AUX
        DXNEW(3)=RFSLO3(I1)/AUX
C       Slowness:
        PNEW(1)=RFSLO1(I1)
        PNEW(2)=RFSLO2(I1)
        PNEW(3)=RFSLO3(I1)
C       dpi/dt=-1/v*dv/dxi:
        AUX=SQRT(AUX)
        DPNEW(1)=-RFVS1(I1)/AUX
        DPNEW(2)=-RFVS2(I1)/AUX
        DPNEW(3)=-RFVS3(I1)/AUX
C       Polarization:
        ENEW(1)=RFPOL1(I1)
        ENEW(2)=RFPOL2(I1)
        ENEW(3)=RFPOL3(I1)
C       dei/dt=v*sum(dv/dxk ek)*pi
        AUX=AUX*
     *  (RFVS1(I1)*RFPOL1(I1)+RFVS2(I1)*RFPOL2(I1)+RFVS3(I1)*RFPOL3(I1))
        DENEW(1)=AUX*RFSLO1(I1)
        DENEW(2)=AUX*RFSLO2(I1)
        DENEW(3)=AUX*RFSLO3(I1)
C
        IF (ABS(RNEWP(12,1)-RNEWP(13,1)).LT.0.00001) THEN
C         New point is isotropic, qS eigenvalues are the same.
C         No change in eigenvalues and eigenvectors:
          RNEWP(14,1)=ROLDP(14)
          RNEWP(15,1)=ROLDP(15)
          RNEWP(16,1)=ROLDP(16)
          RNEWP(17,1)=ROLDP(17)
          RNEWP(18,1)=ROLDP(18)
          RNEWP(19,1)=ROLDP(19)
          RNEWP(20,1)=ROLDP(20)
          RNEWP(21,1)=ROLDP(21)
          RNEWP(22,1)=ROLDP(22)
          RNEWP(23,1)=ROLDP(23)
          RNEWP(24,1)=ROLDP(24)
          RNEWP(25,1)=ROLDP(25)
          RNEWP(26,1)=ROLDP(26)
          DELTFI=0.
          GOTO 20
        ENDIF
        IF (ABS(ROLDP(12)-ROLDP(13)).LT.0.00001) THEN
C         Old point is isotropic, new point is anisotropic:
          IF (LPTISO.EQ.1) THEN
C           New point is the first anisotropic point on the ray,
C           angular change is not to be computed:
            DELTFI=0.
            LPTISO=0
            GOTO 20
          ENDIF
        ENDIF
        IF (RFTT(NPTSE).EQ.RFTT(NPTSE+1)) THEN
C         The ray is crossing an interface.
          DELTFI=0.
          GOTO 20
        ENDIF
      ENDIF
C
C     Computing the angular change in eigenvectors, adding new points
C     on the ray if necessary:
      DO 10, I1=1,NNEWP
        CALL WAVECR(ROLDP(18),ROLDP(21),RNEWP(18,I1),RNEWP(21,I1),
     *              RNEWP(12,I1),RNEWP(13,I1))
        DELTFI=WAANGL(ROLDP(18),ROLDP(21),RNEWP(18,I1),RNEWP(21,I1))
        DEFI=1./SQRT(RNEWP(12,I1))-1./SQRT(RNEWP(13,I1))
        DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13))
        DEFI=DEFI*PI*SINGLF*RFTT(NPTS)/6.
        DEFI=ABS(DEFI)
        IF (DEFI*PI/8..GT.ERRWAN) THEN
          DEFI=ERRWAN/DEFI
        ELSE
          DEFI=PI/8.
        END IF
        IF (ABS(DELTFI).LE.DEFI) THEN
C         Angular change is less than prescribed limit, the I1-th point
C         of array RNEWP will be used as the next point on the ray:
          NNEWP=I1
          GOTO 20
        ENDIF
   10 CONTINUE
C     Angular change is greater than prescribed limit for all the points
C     of array RNEWP.
   15 CONTINUE
C     Loop for adding new points on the ray until the angular change is
C     small enough:
        IF (NNEWP.GE.MNEWP-2) THEN
C         CRTPFA-17
          CALL ERROR('CRTPFA-17: Maximum number of divisions exceeded')
C         The angular change of eigenvectors in two consecutive points
C         is too big. More divisions than fits into the memory is
C         needed to keep the change under the prescribed limit. Try to
C         decrease the parameter
C         STORE,
C         or to increase the dimension MNEWP or the angular change DEFI.
        ENDIF
C       Adding new point to the ray:
        KNEWP(NNEWP)=KNEWP(NNEWP)+1
        NNEWP=NNEWP+1
        IF (NNEWP.GT.MNEWP) THEN
C         CRTPFA-18
          CALL ERROR('CRTPFA-18: Array KNEWP too small')
C         This error should not appear, error17 should appear instead.
        ENDIF
C       Travel time:
        RNEWP(1,NNEWP)=(ROLDP(1)+RNEWP(1,NNEWP-1))*0.5
C       Coordinates:
        CALL HIVD2(3,TTOLD,XOLD,DXOLD,TTNEW,XNEW,DXNEW,
     *  RNEWP(1,NNEWP),RNEWP(2,NNEWP),DER)
C       Slowness vector:
        CALL HIVD2(3,TTOLD,POLD,DPOLD,TTNEW,PNEW,DPNEW,
     *  RNEWP(1,NNEWP),RNEWP(5,NNEWP),DER)
C       Polarization vector:
        CALL HIVD2(3,TTOLD,EOLD,DEOLD,TTNEW,ENEW,DENEW,
     *  RNEWP(1,NNEWP),RNEWP(8,NNEWP),DER)
C       Matrices Q and P, basis vectors of RCCS:
        DO 17, I1=31,56
          RNEWP(I1,NNEWP)=(ROLDP(I1)+RNEWP(I1,NNEWP-1))*0.5
  17    CONTINUE
C       Material parameters:
        CALL PARM2(IFICB(NPTSE+1),RNEWP(2,NNEWP),UP,US,RHO,QP,QS)
        CALL VELOC(IFICB(NPTSE+1),UP,US,QP,QS,VP,VS,VD,QL)
        RNEWP(27,NNEWP)=VD(1)
        RNEWP(28,NNEWP)=VD(2)
        RNEWP(29,NNEWP)=VD(3)
        RNEWP(30,NNEWP)=VD(4)
        CALL PARM3(IFICB(NPTSE+1),RNEWP(2,NNEWP),AA,RHO,QQ)
C       Christoffel matrix and eigenvalues:
        CALL WACHRI(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP),
     *  AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7),
     *  AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14),
     *  AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21),
     *  RNEWP(11,NNEWP),RNEWP(12,NNEWP),RNEWP(13,NNEWP),EE,DER)
C       Storing the eigenvectors of the Christoffel matrix
        RNEWP(18,NNEWP)=EE(4)
        RNEWP(19,NNEWP)=EE(5)
        RNEWP(20,NNEWP)=EE(6)
        RNEWP(21,NNEWP)=EE(7)
        RNEWP(22,NNEWP)=EE(8)
        RNEWP(23,NNEWP)=EE(9)
        RNEWP(24,NNEWP)=EE(1)
        RNEWP(25,NNEWP)=EE(2)
        RNEWP(26,NNEWP)=EE(3)
C       Projection of the qS eigenvectors to the plane
C       perpendicular to the ray:
        CALL WAPROJ(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP),
     *  RNEWP(8,NNEWP),RNEWP(9,NNEWP),RNEWP(10,NNEWP),
     *  EE(4),EE(5),EE(6),EE(7),EE(8),EE(9),
     *  RNEWP(14,NNEWP),RNEWP(15,NNEWP),RNEWP(16,NNEWP),RNEWP(17,NNEWP))
C       Index of the division:
        KNEWP(NNEWP)=KNEWP(NNEWP-1)
C
        IF (ABS(RNEWP(12,NNEWP)-RNEWP(13,NNEWP)).LT.0.00001) THEN
C         Isotropic case, qS eigenvalues are the same,
C         no change in eigenvectors:
          RNEWP(14,NNEWP)=ROLDP(14)
          RNEWP(15,NNEWP)=ROLDP(15)
          RNEWP(16,NNEWP)=ROLDP(16)
          RNEWP(17,NNEWP)=ROLDP(17)
          RNEWP(18,NNEWP)=ROLDP(18)
          RNEWP(19,NNEWP)=ROLDP(19)
          RNEWP(20,NNEWP)=ROLDP(20)
          RNEWP(21,NNEWP)=ROLDP(21)
          RNEWP(22,NNEWP)=ROLDP(22)
          RNEWP(23,NNEWP)=ROLDP(23)
          RNEWP(24,NNEWP)=ROLDP(24)
          RNEWP(25,NNEWP)=ROLDP(25)
          RNEWP(26,NNEWP)=ROLDP(26)
          DELTFI=0.
        ELSE
          CALL WAVECR(ROLDP(18),ROLDP(21),
     *                RNEWP(18,NNEWP),RNEWP(21,NNEWP),
     *                RNEWP(12,NNEWP),RNEWP(13,NNEWP))
          DELTFI=WAANGL(ROLDP(18),ROLDP(21),
     *                  RNEWP(18,NNEWP),RNEWP(21,NNEWP))
          DEFI=1./SQRT(RNEWP(12,NNEWP))-1./SQRT(RNEWP(13,NNEWP))
          DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13))
          DEFI=DEFI*PI*SINGLF*RFTT(NPTS)/6.
          DEFI=ABS(DEFI)
          IF (DEFI*PI/8..GT.ERRWAN) THEN
            DEFI=ERRWAN/DEFI
          ELSE
            DEFI=PI/8.
          END IF
        ENDIF
        IF (ABS(DELTFI).LE.DEFI) THEN
C         Angular change is less than prescribed limit, this point
C         of array RNEWP will be used as the next point on the ray:
          GOTO 20
        ELSE
C         Angular change is greater than prescribed limit, adding
C         a new point to the ray:
          GOTO 15
        ENDIF
C     End of the loop.
C
   20 CONTINUE
C     Angular change DELTFI for points ROLDP, RNEWP(i,NNEWP) is less
C     than prescribed limit. Recording the computed quantities.
      NPTSE=NPTSE+1
      IF (NNEWP.NE.1) THEN
C       The new point was computed by interpolation.
C       Shifting arrays RF* and IF*:
        CALL PFSHIF(NPTSE,NPTSE+1)
C       Recording interpolated quantities:
        IFIND(   NPTSE)=0
        RFTT(    NPTSE)=RNEWP( 1,NNEWP)
        RFCOOR(1,NPTSE)=RNEWP( 2,NNEWP)
        RFCOOR(2,NPTSE)=RNEWP( 3,NNEWP)
        RFCOOR(3,NPTSE)=RNEWP( 4,NNEWP)
        RFSLO1(  NPTSE)=RNEWP( 5,NNEWP)
        RFSLO2(  NPTSE)=RNEWP( 6,NNEWP)
        RFSLO3(  NPTSE)=RNEWP( 7,NNEWP)
        RFPOL1(  NPTSE)=RNEWP( 8,NNEWP)
        RFPOL2(  NPTSE)=RNEWP( 9,NNEWP)
        RFPOL3(  NPTSE)=RNEWP(10,NNEWP)
        IFICB(   NPTSE)=IFICB(NPTSE-1)
        RFVS(    NPTSE)=RNEWP(27,NNEWP)
        RFVS1(   NPTSE)=RNEWP(28,NNEWP)
        RFVS2(   NPTSE)=RNEWP(29,NNEWP)
        RFVS3(   NPTSE)=RNEWP(30,NNEWP)
        RFQMAT(1,NPTSE)=RNEWP(31,NNEWP)
        RFQMAT(2,NPTSE)=RNEWP(32,NNEWP)
        RFQMAT(3,NPTSE)=RNEWP(33,NNEWP)
        RFQMAT(4,NPTSE)=RNEWP(34,NNEWP)
        RFPMAT(1,NPTSE)=RNEWP(35,NNEWP)
        RFPMAT(2,NPTSE)=RNEWP(36,NNEWP)
        RFPMAT(3,NPTSE)=RNEWP(37,NNEWP)
        RFPMAT(4,NPTSE)=RNEWP(38,NNEWP)
        IFISRF(  NPTSE)=0
        RFH( 1,NPTSE)=RNEWP(39,NNEWP)
        RFH( 2,NPTSE)=RNEWP(40,NNEWP)
        RFH( 3,NPTSE)=RNEWP(41,NNEWP)
        RFH( 4,NPTSE)=RNEWP(42,NNEWP)
        RFH( 5,NPTSE)=RNEWP(43,NNEWP)
        RFH( 6,NPTSE)=RNEWP(44,NNEWP)
        RFH( 7,NPTSE)=RNEWP(45,NNEWP)
        RFH( 8,NPTSE)=RNEWP(46,NNEWP)
        RFH( 9,NPTSE)=RNEWP(47,NNEWP)
        RFH(10,NPTSE)=RNEWP(48,NNEWP)
        RFH(11,NPTSE)=RNEWP(49,NNEWP)
        RFH(12,NPTSE)=RNEWP(50,NNEWP)
        RFH(13,NPTSE)=RNEWP(51,NNEWP)
        RFH(14,NPTSE)=RNEWP(52,NNEWP)
        RFH(15,NPTSE)=RNEWP(53,NNEWP)
        RFH(16,NPTSE)=RNEWP(54,NNEWP)
        RFH(17,NPTSE)=RNEWP(55,NNEWP)
        RFH(18,NPTSE)=RNEWP(56,NNEWP)
      ENDIF
C     Recording quantities for computation of anisotropic corrections:
      RFGP(NPTSE)=RNEWP(11,NNEWP)
      RFGS1(NPTSE)=RNEWP(12,NNEWP)
      RFGS2(NPTSE)=RNEWP(13,NNEWP)
      RFES1P(1,NPTSE)=RNEWP(14,NNEWP)
      RFES1P(2,NPTSE)=RNEWP(15,NNEWP)
      RFES2P(1,NPTSE)=RNEWP(16,NNEWP)
      RFES2P(2,NPTSE)=RNEWP(17,NNEWP)
      RFANGL(NPTSE)=DELTFI
      RFTINT(NPTSE)=
     *     (0.5/SQRT(RNEWP(13,NNEWP))-0.5/SQRT(RNEWP(12,NNEWP))
     *     +0.5/SQRT(ROLDP(13))      -0.5/SQRT(ROLDP(12)))
     *    *(RNEWP(1,NNEWP)-ROLDP(1))
      RFEVEC(1,NPTSE)=RNEWP(18,NNEWP)
      RFEVEC(2,NPTSE)=RNEWP(19,NNEWP)
      RFEVEC(3,NPTSE)=RNEWP(20,NNEWP)
      RFEVEC(4,NPTSE)=RNEWP(21,NNEWP)
      RFEVEC(5,NPTSE)=RNEWP(22,NNEWP)
      RFEVEC(6,NPTSE)=RNEWP(23,NNEWP)
      RFEVEC(7,NPTSE)=RNEWP(24,NNEWP)
      RFEVEC(8,NPTSE)=RNEWP(25,NNEWP)
      RFEVEC(9,NPTSE)=RNEWP(26,NNEWP)
C
      IF (NPTSE.LE.NPTS) THEN
C       Continuing with the next point on the ray:
        DO 100, I1=1,MQANT
          ROLDP(I1)=RNEWP(I1,NNEWP)
  100   CONTINUE
        KNEWP(NNEWP)=0
        NNEWP=NNEWP-1
        IF (NNEWP.LE.0) THEN
          TTOLD=TTNEW
          DO 101, I1=1,3
            XOLD(I1)=XNEW(I1)
            DXOLD(I1)=DXNEW(I1)
            POLD(I1)=PNEW(I1)
            DPOLD(I1)=DPNEW(I1)
            EOLD(I1)=ENEW(I1)
            DEOLD(I1)=DENEW(I1)
  101     CONTINUE
        ENDIF
        IF (NPTSE.LT.NPTS) GOTO 5
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE PFPIGE(PIGET1,PIGET2,PIGR1,PIGI1,PIGR2,PIGI2,DPIGD)
C
C-----------------------------------------------------------------------
C Subroutine to calculate the propagator matrix
C
      REAL PIGET1,PIGET2,PIGR1(2,2),PIGI1(2,2),PIGR2(2,2),PIGI2(2,2),
     *     DPIGD
C No input.
C Output:
C     PIGET1,PIGET2 ... Travel times of two coupling rays.
C     PIGR1,PIGI1,PIGR2,PIGI2 ... Real and imaginary part of the
C             propagator matrix of the first and second coupling ray.
C     DPIGD.. Derivative of travel-time difference D with respect
C             to the reference travel time.
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C ...........................
C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points
C of rays and to store the input SEP parameters:
      INCLUDE 'crtpfa.inc'
C ...........................
C External functions directly called:
      EXTERNAL ERROR,WARN,AP21,APYYF,PFPERT,WAMAT,LENGTH
      INTEGER LENGTH
C     ERROR,WARN ... File error.for.
C     AP21,APYYF ... File ap.for.
C     PFPERT..This file.
C     WAMAT.. File wanpfa.for.
C     LENGTH..File length.for.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER I1,I2,I3,IPTTMP,ISRFFT,NYFTMP,ISINGL,NSINGL
      REAL PI
      PARAMETER (PI=3.1415926)
      REAL PIGR(2,2),PIGI(2,2),DPIGR(2,2),DPIGI(2,2),
     *     PIGRA(2,2),PIGIA(2,2),DPIGRA(2,2),DPIGIA(2,2),
     *     PIGRB(2,2),PIGIB(2,2),DPIG,C,S,
     *     DPIGRO(2,2),DPIGIO(2,2),DPIGO,
     *     DPIGRN(2,2),DPIGIN(2,2),DPIGN,
     *     TTCOR,DTAU(7),DTLIN,STEP,FREQ,YI1TMP,YF1TMP,GAMA,
     *     AUX,AUX0,AUX1,AUX2,AUX3,AUX4
      CHARACTER*260 TXTERR
      LOGICAL LPWAVE,LSWAVE
      SAVE LPWAVE,LSWAVE,TTCOR,DTAU,DTLIN,STEP,FREQ,PIGR,PIGI,
     *     DPIGR,DPIGI,DPIGRO,DPIGIO,DPIGO
C-----------------------------------------------------------------------
      IF (IFIND(1).EQ.1) THEN
C       Initial point of the ray:
        LPWAVE=.TRUE.
        LSWAVE=.TRUE.
        TTCOR=0.
        DTAU(1)=0.
        DTAU(3)=0.
        DTLIN=0.
        DO 10, I1=1,NPTS
          IF (IFIND(I1).EQ.2) THEN
            STEP=RFTT(I1)-RFTT(1)
          ENDIF
  10    CONTINUE
        FREQ=SINGLF
        DPIGR(1,1)=0.
        DPIGR(2,1)=0.
        DPIGR(1,2)=0.
        DPIGR(2,2)=0.
        DPIGI(1,1)=0.
        DPIGI(2,1)=0.
        DPIGI(1,2)=0.
        DPIGI(2,2)=0.
        PIGR(1,1)=1.
        PIGR(2,1)=0.
        PIGR(1,2)=0.
        PIGR(2,2)=1.
        PIGI(1,1)=0.
        PIGI(2,1)=0.
        PIGI(1,2)=0.
        PIGI(2,2)=0.
      ENDIF
C     Computing the values of travel time corrections along the ray:
      DO 30, I1=1,NPTS-1
        I2=I1+1
        IF (IFICB(I2).GT.0) THEN
C         CRTPFA-19
          CALL ERROR('CRTPFA-19: P wave encountered')
C         P wave is not yet coded.
CP        TTCOR=TTCOR+(1./SQRT(RFGP(I2))+1./SQRT(RFGP(I1)))*0.5*
CP   *          (RFTT(I2)-RFTT(I1))
CP        LSWAVE=.FALSE.
        ELSE
C         S wave:
            TTCOR=TTCOR+(0.25/SQRT(RFGS1(I2))  +0.25/SQRT(RFGS2(I2))
     *                  +0.25/SQRT(RFGS1(I1))  +0.25/SQRT(RFGS2(I1)))
     *                 *(RFTT(I2)  -RFTT(I1))
          LPWAVE=.FALSE.
        ENDIF
   30 CONTINUE
C
CP    IF (LPWAVE) THEN
CP      P wave along the ray from the source to the current point:
CP      Y(1)=TTCOR+YI(1)
CP      CALL AP21(0,GREEN)
CP      NQ=32
CP      GOTO 91 (to the end of subroutine)
CP    ENDIF
C
      IF (LSWAVE.AND.(KQIRAY.NE.0)) THEN
C       S wave along the ray from the source to the current point,
C       calculating second-order travel-time perturbations:
        TTCOR=TTCOR-0.5*(DTAU(3)+DTAU(1))
        CALL APYYF
        IPTTMP=IPT
        YI1TMP=YI(1)
        NYFTMP=NYF
        ISRFFT=ISRFF
        YF1TMP=YF(1)
C
        YI(1) =RFTT(0)
        NYF   =1
        DO 95, I1=1,NPTS
          IF (IFIND(I1).EQ.0) THEN
            STEP=0.
            GOTO 96
          ENDIF
  95    CONTINUE
  96    CONTINUE
        IF (IFIND(1).EQ.1) THEN
          I2=1
        ELSE
          I2=2
        ENDIF
        DO 100, I1=I2,NPTS
          IF (I1.NE.1) DTLIN=DTLIN+RFTINT(I1)
          IPT   =IFIND(I1)-1
          ISRFF =IFISRF(I1)
          YF(1) =RFTT(I1)
          RFTINT(I1)=RFTINT(I1)+DTAU(1)-DTAU(3)
          CALL PFPERT(I1,STEP,DTAU)
          RFTINT(I1)=RFTINT(I1)+DTAU(3)-DTAU(1)
 100    CONTINUE
C
        IPT   =IPTTMP
        YI(1) =YI1TMP
        NYF   =NYFTMP
        ISRFF =ISRFFT
        YF(1) =YF1TMP
        CALL APYYF
C
        TTCOR=TTCOR+0.5*(DTAU(3)+DTAU(1))
      ELSEIF ((.NOT.LSWAVE).AND.(KQIRAY.NE.0)) THEN
C       CRTPFA-20
        WRITE(TXTERR,'(2A,1I6,A)')
     *  'CRTPFA-20: Second-order perturbations of anisotropic travel ',
     *  'times may not be calculated for the ray ',IRAY,
     *  ', the ray crossed an interface.'
        CALL WARN(TXTERR(1:LENGTH(TXTERR)))
C       The perturbations may be calculated only along the rays
C       which do not cross any interface.
      ENDIF
C
C     S wave in some part of the ray:
C     Computing the propagator matrix PiGe along the ray:
C     Loop along points on the ray:
      DO 40, I1=2,NPTS
        I3=I1-1
        IF ((RFTT(I1).EQ.RFTT(I3)).AND.(IFICB(I1).NE.IFICB(I3))) THEN
C         The ray crossed an interface:
C         CRTPFA-21
          CALL ERROR('CRTPFA-21: Ray crossed an interface.')
C         Not yet coded.
        ENDIF
C       Element of the ray inside a complex block:
        IF (IFICB(I1).LT.0) THEN
C         S wave:
          GAMA=PI*FREQ*RFTINT(I1)
          AUX0=SQRT(RFANGL(I1)**2 + GAMA**2)
          AUX1=COS(AUX0)
          IF (AUX0.EQ.0.) THEN
            AUX2=1.
            AUX3=-1./3.
          ELSE
            AUX2=SIN(AUX0)/AUX0
            AUX3=(COS(AUX0)-AUX2)/AUX0**2
          ENDIF
          AUX4=GAMA/(2.*PI*FREQ)
C         Matrix for this step along the ray:
C         R22, Eq.(49):
          PIGRA(1,1)= AUX1
          PIGRA(2,1)=-RFANGL(I1)*AUX2
          PIGRA(1,2)=-PIGRA(2,1)
          PIGRA(2,2)= AUX1
          PIGIA(1,1)=-GAMA*AUX2
          PIGIA(2,1)= 0.
          PIGIA(1,2)= 0.
          PIGIA(2,2)=-PIGIA(1,1)
C         R22, Eq.(56):
          DPIGRA(1,1)=-AUX2*GAMA*AUX4
          DPIGRA(2,1)=-AUX3*RFANGL(I1)*GAMA*AUX4
          DPIGRA(1,2)=-DPIGRA(2,1)
          DPIGRA(2,2)= DPIGRA(1,1)
          DPIGIA(1,1)=-(AUX3*GAMA**2+AUX2)*AUX4
          DPIGIA(2,1)= 0.
          DPIGIA(1,2)= 0.
          DPIGIA(2,2)=-DPIGIA(1,1)
        ELSE
C         P wave:
          GOTO 39
        ENDIF
C       Matrix for all the steps along the ray to current point:
C       R22, Eq.(53):
        CALL WAMAT(PIGRA,PIGIA,DPIGR,DPIGI)
        PIGRB(1,1)=PIGR(1,1)
        PIGRB(2,1)=PIGR(2,1)
        PIGRB(1,2)=PIGR(1,2)
        PIGRB(2,2)=PIGR(2,2)
        PIGIB(1,1)=PIGI(1,1)
        PIGIB(2,1)=PIGI(2,1)
        PIGIB(1,2)=PIGI(1,2)
        PIGIB(2,2)=PIGI(2,2)
        CALL WAMAT(DPIGRA,DPIGIA,PIGRB,PIGIB)
        DPIGR(1,1)=DPIGR(1,1)+PIGRB(1,1)
        DPIGR(2,1)=DPIGR(2,1)+PIGRB(2,1)
        DPIGR(1,2)=DPIGR(1,2)+PIGRB(1,2)
        DPIGR(2,2)=DPIGR(2,2)+PIGRB(2,2)
        DPIGI(1,1)=DPIGI(1,1)+PIGIB(1,1)
        DPIGI(2,1)=DPIGI(2,1)+PIGIB(2,1)
        DPIGI(1,2)=DPIGI(1,2)+PIGIB(1,2)
        DPIGI(2,2)=DPIGI(2,2)+PIGIB(2,2)
C       R22, Eq.(47):
        CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI)
   39   CONTINUE
   40 CONTINUE
C
C     Single-frequency approximation
C     R22, Eq.(21):
      DPIG=DPIGR(1,1)**2+DPIGI(1,1)**2+DPIGR(2,1)**2+DPIGI(2,1)**2
      DPIG=SQRT(DPIG)
      DPIGRN(1,1)=DPIGR(1,1)
      DPIGIN(1,1)=DPIGI(1,1)
      DPIGRN(2,1)=DPIGR(2,1)
      DPIGIN(2,1)=DPIGI(2,1)
      DPIGRN(1,2)=DPIGR(1,2)
      DPIGIN(1,2)=DPIGI(1,2)
      DPIGRN(2,2)=DPIGR(2,2)
      DPIGIN(2,2)=DPIGI(2,2)
      DPIGN=DPIG
      IF (DPIG.EQ.0.) THEN
        DPIGRN(1,1)= PIGI(1,1)
        DPIGIN(1,1)=-PIGR(1,1)
        DPIGRN(2,1)=-PIGI(2,1)
        DPIGIN(2,1)= PIGR(2,1)
        DPIGRN(1,2)= PIGI(1,2)
        DPIGIN(1,2)=-PIGR(1,2)
        DPIGRN(2,2)=-PIGI(2,2)
        DPIGIN(2,2)= PIGR(2,2)
        DPIGN=1.
      END IF
      IF (IFIND(1).NE.1) THEN
        AUX=DPIGRN(1,1)*DPIGRO(1,1)+DPIGIN(1,1)*DPIGIO(1,1)
     *     +DPIGRN(2,1)*DPIGRO(2,1)+DPIGIN(2,1)*DPIGIO(2,1)
        DPIG= SIGN(DPIG,AUX*DPIGO)
        DPIGN=SIGN(DPIGN,AUX*DPIGO)
      ENDIF
      PIGIA(1,1)=PIGI(1,1)
      PIGRA(1,1)=PIGR(1,1)
      PIGIA(2,1)=PIGI(2,1)
      PIGRA(2,1)=PIGR(2,1)
      PIGIA(1,2)=PIGI(1,2)
      PIGRA(1,2)=PIGR(1,2)
      PIGIA(2,2)=PIGI(2,2)
      PIGRA(2,2)=PIGR(2,2)
      NSINGL=32
C
C     Computing the propagator matrix:
      DO 80 ISINGL=0,NSINGL,32
C       Transforming the propagator matrix from eigenvectors
C       to polarization vectors:
        IF (ISINGL.EQ.0) THEN
C         R22, Eq.(41):
          PIGET1=TTCOR+YI(1)-DPIG
C         R22, Eq.(22):
          PIGIB(1,1)=(PIGIA(1,1)+DPIGRN(1,1)/DPIGN)/2.
          PIGRB(1,1)=(PIGRA(1,1)-DPIGIN(1,1)/DPIGN)/2.
          PIGIB(2,1)=(PIGIA(2,1)+DPIGRN(2,1)/DPIGN)/2.
          PIGRB(2,1)=(PIGRA(2,1)-DPIGIN(2,1)/DPIGN)/2.
          PIGIB(1,2)=(PIGIA(1,2)+DPIGRN(1,2)/DPIGN)/2.
          PIGRB(1,2)=(PIGRA(1,2)-DPIGIN(1,2)/DPIGN)/2.
          PIGIB(2,2)=(PIGIA(2,2)+DPIGRN(2,2)/DPIGN)/2.
          PIGRB(2,2)=(PIGRA(2,2)-DPIGIN(2,2)/DPIGN)/2.
          C=COS(2.*PI*SINGLF*DPIG)
          S=SIN(2.*PI*SINGLF*DPIG)
C         R22, Eq.(43):
          PIGI1(1,1)=C*PIGIB(1,1)+S*PIGRB(1,1)
          PIGR1(1,1)=C*PIGRB(1,1)-S*PIGIB(1,1)
          PIGI1(2,1)=C*PIGIB(2,1)+S*PIGRB(2,1)
          PIGR1(2,1)=C*PIGRB(2,1)-S*PIGIB(2,1)
          PIGI1(1,2)=C*PIGIB(1,2)+S*PIGRB(1,2)
          PIGR1(1,2)=C*PIGRB(1,2)-S*PIGIB(1,2)
          PIGI1(2,2)=C*PIGIB(2,2)+S*PIGRB(2,2)
          PIGR1(2,2)=C*PIGRB(2,2)-S*PIGIB(2,2)
        ELSE
C         R22, Eq.(42):
          PIGET2=TTCOR+YI(1)+DPIG
C         R22, Eq.(23):
          PIGIB(1,1)=(PIGIA(1,1)-DPIGRN(1,1)/DPIGN)/2.
          PIGRB(1,1)=(PIGRA(1,1)+DPIGIN(1,1)/DPIGN)/2.
          PIGIB(2,1)=(PIGIA(2,1)-DPIGRN(2,1)/DPIGN)/2.
          PIGRB(2,1)=(PIGRA(2,1)+DPIGIN(2,1)/DPIGN)/2.
          PIGIB(1,2)=(PIGIA(1,2)-DPIGRN(1,2)/DPIGN)/2.
          PIGRB(1,2)=(PIGRA(1,2)+DPIGIN(1,2)/DPIGN)/2.
          PIGIB(2,2)=(PIGIA(2,2)-DPIGRN(2,2)/DPIGN)/2.
          PIGRB(2,2)=(PIGRA(2,2)+DPIGIN(2,2)/DPIGN)/2.
          S=-S
C         R22, Eq.(43):
          PIGI2(1,1)=C*PIGIB(1,1)+S*PIGRB(1,1)
          PIGR2(1,1)=C*PIGRB(1,1)-S*PIGIB(1,1)
          PIGI2(2,1)=C*PIGIB(2,1)+S*PIGRB(2,1)
          PIGR2(2,1)=C*PIGRB(2,1)-S*PIGIB(2,1)
          PIGI2(1,2)=C*PIGIB(1,2)+S*PIGRB(1,2)
          PIGR2(1,2)=C*PIGRB(1,2)-S*PIGIB(1,2)
          PIGI2(2,2)=C*PIGIB(2,2)+S*PIGRB(2,2)
          PIGR2(2,2)=C*PIGRB(2,2)-S*PIGIB(2,2)
        ENDIF
   80 CONTINUE
C
C     Calculating the length correction of the reference slowness
C     vector (R23 eq.37):
      DPIGD=(PIGIA(1,1)*DPIGRN(1,1)-PIGRA(1,1)*DPIGIN(1,1)
     *      +PIGIA(1,2)*DPIGRN(1,2)-PIGRA(1,2)*DPIGIN(1,2))/DPIGN
      DPIGD=0.5*DPIGD*(1./SQRT(RFGS2(NPTS))-1./SQRT(RFGS1(NPTS)))
C
C     Saving quantities for future invocation of the subroutine:
      DPIGRO(1,1)=DPIGRN(1,1)
      DPIGRO(2,1)=DPIGRN(2,1)
      DPIGRO(1,2)=DPIGRN(1,2)
      DPIGRO(2,2)=DPIGRN(2,2)
      DPIGIO(1,1)=DPIGIN(1,1)
      DPIGIO(2,1)=DPIGIN(2,1)
      DPIGIO(1,2)=DPIGIN(1,2)
      DPIGIO(2,2)=DPIGIN(2,2)
      DPIGO      =DPIGN
C
      RETURN
      END
C
C=======================================================================
C
       SUBROUTINE PFPERT(I1,STEP,DTAU)
C
C----------------------------------------------------------------------
C Subroutine to compute second-order travel-time corrections
      INTEGER I1
      REAL    STEP,DTAU(7)
C
C Input:
C     I1 ...  Position of the point in arrays of common block PFPTS.
C     STEP... Step of the independent variable along the ray, see
C             the CRT input parameter STORE.
C Output:
C     DTAU... Second-order travel-time corrections.
C
C.......................................................................
C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points
C of rays and to store the input SEP parameters:
      INCLUDE 'crtpfa.inc'
C ...........................
C Subroutines and external functions directly called:
      EXTERNAL ERROR,AP28,PARM3,HDER,WAIJKL,WASUM,WASUM3,WASUM4,WASUM5
      REAL WASUM,WASUM3,WASUM4,WASUM5
C     ERROR.. File error.for.
C     AP28... File ap.for.
C     PARM3.. File parm.for.
C     HDER... File hder.for.
C     WAIJKL,WASUM,WASUM3,WASUM4,WASUM5 ... File wanpfa.for.
C.......................................................................
C Auxiliary storage locations:
      REAL AA(10,21),RHO,QQ(21),
     *     AAA(3,3,3,3),AAA1(3,3,3,3),AAA2(3,3,3,3),AAA3(3,3,3,3)
      REAL V0,VV0(3),P(3),E1(3),E2(3),H1(3),H2(3),H3(3),
     *     HC1(3),HC2(3),HC3(3),
     *     HHC0(3),HHC0SP(3,3),HHC0S(3,3),HHTC0S(2,2),
     *     QQT(4),PPT(4),QIT(4),GIJT(2,2),
     *     HH0(3),HHT0(3),HH1,HH2,HH3,
     *     HH1L(3),HH2L(3),HH3L(3),HH1U(3),HH2U(3),HH3U(3),
     *     HHT1L(2),HHT2L(2),HHT3L(2),HHT1U(2),HHT2U(2),HHT3U(2),
     *     QTKT(7),TT(7),TAUT(4),TOUT(2,3),VK(6),T(6)
      REAL HDEP,HDES,HDE,HDE1,HDE2,HDE3,HDE4,HDE5,HDE6,
     *     HDE11,HDE12,HDE22,HDE13,HDE23,HDE33,HDE14,HDE24,HDE34,
     *     HDE44,HDE15,HDE25,HDE35,HDE45,HDE55,HDE16,HDE26,HDE36,
     *     HDE46,HDE56,HDE66,EE(9)
      INTEGER ISRF1T,IFUN1T(7),IFUN2T(7),ISRF1K,IFUN1K(6),IFUN2K(6),
     *     NSUM1,NSUM2,IX,NDER,NFUN11,NFUN12,NFUN21,NFUN22
      REAL AUX,PIPI,X1T,FUN1T(7),X1K,FUN1K(6)
      SAVE NSUM1,TT,IX,NDER,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T,NFUN21,
     *     IFUN2T,QTKT,
     *     NSUM2,T,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K,NFUN22,IFUN2K,VK
C     AA ...  Values, first and second partial derivatives of real
C             parts of 21 reduced (divided by the density) elastic
C             parameters.  The order of the value, first and second
C             partial derivatives of each parameter Aij is:
C             Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33.
C             The order of parameters (second array index) is:
C             A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45,
C             A55,A16,A26,A36,A46,A56,A66.
C     AAA...  3*3*3*3 matrix of reduced elastic parameters.
C     AAA1... 3*3*3*3 matrix of first partial derivatives of reduced
C             elastic parameters according to X1.
C     AAA2... 3*3*3*3 matrix of first partial derivatives of reduced
C             elastic parameters according to X2.
C     AAA3... 3*3*3*3 matrix of first partial derivatives of reduced
C             elastic parameters according to X3.
C-----------------------------------------------------------------------
C
      V0=RFVS(I1)
      VV0(1)=RFVS1(I1)
      VV0(2)=RFVS2(I1)
      VV0(3)=RFVS3(I1)
      P(1)=RFSLO1(I1)
      P(2)=RFSLO2(I1)
      P(3)=RFSLO3(I1)
      E1(1)=RFEVEC(1,I1)
      E1(2)=RFEVEC(2,I1)
      E1(3)=RFEVEC(3,I1)
      E2(1)=RFEVEC(4,I1)
      E2(2)=RFEVEC(5,I1)
      E2(3)=RFEVEC(6,I1)
      QQT(1)=RFQMAT(1,I1)
      QQT(2)=RFQMAT(2,I1)
      QQT(3)=RFQMAT(3,I1)
      QQT(4)=RFQMAT(4,I1)
      PPT(1)=RFPMAT(1,I1)
      PPT(2)=RFPMAT(2,I1)
      PPT(3)=RFPMAT(3,I1)
      PPT(4)=RFPMAT(4,I1)
      CALL PARM3(IFICB(I1),RFCOOR(1,I1),AA,RHO,QQ)
      CALL WAIJKL(AA,AAA,1)
      CALL WAIJKL(AA,AAA1,2)
      CALL WAIJKL(AA,AAA2,3)
      CALL WAIJKL(AA,AAA3,4)
      H1(1)=RFH(1,I1)
      H1(2)=RFH(2,I1)
      H1(3)=RFH(3,I1)
      H2(1)=RFH(4,I1)
      H2(2)=RFH(5,I1)
      H2(3)=RFH(6,I1)
      H3(1)=RFH(7,I1)
      H3(2)=RFH(8,I1)
      H3(3)=RFH(9,I1)
C
      IF (CRTANI.EQ.0.) THEN
        HH0(1)=VV0(1)/V0
        HH0(2)=VV0(2)/V0
        HH0(3)=VV0(3)/V0
        HHT0(1)=WASUM(HH0,H1)
        HHT0(2)=WASUM(HH0,H2)
        HHT0(3)=WASUM(HH0,H3)
        HH1=-1./SQRT(WASUM5(AAA,E1,P,E1,P))
        HH2=-1./SQRT(WASUM5(AAA,E2,P,E2,P))
        HH1L(1)=-0.5*(HH1**3)*WASUM5(AAA1,E1,P,E1,P)
        HH1L(2)=-0.5*(HH1**3)*WASUM5(AAA2,E1,P,E1,P)
        HH1L(3)=-0.5*(HH1**3)*WASUM5(AAA3,E1,P,E1,P)
        HH2L(1)=-0.5*(HH2**3)*WASUM5(AAA1,E2,P,E2,P)
        HH2L(2)=-0.5*(HH2**3)*WASUM5(AAA2,E2,P,E2,P)
        HH2L(3)=-0.5*(HH2**3)*WASUM5(AAA3,E2,P,E2,P)
        HH1U(1)=-(HH1**3)*WASUM4(AAA,E1,E1,P,1)
        HH1U(2)=-(HH1**3)*WASUM4(AAA,E1,E1,P,2)
        HH1U(3)=-(HH1**3)*WASUM4(AAA,E1,E1,P,3)
        HH2U(1)=-(HH2**3)*WASUM4(AAA,E2,E2,P,1)
        HH2U(2)=-(HH2**3)*WASUM4(AAA,E2,E2,P,2)
        HH2U(3)=-(HH2**3)*WASUM4(AAA,E2,E2,P,3)
        HHT1L(1)=WASUM(HH1L,H1)
        HHT1L(2)=WASUM(HH1L,H2)
        HHT2L(1)=WASUM(HH2L,H1)
        HHT2L(2)=WASUM(HH2L,H2)
        HHT1U(1)=WASUM(HH1U,H1)
        HHT1U(2)=WASUM(HH1U,H2)
        HHT2U(1)=WASUM(HH2U,H1)
        HHT2U(2)=WASUM(HH2U,H2)
        QTKT(1)=-((HHT1L(1)+HH1*HHT0(1))*QQT(1)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(2)+
     *            HHT1U(1)*PPT(1)+HHT1U(2)*PPT(2))
        QTKT(2)=-((HHT2L(1)+HH2*HHT0(1))*QQT(1)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(2)+
     *            HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2))
        QTKT(3)=-((HHT1L(1)+HH1*HHT0(1))*QQT(3)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(4)+
     *            HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4))
        QTKT(4)=-((HHT2L(1)+HH2*HHT0(1))*QQT(3)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(4)+
     *            HHT2U(1)*PPT(3)+HHT2U(2)*PPT(4))
        IF (IFIND(I1).EQ.1) THEN
C         Initial point of the ray:
C         Initiating first integration:
          NSUM1=4
          IX=1
          NDER=1
          NFUN11=4
          NFUN21=4
          IFUN1T(1)=1
          IFUN1T(2)=3
          IFUN1T(3)=2
          IFUN1T(4)=4
          FUN1T(1)=QTKT(1)
          FUN1T(2)=QTKT(2)
          FUN1T(3)=QTKT(3)
          FUN1T(4)=QTKT(4)
          IFUN2T(1)=1
          IFUN2T(2)=3
          IFUN2T(3)=2
          IFUN2T(4)=4
C         Initiating second integration:
          NSUM2=3
          NFUN12=3
          NFUN22=3
          IFUN1K(1)=1
          IFUN1K(2)=2
          IFUN1K(3)=3
          AUX=V0**2
          TAUT(1)=-HHT1U(1)/AUX
          TAUT(2)=-HHT1U(2)/AUX
          TAUT(3)=-HHT2U(1)/AUX
          TAUT(4)=-HHT2U(2)/AUX
          FUN1K(1)=-(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) +
     *               (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*AUX)
          FUN1K(2)=-(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) +
     *               (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*AUX)
          FUN1K(3)=-(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) +
     *               HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) +
     *               (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4))*AUX)
          IFUN2K(1)=1
          IFUN2K(2)=2
          IFUN2K(3)=3
C         Output values:
          DTAU(1)=0.
          DTAU(2)=0.
          DTAU(3)=0.
          DTAU(4)=0.
          DTAU(5)=0.
          DTAU(6)=0.
          DTAU(7)=0.
          RETURN
        ENDIF
        CALL AP28(NSUM1,TT,IX,NDER,STEP,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T,
     *            NFUN21,IFUN2T,QTKT)
        AUX=QQT(1)*QQT(4)-QQT(2)*QQT(3)
        IF (AUX.EQ.0) THEN
C         CRTPFA-22
          CALL ERROR('CRTPFA-22: Zero determinant of geom. spread.')
C         This error should not appear.
        ENDIF
        QIT(1)= QQT(4)/AUX
        QIT(4)= QQT(1)/AUX
        QIT(2)=-QQT(2)/AUX
        QIT(3)=-QQT(3)/AUX
        TAUT(1)=QIT(1)*TT(1)+QIT(2)*TT(2)
        TAUT(2)=QIT(3)*TT(1)+QIT(4)*TT(2)
        TAUT(3)=QIT(1)*TT(3)+QIT(2)*TT(4)
        TAUT(4)=QIT(3)*TT(3)+QIT(4)*TT(4)
        VK(1)=-(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) +
     *          (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*V0**2)
        VK(2)=-(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) +
     *          (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*V0**2)
        VK(3)=-(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) +
     *          HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) +
     *          (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4))*V0**2)
        CALL AP28(NSUM2,T,IX,NDER,STEP,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K,
     *            NFUN22,IFUN2K,VK)
C       Second order perturbations towards anisotropic rays:
        DTAU(1)=T(1)/2.
        DTAU(3)=T(2)/2.
C       Estimation of the two equal second order perturbations from
C       anisotropic common ray towards anisotropic rays:
        DTAU(2)=T(3)
        DTAU(4)=0.
        DTAU(5)=0.
        DTAU(6)=0.
        DTAU(7)=0.
        RETURN
      ENDIF
C
      IF (CRTANI.NE.0.) THEN
        HC1(1)=RFH(10,I1)
        HC1(2)=RFH(11,I1)
        HC1(3)=RFH(12,I1)
        HC2(1)=RFH(13,I1)
        HC2(2)=RFH(14,I1)
        HC2(3)=RFH(15,I1)
        HC3(1)=RFH(16,I1)
        HC3(2)=RFH(17,I1)
        HC3(3)=RFH(18,I1)
        CALL HDER(-1,AA,P(1),P(2),P(3),HDEP,HDES,HDE,
     *            HDE1,HDE2,HDE3,HDE4,HDE5,HDE6,
     *            HDE11,HDE12,HDE22,HDE13,HDE23,HDE33,HDE14,HDE24,HDE34,
     *            HDE44,HDE15,HDE25,HDE35,HDE45,HDE55,HDE16,HDE26,HDE36,
     *            HDE46,HDE56,HDE66,EE)
        HH0(1)=HDE1
        HH0(2)=HDE2
        HH0(3)=HDE3
        HHC0(1)=HDE4
        HHC0(2)=HDE5
        HHC0(3)=HDE6
        HHC0SP(1,1)=HDE44
        HHC0SP(2,1)=HDE45
        HHC0SP(3,1)=HDE46
        HHC0SP(1,2)=HDE45
        HHC0SP(2,2)=HDE55
        HHC0SP(3,2)=HDE56
        HHC0SP(1,3)=HDE46
        HHC0SP(2,3)=HDE56
        HHC0SP(3,3)=HDE66
C       Eq. 10:
        HHC0S(1,1)=HHC0SP(1,1)-3.*HHC0(1)*HHC0(1)
        HHC0S(2,1)=HHC0SP(2,1)-3.*HHC0(2)*HHC0(1)
        HHC0S(3,1)=HHC0SP(3,1)-3.*HHC0(3)*HHC0(1)
        HHC0S(1,2)=HHC0SP(1,2)-3.*HHC0(1)*HHC0(2)
        HHC0S(2,2)=HHC0SP(2,2)-3.*HHC0(2)*HHC0(2)
        HHC0S(3,2)=HHC0SP(3,2)-3.*HHC0(3)*HHC0(2)
        HHC0S(1,3)=HHC0SP(1,3)-3.*HHC0(1)*HHC0(3)
        HHC0S(2,3)=HHC0SP(2,3)-3.*HHC0(2)*HHC0(3)
        HHC0S(3,3)=HHC0SP(3,3)-3.*HHC0(3)*HHC0(3)
C       Eq. 39:
        HHT0(1)=WASUM(HH0,HC1)
        HHT0(2)=WASUM(HH0,HC2)
        HHT0(3)=WASUM(HH0,HC3)
C       Eq. 41:
        HHTC0S(1,1)=WASUM3(HHC0S,H1,H1)
        HHTC0S(1,2)=WASUM3(HHC0S,H1,H2)
        HHTC0S(2,1)=WASUM3(HHC0S,H2,H1)
        HHTC0S(2,2)=WASUM3(HHC0S,H2,H2)
C       Eq. 11:
        HH1=-1./SQRT(WASUM5(AAA,E1,P,E1,P))
        HH2=-1./SQRT(WASUM5(AAA,E2,P,E2,P))
C       Eq. 15:
        PIPI=(P(1)**2+P(2)**2+P(3)**2)
        HH3=-1./V0/SQRT(PIPI)
C       Eq. 12:
        HH1L(1)=-0.5*(HH1**3)*WASUM5(AAA1,E1,P,E1,P)
        HH1L(2)=-0.5*(HH1**3)*WASUM5(AAA2,E1,P,E1,P)
        HH1L(3)=-0.5*(HH1**3)*WASUM5(AAA3,E1,P,E1,P)
        HH2L(1)=-0.5*(HH2**3)*WASUM5(AAA1,E2,P,E2,P)
        HH2L(2)=-0.5*(HH2**3)*WASUM5(AAA2,E2,P,E2,P)
        HH2L(3)=-0.5*(HH2**3)*WASUM5(AAA3,E2,P,E2,P)
C       Eq. 16:
        HH3L(1)=VV0(1)/V0**2/SQRT(PIPI)
        HH3L(2)=VV0(2)/V0**2/SQRT(PIPI)
        HH3L(3)=VV0(3)/V0**2/SQRT(PIPI)
C       Eq. 13:
        HH1U(1)=-(HH1**3)*WASUM4(AAA,E1,E1,P,1)
        HH1U(2)=-(HH1**3)*WASUM4(AAA,E1,E1,P,2)
        HH1U(3)=-(HH1**3)*WASUM4(AAA,E1,E1,P,3)
        HH2U(1)=-(HH2**3)*WASUM4(AAA,E2,E2,P,1)
        HH2U(2)=-(HH2**3)*WASUM4(AAA,E2,E2,P,2)
        HH2U(3)=-(HH2**3)*WASUM4(AAA,E2,E2,P,3)
C       Eq. 17:
        HH3U(1)=1./V0/SQRT(PIPI**3)*P(1)
        HH3U(2)=1./V0/SQRT(PIPI**3)*P(2)
        HH3U(3)=1./V0/SQRT(PIPI**3)*P(3)
C       Eq. 42:
        HHT1L(1)=WASUM(HH1L,HC1)
        HHT1L(2)=WASUM(HH1L,HC2)
        HHT2L(1)=WASUM(HH2L,HC1)
        HHT2L(2)=WASUM(HH2L,HC2)
        HHT3L(1)=WASUM(HH3L,HC1)
        HHT3L(2)=WASUM(HH3L,HC2)
C       Eq. 43:
        HHT1U(1)=WASUM(HH1U,H1)
        HHT1U(2)=WASUM(HH1U,H2)
        HHT2U(1)=WASUM(HH2U,H1)
        HHT2U(2)=WASUM(HH2U,H2)
        HHT3U(1)=WASUM(HH3U,H1)
        HHT3U(2)=WASUM(HH3U,H2)
C       Eq. 55:
        QTKT(1)=-((HHT1L(1)+HH1*HHT0(1))*QQT(1)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(2)+
     *            HHT1U(1)*PPT(1)+HHT1U(2)*PPT(2))
        QTKT(2)=-((HHT2L(1)+HH2*HHT0(1))*QQT(1)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(2)+
     *            HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2))
        QTKT(3)=-((HHT3L(1)+HH3*HHT0(1))*QQT(1)+
     *            (HHT3L(2)+HH3*HHT0(2))*QQT(2)+
     *            HHT3U(1)*PPT(1)+HHT3U(2)*PPT(2))
        QTKT(4)=-((HHT1L(1)+HH1*HHT0(1))*QQT(3)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(4)+
     *            HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4))
        QTKT(5)=-((HHT2L(1)+HH2*HHT0(1))*QQT(3)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(4)+
     *            HHT2U(1)*PPT(3)+HHT2U(2)*PPT(4))
        QTKT(6)=-((HHT3L(1)+HH3*HHT0(1))*QQT(3)+
     *            (HHT3L(2)+HH3*HHT0(2))*QQT(4)+
     *            HHT3U(1)*PPT(3)+HHT3U(2)*PPT(4))
        QTKT(7)=-HH3
        IF (IFIND(I1).EQ.1) THEN
C         Initial point of the ray:
C         Initiating first integration:
          NSUM1=7
          IX=1
          NDER=1
          NFUN11=7
          NFUN21=7
          IFUN1T(1)=1
          IFUN1T(2)=3
          IFUN1T(3)=5
          IFUN1T(4)=2
          IFUN1T(5)=4
          IFUN1T(6)=6
          IFUN1T(7)=7
          FUN1T(1)=QTKT(1)
          FUN1T(2)=QTKT(2)
          FUN1T(3)=QTKT(3)
          FUN1T(4)=QTKT(4)
          FUN1T(5)=QTKT(5)
          FUN1T(6)=QTKT(6)
          FUN1T(7)=QTKT(7)
          IFUN2T(1)=1
          IFUN2T(2)=3
          IFUN2T(3)=5
          IFUN2T(4)=2
          IFUN2T(5)=4
          IFUN2T(6)=6
          IFUN2T(7)=7
C         Initiating second integration:
          NSUM2=6
          NFUN12=6
          NFUN22=6
          IFUN1K(1)=1
          IFUN1K(2)=2
          IFUN1K(3)=3
          IFUN1K(4)=4
          IFUN1K(5)=5
          IFUN1K(6)=6
          AUX=HHTC0S(1,1)*HHTC0S(2,2)-HHTC0S(2,1)*HHTC0S(1,2)
          IF (AUX.EQ.0) THEN
C           CRTPFA-23
            CALL ERROR
     *           ('CRTPFA-23: Zero determinant of second der. of Ham.')
C           This error should not appear.
          ENDIF
          GIJT(1,1)= HHTC0S(2,2)/AUX
          GIJT(2,2)= HHTC0S(1,1)/AUX
          GIJT(2,1)=-HHTC0S(2,1)/AUX
          GIJT(1,2)=-HHTC0S(1,2)/AUX
          TOUT(1,1)=-GIJT(1,1)*HHT1U(1)-GIJT(1,2)*HHT1U(2)
          TOUT(2,1)=-GIJT(2,1)*HHT1U(1)-GIJT(2,2)*HHT1U(2)
          TOUT(1,2)=-GIJT(1,1)*HHT2U(1)-GIJT(1,2)*HHT2U(2)
          TOUT(2,2)=-GIJT(2,1)*HHT2U(1)-GIJT(2,2)*HHT2U(2)
          TOUT(1,3)=-GIJT(1,1)*HHT3U(1)-GIJT(1,2)*HHT3U(2)
          TOUT(2,3)=-GIJT(2,1)*HHT3U(1)-GIJT(2,2)*HHT3U(2)
          FUN1K(1)=-(HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *               HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *    HHTC0S(1,1)*TOUT(1,1)*TOUT(1,1)+
     *    HHTC0S(2,1)*TOUT(2,1)*TOUT(1,1)+
     *    HHTC0S(1,2)*TOUT(1,1)*TOUT(2,1)+
     *    HHTC0S(2,2)*TOUT(2,1)*TOUT(2,1))
          FUN1K(2)=-(HHT1U(1)*TOUT(1,2)+HHT1U(2)*TOUT(2,2) +
     *               HHT2U(1)*TOUT(1,1)+HHT2U(2)*TOUT(2,1) +
     *    HHTC0S(1,1)*TOUT(1,1)*TOUT(1,2)+
     *    HHTC0S(2,1)*TOUT(2,1)*TOUT(1,2)+
     *    HHTC0S(1,2)*TOUT(1,1)*TOUT(2,2)+
     *    HHTC0S(2,2)*TOUT(2,1)*TOUT(2,2))
          FUN1K(3)=-(HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *               HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *    HHTC0S(1,1)*TOUT(1,2)*TOUT(1,2)+
     *    HHTC0S(2,1)*TOUT(2,2)*TOUT(1,2)+
     *    HHTC0S(1,2)*TOUT(1,2)*TOUT(2,2)+
     *    HHTC0S(2,2)*TOUT(2,2)*TOUT(2,2))
          FUN1K(4)=-(HHT1U(1)*TOUT(1,3)+HHT1U(2)*TOUT(2,3) +
     *               HHT3U(1)*TOUT(1,1)+HHT3U(2)*TOUT(2,1) +
     *    HHTC0S(1,1)*TOUT(1,1)*TOUT(1,3)+
     *    HHTC0S(2,1)*TOUT(2,1)*TOUT(1,3)+
     *    HHTC0S(1,2)*TOUT(1,1)*TOUT(2,3)+
     *    HHTC0S(2,2)*TOUT(2,1)*TOUT(2,3))
          FUN1K(5)=-(HHT2U(1)*TOUT(1,3)+HHT2U(2)*TOUT(2,3) +
     *               HHT3U(1)*TOUT(1,2)+HHT3U(2)*TOUT(2,2) +
     *    HHTC0S(1,1)*TOUT(1,2)*TOUT(1,3)+
     *    HHTC0S(2,1)*TOUT(2,2)*TOUT(1,3)+
     *    HHTC0S(1,2)*TOUT(1,2)*TOUT(2,3)+
     *    HHTC0S(2,2)*TOUT(2,2)*TOUT(2,3))
          FUN1K(6)=-(HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *               HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *    HHTC0S(1,1)*TOUT(1,3)*TOUT(1,3)+
     *    HHTC0S(2,1)*TOUT(2,3)*TOUT(1,3)+
     *    HHTC0S(1,2)*TOUT(1,3)*TOUT(2,3)+
     *    HHTC0S(2,2)*TOUT(2,3)*TOUT(2,3))
          IFUN2K(1)=1
          IFUN2K(2)=2
          IFUN2K(3)=3
          IFUN2K(4)=4
          IFUN2K(5)=5
          IFUN2K(6)=6
C         Output values:
          DTAU(1)=0.
          DTAU(2)=0.
          DTAU(3)=0.
          DTAU(4)=0.
          DTAU(5)=0.
          DTAU(6)=0.
          DTAU(7)=0.
          RETURN
        ENDIF
C       Eq. 51:
        CALL AP28(NSUM1,TT,IX,NDER,STEP,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T,
     *            NFUN21,IFUN2T,QTKT)
C       Eq. 57:
        AUX=QQT(1)*QQT(4)-QQT(2)*QQT(3)
        IF (AUX.EQ.0) THEN
C         CRTPFA-24
          CALL ERROR('CRTPFA-24: Zero determinant of geom. spread.')
C         This error should not appear.
        ENDIF
        QIT(1)= QQT(4)/AUX
        QIT(4)= QQT(1)/AUX
        QIT(2)=-QQT(2)/AUX
        QIT(3)=-QQT(3)/AUX
        TOUT(1,1)=QIT(1)*TT(1)+QIT(2)*TT(2)
        TOUT(2,1)=QIT(3)*TT(1)+QIT(4)*TT(2)
        TOUT(1,2)=QIT(1)*TT(3)+QIT(2)*TT(4)
        TOUT(2,2)=QIT(3)*TT(3)+QIT(4)*TT(4)
        TOUT(1,3)=QIT(1)*TT(5)+QIT(2)*TT(6)
        TOUT(2,3)=QIT(3)*TT(5)+QIT(4)*TT(6)
C       Eq. 62:
        VK(1)=-(HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *          HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *  HHTC0S(1,1)*TOUT(1,1)*TOUT(1,1)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,1)+
     *  HHTC0S(1,2)*TOUT(1,1)*TOUT(2,1)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,1))
        VK(2)=-(HHT1U(1)*TOUT(1,2)+HHT1U(2)*TOUT(2,2) +
     *          HHT2U(1)*TOUT(1,1)+HHT2U(2)*TOUT(2,1) +
     *  HHTC0S(1,1)*TOUT(1,1)*TOUT(1,2)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,2)+
     *  HHTC0S(1,2)*TOUT(1,1)*TOUT(2,2)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,2))
        VK(3)=-(HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *          HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *  HHTC0S(1,1)*TOUT(1,2)*TOUT(1,2)+HHTC0S(2,1)*TOUT(2,2)*TOUT(1,2)+
     *  HHTC0S(1,2)*TOUT(1,2)*TOUT(2,2)+HHTC0S(2,2)*TOUT(2,2)*TOUT(2,2))
        VK(4)=-(HHT1U(1)*TOUT(1,3)+HHT1U(2)*TOUT(2,3) +
     *          HHT3U(1)*TOUT(1,1)+HHT3U(2)*TOUT(2,1) +
     *  HHTC0S(1,1)*TOUT(1,1)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,3)+
     *  HHTC0S(1,2)*TOUT(1,1)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,3))
        VK(5)=-(HHT2U(1)*TOUT(1,3)+HHT2U(2)*TOUT(2,3) +
     *          HHT3U(1)*TOUT(1,2)+HHT3U(2)*TOUT(2,2) +
     *  HHTC0S(1,1)*TOUT(1,2)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,2)*TOUT(1,3)+
     *  HHTC0S(1,2)*TOUT(1,2)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,2)*TOUT(2,3))
        VK(6)=-(HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *          HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *  HHTC0S(1,1)*TOUT(1,3)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,3)*TOUT(1,3)+
     *  HHTC0S(1,2)*TOUT(1,3)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,3)*TOUT(2,3))
C       Eq. 59:
        CALL AP28(NSUM2,T,IX,NDER,STEP,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K,
     *            NFUN22,IFUN2K,VK)
C       Second order perturbations towards anisotropic rays:
        DTAU(1)=T(1)/2.
        DTAU(2)=T(2)
        DTAU(3)=T(3)/2.
C       Second order perturbation towards isotropic ray:
        DTAU(4)=T(4)
        DTAU(5)=T(5)
        DTAU(6)=T(6)/2.
C       Estimation of the isotropic ray travel time:
        DTAU(7)=TT(7)
        RETURN
      ENDIF
      END
C
C=======================================================================
C
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'an.for'
C     an.for
      INCLUDE 'apw.for'
C     apw.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'model.for'
C     model.for
      INCLUDE 'hder.for'
C     hder.for
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'wanpfa.for'
C     wanpfa.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'eigen.for'
C     eigen.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
C
C