C 
C Program WFSRF to generate wavefronts composed of polygons for display
C purposes
C
C Version: 5.40  !!! Preliminary version, poorly debugged            !!!
C                !!! Tested only within history file len-wf.h.       !!!
C                    len-wf.h
C Date: 2000, May 19
C                                                  
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C
C The post processing program to the CRT program performs linear
C interpolation of wavefronts. Wavefronts are approximated by triangles
C composed of points located at the rays of individual ray tubes. If
C such a triangle would intersect a structural interface, the wavefront
C is locally approximated by a triangle and a tetragon. Similarly in
C more complicated situations.
C
C This version of the program generates a single wavefront specified by
C travel time WFTIME during one invocation of WFSRF.
C
C The interpolation is not performed in demarcation belts between
C different ray histories, similarly as in program
C mtt.for.
C The maximum width of the belts is controlled
C by input parameter AERR of the CRT
C program.  The typical width of the belts, measured in take-off angles,
C is 0.0001 radians in order of magnitude.  The belts may become wide
C in areas of large geometrical spreading.  The division of rays into
C different histories is, by default, very detailed and is controlled
C by input parameter PRM0(3) of the
C CRT program.  Such a behaviour is useful for two-point ray tracing but
C has some awkward consequences on the interpolation.  For example, the
C demarcation belt between the rays incident at different sides of the
C computational volume for ray tracing extends across the whole model,
C causing an artificial gap within the continuous wavefronts.
C If the ray history does not account for the termination of a ray at
C different sides of the computational volume, the gaps corresponding
C to the edges of the computational volume are removed, but the corners
C along the edges are not filled with the ray cells any more, and
C the wavefront may disappear at the corners.  Then the
C computational volume for ray tracing should sufficiently exceed the
C extent of target volume.
C Input parameter PRM0(3) of the CRT
C program may thus considerably influence the results of the
C interpolation. However, gaps due to the demarcation belts
C corresponding to structural interfaces cannot be removed within
C the present interpolation algorithm.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of input files related to ray tracing (output of CRT):
C     CRTOUT='string'... File with the names of the output files of
C             program CRT.  A single set of CRT output files is read
C             from file CRTOUT.  Only the files 'CRT-R', 'CRT-I' and
C             'CRT-T' are read by WFSRF. File 'CRT-R' must contain all
C             rays traced by CRT.
C             For general description of file CRTOUT refer to
C             file writ.for.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C             'CRT-I'='r01i.out', 'CRT-T'='t01.out'
C Wavefront travel time:
C     WFTIME=real ... Travel time corresponding to the wavefront to
C             be covered by polygons.
C             Default: WFTIME=1
C Names of the output files:
C     VRTX='string'... Name of the output file with vertices of the
C             polygons.
C             Form of file VRTX
C             Default: VRTX='vrtx.out'
C     PLGN='string'... Name of the output file describing the polygons.
C             If blank, the file is not generated.
C             Form of file PLGN
C             Default: PLGN='plgn.out'
C     PLGNS='string'... Name of the file describing the polygons in
C             terms of the names of the vertices.
C             Form of file PLGNS
C             Default: PLGNS=' ' (the file is not generated)
C Parameters specifying the quantities to be written into the file VRTX.
C     TEXTP='string' ... First part of names of vertices. The second
C             part of the name of a vertex is formed by number giving
C             its position in the file VRTX.
C             Default: TEXTP='V'
C     COLUMN01 to COLUMN69, POWER01 to POWER69, IVALUE01 to IVALUE69:
C     IVALUEii=integer ... An integer value required for some special
C             values of COLUMNii. Reserved for future extensions.
C     POWERii=real ... Power of the quantity to be written in column ii.
C     COLUMNii='string' ... String which specifies the quantity to be
C             written to the column ii of the file VRTX. First six
C             columns usually contain coordinates of the vertices and
C             the normals. Column zero is reserved for names of the
C             vertices. Following strings are allowed:
C               ' ' (a space) ... Nothing is to be written to the column
C                        ii and to all the following columns.
C               'X1' ... First coordinate of the vertex.
C               'X2' ... Second coordinate of the vertex.
C               'X3' ... Third coordinate of the vertex.
C               'NORM1' ... First component of the normal to the
C                 wavefront at the vertex.  Slowness vectors are used
C                 as the normals in the current version of WFSRF.
C               'NORM2' ... Second component of the normal.
C               'NORM3' ... Third component of the normal.
C               'ICB' ... Index of complex block in which the vertex
C                 is located.
C             All strings may be entered either in uppercase or in
C             lowercase letters.
C             Defaults: COLUMN01='X1', COLUMN02='X2', COLUMN03='X3',
C               COLUMN04='NORM1', COLUMN05='NORM2', COLUMN06='NORM3',
C               COLUMN07 to COLUMN69=' ',
C               POWER01 to POWER69=1, IVALUE01 to IVALUE69=1.
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL WFEROR,WFREAD,WFERAS,WFTUBE,WFINTP,
     *         WFQD,WFQRI,WFQRA,WFSIDE,WFLINE,WFWRPT,
     *ERROR,WARN,RSEP1,RSEP3T,RSEP3R,RSEP3I,INDEXX,AP00,
     *FORM1,LOWER,LENGTH
      INTEGER LENGTH
C     WFEROR,WFREAD,WFERAS,WFTUBE,WFINTP,
C     WFQD,WFQRI,WFQRA,WFSIDE,WFLINE,WFWRPT  ... This file.
C     ERROR,WARN ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C               File sep.for.
C     INDEXX ... File indexx.for.
C     AP00 ... File ap.for.
C     FORM1,LOWER ... forms.for.
C     LENGTH ... length.for.
C
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C wfsrf.inc.
C.......................................................................
C     Auxiliary storage locations:
      INTEGER KRAMP
      PARAMETER (KRAMP=2)
      INTEGER IRAY1,IRAY2,IRAY3
      INTEGER IT1,IIT1
      INTEGER I
      INTEGER LU1,LU2,LU3,LU4,LU5
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5)
      CHARACTER*80 FILSEP,FILCRT,FILE1,FILE2,FILE3,CH
      INTEGER I1,I2,I3,I4,J1,J2,NPOINT,LEN1,LEN2,LENG
      INTEGER MLJPOL
      PARAMETER (MLJPOL=10)
      REAL Z(69),OUTMIN,OUTMAX
      CHARACTER*20 FORMAT,FORMA1,FORMA2,TEXTS(MLJPOL)
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      FILSEP=' '
      WRITE(*,'(A)') '+WFSRF: Enter input filename: '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       WFSRF-01
        CALL ERROR('WFSRF-01: No input file specified')
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
      WRITE(*,'(A)') '+WFSRF: Reading input data.   '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading filenames of the files with computed rays
C     and triangles:
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='t01.out'
      IF (FILCRT.NE.' ') THEN
        OPEN (LU2,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU2,*) FILE1,CH,FILE2,FILE3
        CLOSE (LU2)
      ENDIF
C
C     Opening the CRT output files:
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD')
C
C     Reading travel time(s) of the wavefront(s) to be triangulated,
C     filename(s) of the output file(s),
C     recording which quantities are to be written into the file(s):
      CALL WFQD
C
C     Reading file with computed triangles,
C     sorting the rays in triangles:
      WRITE(*,'(A)') '+WFSRF: Reading the file with triangles.'
      NT=0
      NRAMP=0
      IRAYMI=999999
      IRAYMA=0
  10  CONTINUE
        READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
        IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1
        IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2
        IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3
        IF (IRAY1.LT.IRAYMI)                    IRAYMI=IRAY1
        IF (IRAY2.LT.IRAYMI)                    IRAYMI=IRAY2
        IF ((IRAY3.NE.0).AND.(IRAY3.LT.IRAYMI)) IRAYMI=IRAY3
        IF (IRAY1.LT.IRAY2) THEN
          I=IRAY1
          IRAY1=IRAY2
          IRAY2=I
        ENDIF
        IF (IRAY2.LT.IRAY3) THEN
          I=IRAY2
          IRAY2=IRAY3
          IRAY3=I
        ENDIF
        IF (IRAY1.LT.IRAY2) THEN
          I=IRAY1
          IRAY1=IRAY2
          IRAY2=I
        ENDIF
        IF ((IRAY1.LE.0).OR.(IRAY2.LE.0).OR.(IRAY3.LT.0)) THEN
C         WFSRF-02
          CALL ERROR('WFSRF-02: Disorder in file CRT-T.')
C         The 'triangles' in the input file CRT-T should be written
C         as three nonnegative integers (indices of rays), two of them
C         being positive (the third may equal zero in 2-D computations).
        ENDIF
        IF (NRAMP+3.GT.MRAM) CALL WFEROR(1)
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY1
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY2
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY3
        NT=NT+1
      GOTO 10
  20  CONTINUE
      CLOSE(LU3)
      NR=IRAYMA-IRAYMI+1
C
C     Check for the 'triangles' in case of 2D computation:
      IF (IRAM(3).EQ.0) THEN
        DO 25, I1=6,3*NT,3
          IF (IRAM(I1).NE.0) THEN
C           WFSRF-03
            CALL ERROR('WFSRF-03: Disorder in triangles.')
C           The first 'triangle' is formed by two rays with zero
C           instead of the third ray. This identifies 2-D calculation
C           and all the 'triangles' should be formed by two rays and
C           zero instead of the third ray. This is not the case of the
C           current triangle. Either input file CRT-T
C           is wrong, or there is a bug in the code. In the latter case,
C           please, contact the author.
          ENDIF
  25    CONTINUE
        L3D=.FALSE.
      ELSE
        L3D=.TRUE.
      ENDIF
C
C     Sorting the triangles:
      WRITE(*,'(A)') '+WFSRF: Sorting the triangles.          '
      IF (NRAMP+2*NT.GT.MRAM) CALL WFEROR(1)
      DO 30, I1=NRAMP+NT+1,NRAMP+NT+NT
        RAM(I1)=FLOAT(IRAM((I1-NRAMP-NT-1)*3+1))
  30  CONTINUE
      CALL INDEXX(NT,RAM(NRAMP+NT+1),IRAM(NRAMP+1))
      NRAMP=NRAMP+NT
C
C     Setting the value of MRAMP:
      MRAMP=MRAM-KRAMP*(NWFT*NT*5 + NWFT*NR*NQPOI)
      IF (MRAMP.LE.NRAMP) CALL WFEROR(1)
      NPOI=MRAMP
      NPLGN=MRAM+1
C
C
C     Forming an auxiliary array with information about when can be rays
C     erased from the memory ("deleting array"):
      IF (NRAMP+NR.GT.MRAMP) CALL WFEROR(1)
      DO 50, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  50  CONTINUE
      NRAMP=NRAMP+NR
      ORAYE=IRAYMI-1-4*NT
      DO 60, I2=3*NT+1,4*NT
        I1=(IRAM(I2)-1)*3+1
        IRAM(IRAM(I1  )-ORAYE)=IRAM(I1)
        IRAM(IRAM(I1+1)-ORAYE)=IRAM(I1)
        IF (IRAM(I1+2).NE.0) IRAM(IRAM(I1+2)-ORAYE)=IRAM(I1)
  60  CONTINUE
C
C
C     Forming an auxiliary array with information about addresses
C     of the ends of records for rays in array RAM ("addressing array"):
C     "Ray" IRAYMI-1:
      NRAMP=NRAMP+1
      IF (NRAMP.GT.MRAMP) CALL WFEROR(1)
      IRAM(NRAMP)=NRAMP+NR+NWFT*NR
C     All other rays:
      IF (NRAMP+NR.GT.MRAMP) CALL WFEROR(1)
      DO 70, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  70  CONTINUE
      NRAMP=NRAMP+NR
      ORAYA=IRAYMI-1-4*NT-NR-1
C
C
C     Forming an auxiliary array(s) with information about addresses
C     of the ends of records for points on wavefronts according to rays
C     in array RAM ("points at wavefronts addressing array(s)"):
      IF (NRAMP+NWFT*NR.GT.MRAMP) CALL WFEROR(1)
      DO 80, I1=NRAMP+1,NRAMP+NWFT*NR
        IRAM(I1)=0
  80  CONTINUE
      NRAMP=NRAMP+NWFT*NR
      OWFT(1)=IRAYMI-1-4*NT-NR-NR-1
C     OWFT(2)=IRAYMI-1-4*NT-NR-NR-1-NR
C     OWFT(3)=IRAYMI-1-4*NT-NR-NR-1-NR-NR
C
C
C     Loop for all the triangles:
      IIT1=1
      I1=-1
  100 CONTINUE
        IT1=(IRAM(3*NT+IIT1)-1)*3+1
        I2=NINT((100.*IIT1)/(NT))
        IF (I2.NE.I1) THEN
          WRITE(*,'(''+'',79('' ''))')
          WRITE(*,'(A,1I4,A)')
     *    '+WFSRF: Interpolating. ',I2,'% of ray tubes completed.'
          I1=I2
        ENDIF
C
C       If necessary, reading new rays:
        IF
     *  (((IRAM(IRAM(IT1)-ORAYA+1).EQ.0).AND.(IRAM(IT1).NE.IRAYMA)) .OR.
     *   ((IRAM(IRAM(IT1)-ORAYA  ).EQ.0).AND.(IRAM(IT1).EQ.IRAYMA)))
     *  CALL WFREAD(LU1,LU2,IT1)
C
C       Empting the array RAM to enable writing of new points
C       and polygons at wavefronts:
        IF ((MRAMP-NRAMP).LT.(NINT(MRAM/10.))) CALL WFERAS
C
C       Creating polygons at wavefronts along the rays of the triangle:
        CALL WFTUBE(IT1)
C
      IIT1=IIT1+1
      IF (IIT1.LE.NT) GOTO 100
C     End of the loop for all the triangles.
      CLOSE(LU1)
      CLOSE(LU2)
C
C
C     Storing points and polygons on wavefronts:
      WRITE(*,'(A)')
     *'+WFSRF: Writing.                                     '
C
C     Storing points:
      OPEN(LU1,FILE=VRTX,FORM='FORMATTED')
      WRITE(LU1,'(A)') '/'
      NPOINT=(NPOI-MRAMP)/NQPOI
C     Length of the names of points:
      LEN1=LENGTH(TEXTP)
      LEN2=0
 202  CONTINUE
        IF (NPOINT.GE.10**LEN2) THEN
          LEN2=LEN2+1
          GOTO 202
        ENDIF
      LENG=LEN1+LEN2
C     Loop over points:
      DO 218, I1=1,NPOINT
C       Address in RAM just before the current point:
        J1=MRAMP+(I1-1)*NQPOI
C       Name of the point:
        DO 204, I2=0,LEN2-1
          TEXTP(LENG-I2:LENG-I2)=
     *    CHAR(ICHAR('0')+MOD(I1,10**(I2+1))/10**I2)
 204    CONTINUE
C       The quantities at the point:
        J2=0
C       Loop over the quantities to be stored:
 206    CONTINUE
          IF ((TEXTC(J2+1).NE.' ').AND.(J2+1.LE.69)) THEN
            J2=J2+1
            IF     (TEXTC(J2).EQ.'x1') THEN
C             First coordinate of the point:
              IF (POWER(J2).EQ.1.) THEN
                Z(J2)=RAM(J1+2)
              ELSE
                Z(J2)=RAM(J1+2)**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'x2') THEN
C             Second coordinate of the point:
              IF (POWER(J2).EQ.1.) THEN
                Z(J2)=RAM(J1+3)
              ELSE
                Z(J2)=RAM(J1+3)**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'x3') THEN
C             Third coordinate of the point:
              IF (POWER(J2).EQ.1.) THEN
                Z(J2)=RAM(J1+4)
              ELSE
                Z(J2)=RAM(J1+4)**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'norm1') THEN
C             First component of the normal:
              IF (POWER(J2).EQ.1.) THEN
                Z(J2)=RAM(J1+5)
              ELSE
                Z(J2)=RAM(J1+5)**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'norm2') THEN
C             Second component of the normal:
              IF (POWER(J2).EQ.1.) THEN
                Z(J2)=RAM(J1+6)
              ELSE
                Z(J2)=RAM(J1+6)**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'norm3') THEN
C             Third component of the normal:
              IF (POWER(J2).EQ.1.) THEN
                Z(J2)=RAM(J1+7)
              ELSE
                Z(J2)=RAM(J1+7)**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'icb') THEN
              Z(J2)=FLOAT(IRAM(J1+1))**POWER(J2)
            ENDIF
            GOTO 206
          ENDIF
C       End of the loop for quantities to be stored.
C       Writing the quantities at the point:
C       Setting output format for the array:
        FORMAT='(3A,00(F00.0,1X),A)'
        FORMAT(6:6)=CHAR(ICHAR('0')+MOD(J2,10))
        FORMAT(5:5)=CHAR(ICHAR('0')+J2/10)
        OUTMIN=0.
        OUTMAX=0.
        DO 214, I2=1,J2
          IF (OUTMIN.GT.Z(I2)) OUTMIN=Z(I2)
          IF (OUTMAX.LT.Z(I2)) OUTMAX=Z(I2)
 214    CONTINUE
        CALL FORM1(OUTMIN,OUTMAX,FORMAT(8:15))
        FORMAT(14:17)= '1X),'
C       Output format is set.
        WRITE(LU1,FORMAT) '''',TEXTP(1:LENG),''' ',(Z(I2),I2=1,J2),'/'
 218  CONTINUE
C     End of the loop over all points.
      WRITE(LU1,'(A)') '/'
      CLOSE(LU1)
C
C     Storing polygons:
      IF(PLGN.NE.' ') OPEN(LU1,FILE=PLGN,FORM='FORMATTED')
      IF(PLGNS.NE.' ')OPEN(LU2,FILE=PLGNS,FORM='FORMATTED')
      FORMA1='(00(I8,1X),A)'
      FORMA2='(00(3A,1X),A)'
      I1=MRAM-1
C     Loop over polygons:
 220  CONTINUE
      IF (I1.GT.NPLGN) THEN
        J1=IRAM(I1)
        I3=MOD(J1,10)
        FORMA1(3:3)=CHAR(ICHAR('0')+I3)
        FORMA2(3:3)=CHAR(ICHAR('0')+I3)
        I3=J1/10
        FORMA1(2:2)=CHAR(ICHAR('0')+I3)
        FORMA2(2:2)=CHAR(ICHAR('0')+I3)
        J2=IRAM(I1-1)
        IF (PLGN.NE.' ') THEN
          WRITE(LU1,FORMA1)
     *         ((IRAM(I2)-MRAMP)/NQPOI,I2=I1-1,I1-J1,-1),'/'
        ENDIF
        IF (PLGNS.NE.' ') THEN
          IF (J1.GT.MLJPOL) THEN
C           WFSRF-04
            CALL ERROR ('WFSRF-04: Small array for names of polygons.')
C           This error may be caused by small dimension MLJPOL of
C           array TEXTS defined in this program.
          ENDIF
          DO 224, I2=I1-1,I1-J1,-1
C           Index of the point (from 1 to NPOINT):
            J2=(IRAM(I2)-MRAMP)/NQPOI
C           Name of the point:
            I4=I1-I2
            DO 222, I3=0,LEN2-1
              TEXTS(I4)(1:LEN1)=TEXTP(1:LEN1)
              TEXTS(I4)(LENG-I3:LENG-I3)=
     *        CHAR(ICHAR('0')+MOD(J2,10**(I3+1))/10**I3)
 222        CONTINUE
 224      CONTINUE
          WRITE(LU2,FORMA2)
     *         ('''',TEXTS(I2)(1:LENG),'''',I2=1,J1),'/'
        ENDIF
        I1=I1-(J1+2)
        GOTO 220
      ENDIF
C     End of the loop over polygons.
      IF (PLGN.NE.' ') THEN
        WRITE(LU1,'(A)') '/'
        CLOSE(LU1)
      ENDIF
      IF (PLGNS.NE.' ') THEN
        WRITE(LU2,'(A)') '/'
        CLOSE(LU2)
      ENDIF
C
      WRITE(*,'(A)')
     *'+WFSRF: Done.                                        '
      STOP
      END
C
C
C=======================================================================
C
      SUBROUTINE WFREAD(LU1,LU2,IT1)
C
C----------------------------------------------------------------------
C Subroutine to read the unformatted output of program CRT and
C to write it into array (I)RAM used in program WFSRF.
C Reading the output files is completed by a simple invocation of
C subroutine AP00 of file 'ap.for'.
C
      INTEGER LU1,LU2,IT1
C Input:
C     LU1 ... Number of logical unit corresponding to the file with
C             the quantities stored along rays (see C.R.T.5.5.1).
C     LU2 ... Number of logical unit corresponding to the file with
C             the quantities at the initial points of rays,
C             corresponding to file 'RAYPOINTS' (see C.R.T.6.1).
C     IT1 ... Position of the first ray of the actually processed
C             triangle in the array IRAM.
C No output.
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C None of the storage locations of the common block are altered.
C ...........................
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C wfsrf.inc.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY0,I1
C     IRAY0.. Index of the last ray read in into the array RAM.
C
C-----------------------------------------------------------------------
C     Loop for the points of rays:
   10 CONTINUE
        IF ((NRAMP+2*NQ).GT.MRAMP) THEN
C         Empting the memory:
          CALL WFERAS
          IF ((NRAMP+2*NQ).GT.MRAMP) CALL WFEROR(1)
        ENDIF
C       Reading the results of the complete ray tracing:
        IRAY0=IRAY
        CALL AP00(0,LU1,LU2)
        IF (IWAVE.LT.1) THEN
C         End of rays:
          IF (IRAY0.LT.IRAYMA) THEN
C           Some rays at the end of the CRT output file missing:
            DO 12, I1=IRAY0+1,IRAYMA
              IF (I1.GE.IRAYMI) THEN
                IRAM(I1-ORAYA)=NRAMP
              ENDIF
  12        CONTINUE
          ENDIF
          RETURN
        ENDIF
        IF (IRAY.LT.IRAYMI) GOTO 10
        IF (IRAY.GT.IRAYMA) GOTO 10
        IF ((IRAY-IRAY0).GE.2) THEN
C         Some rays skipped by AP00:
          DO 15, I1=IRAY0+1,IRAY-1
            IF (I1.GE.IRAYMI) THEN
              IRAM(I1-ORAYA)=IRAM(IRAY0-ORAYA)
            ENDIF
  15      CONTINUE
        ENDIF
        IF (IRAM(IRAY-ORAYE).NE.0) THEN
C         Writing the results of the complete ray tracing - recording
C         new point on a ray:
          IF (IPT.LE.1) THEN
            IF ((ISHEET.EQ.0).AND.(IRAY.EQ.IRAYMI)) THEN
C             WFSRF-05
              CALL WARN
     *        ('WFSRF-05: a ray with history equal to 0 was observed')
C              All the rays are probably of the same history 0. Only the
C             initial-value ray tracing was performed.  Width and shape
C             of ray tubes were not controlled.  The interpolation
C             in such a case makes sense only in smooth models.
C              Check the value of the parameter IPOINT in CRT input data
C             RPAR (4).
            ENDIF
C           New ray - recording the initial point:
C           Coordinates:
            RAM(NRAMP+1)=YI(3)
            RAM(NRAMP+2)=YI(4)
            RAM(NRAMP+3)=YI(5)
C           Index of surface:
            IRAM(NRAMP+4)=0
C           Index of complex block:
            IRAM(NRAMP+5)=ICB1I
C           Travel time, slowness vector,
C           and quantities to be interpolated:
            CALL WFQRI
          ENDIF
C         Recording the new point:
C         Coordinates:
          RAM(NRAMP+1)=Y(3)
          RAM(NRAMP+2)=Y(4)
          RAM(NRAMP+3)=Y(5)
C         Index of surface:
          IRAM(NRAMP+4)=ISRF
C         Index of complex block:
          IRAM(NRAMP+5)=ICB1
C         Travel time, slowness vector,
C         and quantities to be interpolated:
          CALL WFQRA
        ENDIF
        IRAM(IRAY-ORAYA)=NRAMP
        IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA))
     *     RETURN
      GOTO 10
C
      END
C
C
C=======================================================================
C
      SUBROUTINE WFERAS
C
C----------------------------------------------------------------------
C Subroutine for empting the array (I)RAM. All the parameters
C of all the rays, which will no more be used, are erased.
C
C No input.
C No output.
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     IRAY .. Index of the ray being actually read in by CIREAD.
C             This procedure supposes, that any ray with higher
C             index than IRAY was not read in.
C None of the storage locations of the common block are altered.
C ...........................
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C wfsrf.inc.
C.......................................................................
C     Auxiliary storage locations:
      INTEGER I1,I2,J1
      INTEGER IADDRP
C     I1 ... Controls main loop over rays.
C     I2 ... Controls the loop over parameters of ray I1.
C     J1 ... address of the last used record of array RAM.
C
C-----------------------------------------------------------------------
      J1=IRAM(IRAYMI-1-ORAYA)
      IADDRP=J1
C     Loop for the rays:
      DO 20, I1=IRAYMI,IRAY
        IF (IRAM(I1-ORAYE).GE.(IRAY-1)) THEN
C         This ray is not to be erased:
          DO 10, I2=IADDRP+1,IRAM(I1-ORAYA)
            J1=J1+1
            IRAM(J1)=IRAM(I2)
   10     CONTINUE
        ELSE
C         This ray is to be erased:
          IRAM(I1-ORAYE)=0
        ENDIF
        IADDRP=IRAM(I1-ORAYA)
        IRAM(I1-ORAYA)=J1
   20 CONTINUE
      NRAMP=J1
      RETURN
      END
C=======================================================================
C
      SUBROUTINE WFTUBE(IT1)
C
C----------------------------------------------------------------------
C Subroutine for interpolation within ray tube formed by the rays
C IRAM(IT1), IRAM(IT1+1), IRAM(IT1+2).
C
      INTEGER IT1
C Input:
C     IT1 ... The address of the index of the first ray of the ray tube,
C             in which the interpolation is to be performed.
C No output.
C
C ...........................
C ...........................
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C wfsrf.inc.
C.......................................................................
      INTEGER IRAY1,IRAY2,IRAY3
      INTEGER I1MI,I2MI,I3MI,I1MA,I2MA,I3MA,I1,I2,J1,J2
      INTEGER KRAY1,KRAY2,KRAY3,ICBR1,ICBR2,ICBR3,NEXT
      INTEGER NPOL,MPOL
      PARAMETER (MPOL=10)
      INTEGER KPOL(MPOL)
      REAL X1,X2,X3,P1,P2,P3
      LOGICAL LSHIFT
C     J1 ...
C     J2 ...
C
C-----------------------------------------------------------------------
      IRAY1=IRAM(IT1  )
      IRAY2=IRAM(IT1+1)
      IRAY3=IRAM(IT1+2)
      IF (         (IRAM(IRAY1-ORAYA).EQ.0).OR.
     *             (IRAM(IRAY2-ORAYA).EQ.0).OR.
     *    (L3D.AND.(IRAM(IRAY3-ORAYA).EQ.0))) THEN
C       WFSRF-06
        CALL ERROR
     *     ('WFSRF-06: Parameters of a ray not found in WFTUBE.')
C       This error may be caused by
C       K2P
C       not equal to zero, then only two-point rays are stored in
C       output files of CRT. We recommend to run CRT with file
C       writall.dat.
      ENDIF
C
      I1MI=         IRAM(IRAY1-1-ORAYA)
      I2MI=         IRAM(IRAY2-1-ORAYA)
      IF (L3D) I3MI=IRAM(IRAY3-1-ORAYA)
      I1MA=         IRAM(IRAY1-ORAYA)-NQ
      I2MA=         IRAM(IRAY2-ORAYA)-NQ
      IF (L3D) I3MA=IRAM(IRAY3-ORAYA)-NQ
      IF ((I1MI.GT.I1MA).OR.(I2MI.GT.I2MA).OR.(L3D.AND.(I3MI.GT.I3MA)))
     *  RETURN
C
C     Loop for wavefronts, searching for a polygon
C     which is an intersection of the tube with the wavefront:
      DO 100, I1=1,NWFT
C
        KRAY1=0
        KRAY2=0
        KRAY3=0
        ICBR1=0
        ICBR2=0
        ICBR3=0
C
C       Computing point on the first ray:
        IF (IRAM(IRAY1-OWFT(I1)).EQ.0) THEN
C         The point is not computed, computing the point:
          IF ((RAM(I1MI+6).LE.WFTIME(I1)).AND.
     *        (RAM(I1MA+6).GE.WFTIME(I1))) THEN
C           The point exists, computing the point:
            DO 10, I2=I1MI+NQ,I1MA,NQ
              IF (RAM(I2+6).GE.WFTIME(I1)) THEN
                CALL WFINTP(I2-NQ,I2,WFTIME(1),
     *                      X1,X2,X3,P1,P2,P3)
                CALL WFWRPT(IRAM(I2+5),X1,X2,X3,P1,P2,P3)
                KRAY1=NPOI
                ICBR1=IRAM(KRAY1-NQPOI+1)
                IRAM(IRAY1-OWFT(I1))=NPOI
                GOTO 11
              ENDIF
  10        CONTINUE
  11        CONTINUE
          ELSE
C           The point does not exist:
            IRAM(IRAY1-OWFT(I1))=-1
          ENDIF
        ELSEIF (IRAM(IRAY1-OWFT(I1)).EQ.-1) THEN
C         The point does not exist:
          CONTINUE
        ELSE
C         The point is already computed:
          KRAY1=IRAM(IRAY1-OWFT(I1))
          ICBR1=IRAM(KRAY1-NQPOI+1)
        ENDIF
C
C       Computing point on the second ray:
        IF (IRAM(IRAY2-OWFT(I1)).EQ.0) THEN
C         The point is not computed, computing the point:
          IF ((RAM(I2MI+6).LE.WFTIME(I1)).AND.
     *        (RAM(I2MA+6).GE.WFTIME(I1))) THEN
C           The point exists, computing the point:
            DO 20, I2=I2MI+NQ,I2MA,NQ
              IF (RAM(I2+6).GE.WFTIME(I1)) THEN
                CALL WFINTP(I2-NQ,I2,WFTIME(1),
     *                      X1,X2,X3,P1,P2,P3)
                CALL WFWRPT(IRAM(I2+5),X1,X2,X3,P1,P2,P3)
                KRAY2=NPOI
                ICBR2=IRAM(KRAY2-NQPOI+1)
                IRAM(IRAY2-OWFT(I1))=NPOI
                GOTO 21
              ENDIF
  20        CONTINUE
  21        CONTINUE
          ELSE
C           The point does not exist:
            IRAM(IRAY2-OWFT(I1))=-1
          ENDIF
        ELSEIF (IRAM(IRAY2-OWFT(I1)).EQ.-1) THEN
C         The point does not exist:
          CONTINUE
        ELSE
C         The point is already computed:
          KRAY2=IRAM(IRAY2-OWFT(I1))
          ICBR2=IRAM(KRAY2-NQPOI+1)
        ENDIF
C
C       Computing point on the third ray:
        IF (L3D) THEN
          IF (IRAM(IRAY3-OWFT(I1)).EQ.0) THEN
C           The point is not computed, computing the point:
            IF ((RAM(I3MI+6).LE.WFTIME(I1)).AND.
     *          (RAM(I3MA+6).GE.WFTIME(I1))) THEN
C             The point exists, computing the point:
              DO 30, I2=I3MI+NQ,I3MA,NQ
                IF (RAM(I2+6).GE.WFTIME(I1)) THEN
                  CALL WFINTP(I2-NQ,I2,WFTIME(1),
     *                        X1,X2,X3,P1,P2,P3)
                  CALL WFWRPT(IRAM(I2+5),X1,X2,X3,P1,P2,P3)
                  KRAY3=NPOI
                  ICBR3=IRAM(KRAY3-NQPOI+1)
                  IRAM(IRAY3-OWFT(I1))=NPOI
                  GOTO 31
                ENDIF
  30          CONTINUE
  31          CONTINUE
            ELSE
C             The point does not exist:
              IRAM(IRAY3-OWFT(I1))=-1
            ENDIF
          ELSEIF (IRAM(IRAY3-OWFT(I1)).EQ.-1) THEN
C           The point does not exist:
            CONTINUE
          ELSE
C           The point is already computed:
            KRAY3=IRAM(IRAY3-OWFT(I1))
            ICBR3=IRAM(KRAY3-NQPOI+1)
          ENDIF
        ENDIF
C
C       Recording point on the first ray:
        IF (KRAY1.NE.0) THEN
          KPOL(1)=KRAY1
          NPOL=1
        ELSE
          NPOL=0
        ENDIF
C
C       Computing and recording points
C       between the first and the second ray:
        IF (ICBR1.NE.ICBR2) THEN
          NEXT=0
  40      CONTINUE
            CALL WFSIDE(IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2,
     *                  I1MI,I1MA,I2MI,I2MA,WFTIME(I1),NEXT)
            NPOL=NPOL+1
            IF (NPOL.GT.MPOL) CALL WFEROR(4)
            KPOL(NPOL)=NPOI
          IF (NEXT.NE.0) GOTO 40
        ENDIF
C
C       Recording point on the second ray:
        IF (KRAY2.NE.0) THEN
          NPOL=NPOL+1
          IF (NPOL.GT.MPOL) CALL WFEROR(4)
          KPOL(NPOL)=KRAY2
        ENDIF
C
        IF (L3D) THEN
C         Computing and recording points
C         between the second and the third ray:
          IF (ICBR2.NE.ICBR3) THEN
            NEXT=0
  50        CONTINUE
              CALL WFSIDE(IRAY2,KRAY2,ICBR2,IRAY3,KRAY3,ICBR3,
     *                    I2MI,I2MA,I3MI,I3MA,WFTIME(I1),NEXT)
              NPOL=NPOL+1
              IF (NPOL.GT.MPOL) CALL WFEROR(4)
              KPOL(NPOL)=NPOI
            IF (NEXT.NE.0) GOTO 50
          ENDIF
C
C         Recording point on the third ray:
          IF (KRAY3.NE.0) THEN
            NPOL=NPOL+1
            IF (NPOL.GT.MPOL) CALL WFEROR(4)
            KPOL(NPOL)=KRAY3
          ENDIF
C
C         Computing and recording points
C         between the third and the first ray:
          IF (ICBR3.NE.ICBR1) THEN
            NEXT=0
  60        CONTINUE
              CALL WFSIDE(IRAY3,KRAY3,ICBR3,IRAY1,KRAY1,ICBR1,
     *                    I3MI,I3MA,I1MI,I1MA,WFTIME(I1),NEXT)
              NPOL=NPOL+1
              IF (NPOL.GT.MPOL) CALL WFEROR(4)
              KPOL(NPOL)=NPOI
            IF (NEXT.NE.0) GOTO 60
          ENDIF
        ENDIF
C
C       The polygon is found, shifting it in case that it contains
C       points in more than one complex block:
        LSHIFT=.FALSE.
        IF (L3D) THEN
          DO 70, I2=2,NPOL
            IF (IRAM(KPOL(1)-NQPOI+1).NE.IRAM(KPOL(I2)-NQPOI+1))
     *        LSHIFT=.TRUE.
  70      CONTINUE
        ENDIF
        IF (LSHIFT) THEN
  72      CONTINUE
            IF (IRAM(KPOL(1)-NQPOI+1).EQ.IRAM(KPOL(NPOL)-NQPOI+1)) THEN
              IF (NPOL+1.GT.MPOL) CALL WFEROR(4)
              KPOL(NPOL+1)=KPOL(1)
              DO 74, I2=1,NPOL
                KPOL(I2)=KPOL(I2+1)
  74          CONTINUE
              GOTO 72
            ENDIF
        ENDIF
C
        IF (NPOL.GE.2) THEN
C         Recording the polygon to RAM:
          J2=1
C         Loop over polygons:
  80      CONTINUE
            NPLGN=NPLGN-1
            IF (NPLGN.LE.NPOI) CALL WFEROR(3)
            RAM(NPLGN)=WFTIME(I1)
            NPLGN=NPLGN-1
            IF (NPLGN.LE.NPOI) CALL WFEROR(3)
            J1=NPLGN
            IRAM(J1)=1
            NPLGN=NPLGN-1
            IF (NPLGN.LE.NPOI) CALL WFEROR(3)
            IRAM(NPLGN)=KPOL(J2)
C           Loop over points of the polygon:
  82        CONTINUE
              J2=J2+1
              IF (J2.GT.NPOL) GOTO 89
              IF (IRAM(KPOL(J2)-NQPOI+1).NE.IRAM(KPOL(J2-1)-NQPOI+1))
     *           GOTO 80
              IRAM(J1)=IRAM(J1)+1
              NPLGN=NPLGN-1
              IF (NPLGN.LE.NPOI) CALL WFEROR(3)
              IRAM(NPLGN)=KPOL(J2)
            GOTO 82
  89      CONTINUE
        ENDIF
C
 100  CONTINUE
C     End of the loop over wavefronts.
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WFQD
C
C.......................................................................
C
C Subroutine designed to read the input data describing the names
C of the output files with the interpolated quantities, to read the
C quantities computed by CRT, to interpolate and to write the output
C files with the interpolated quantities.
C
C.......................................................................
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C For the description of the parameters stored in the common block
C WFSRFC refer to the file wfsrf.inc.
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C pointc.inc.
C None of the storage locations of the common block are altered.
C.......................................................................
C     Auxiliary storage locations for all the entries:
      INTEGER I1
      CHARACTER*20 FORMA1,FORMA2,FORMA3,TEXT
C.......................................................................
C
C     Reading the wavefront travel time:
      CALL RSEP3R('WFTIME',WFTIME(1),1.)
      NWFT=1
C
C     Recalling the output filenames:
      CALL RSEP3T('VRTX',VRTX,'vrtx.out')
      CALL RSEP3T('PLGN',PLGN,'plgn.out')
      CALL RSEP3T('PLGNS',PLGNS,' ')
C
C     Recalling the first part of names of points in output file VRTX:
      CALL RSEP3T('TEXTP',TEXTP,'V')
C
C     Recalling the data specifying the quantities to be written into
C     the output file with points at structural interfaces:
      FORMA1='COLUMN00'
      FORMA2='POWER00'
      FORMA3='IVALUE00'
      I1=1
   5  CONTINUE
        FORMA1(8:8)=CHAR(ICHAR('0')+MOD(I1,10))
        FORMA2(7:7)=FORMA1(8:8)
        FORMA3(8:8)=FORMA1(8:8)
        FORMA1(7:7)=CHAR(ICHAR('0')+I1/10)
        FORMA2(6:6)=FORMA1(7:7)
        FORMA3(7:7)=FORMA1(7:7)
        IF     (I1.EQ.1) THEN
          CALL RSEP3T(FORMA1,TEXTC(I1),'X1')
        ELSEIF (I1.EQ.2) THEN
          CALL RSEP3T(FORMA1,TEXTC(I1),'X2')
        ELSEIF (I1.EQ.3) THEN
          CALL RSEP3T(FORMA1,TEXTC(I1),'X3')
        ELSEIF (I1.EQ.4) THEN
          CALL RSEP3T(FORMA1,TEXTC(I1),'NORM1')
        ELSEIF (I1.EQ.5) THEN
          CALL RSEP3T(FORMA1,TEXTC(I1),'NORM2')
        ELSEIF (I1.EQ.6) THEN
          CALL RSEP3T(FORMA1,TEXTC(I1),'NORM3')
        ELSE
          CALL RSEP3T(FORMA1,TEXTC(I1),' ')
        ENDIF
        CALL RSEP3R(FORMA2,POWER(I1),1.)
        CALL RSEP3I(FORMA3,IVALUE(I1),1)
        IF (TEXTC(I1).NE.' ') THEN
          CALL LOWER(TEXTC(I1))
          IF ((TEXTC(I1).NE.'x1').AND.(TEXTC(I1).NE.'x2').AND.
     *        (TEXTC(I1).NE.'x3').AND.(TEXTC(I1).NE.'norm1').AND.
     *        (TEXTC(I1).NE.'norm2').AND.(TEXTC(I1).NE.'norm3')
     *        .AND.(TEXTC(I1).NE.'icb')
C    *        (TEXTC(I1).NE.'isrf').AND.
C    *        (TEXTC(I1).NE.'+isb').AND.(TEXTC(I1).NE.'-isb').AND.
C    *        (TEXTC(I1).NE.'+icb').AND.(TEXTC(I1).NE.'-icb').AND.
C    *        (TEXTC(I1).NE.'+vp') .AND.(TEXTC(I1).NE.'-vp') .AND.
C    *        (TEXTC(I1).NE.'+vs') .AND.(TEXTC(I1).NE.'-vs') .AND.
C    *        (TEXTC(I1).NE.'+den').AND.(TEXTC(I1).NE.'-den').AND.
C    *        (TEXTC(I1).NE.'+qp') .AND.(TEXTC(I1).NE.'-qp') .AND.
C    *        (TEXTC(I1).NE.'+qs') .AND.(TEXTC(I1).NE.'-qs') .AND.
C    *        (TEXTC(I1).NE.'+a11').AND.(TEXTC(I1).NE.'-a11').AND.
C    *        (TEXTC(I1).NE.'+a12').AND.(TEXTC(I1).NE.'-a12').AND.
C    *        (TEXTC(I1).NE.'+a22').AND.(TEXTC(I1).NE.'-a22').AND.
C    *        (TEXTC(I1).NE.'+a13').AND.(TEXTC(I1).NE.'-a13').AND.
C    *        (TEXTC(I1).NE.'+a23').AND.(TEXTC(I1).NE.'-a23').AND.
C    *        (TEXTC(I1).NE.'+a33').AND.(TEXTC(I1).NE.'-a33').AND.
C    *        (TEXTC(I1).NE.'+a14').AND.(TEXTC(I1).NE.'-a14').AND.
C    *        (TEXTC(I1).NE.'+a24').AND.(TEXTC(I1).NE.'-a24').AND.
C    *        (TEXTC(I1).NE.'+a34').AND.(TEXTC(I1).NE.'-a34').AND.
C    *        (TEXTC(I1).NE.'+a44').AND.(TEXTC(I1).NE.'-a44').AND.
C    *        (TEXTC(I1).NE.'+a15').AND.(TEXTC(I1).NE.'-a15').AND.
C    *        (TEXTC(I1).NE.'+a25').AND.(TEXTC(I1).NE.'-a25').AND.
C    *        (TEXTC(I1).NE.'+a35').AND.(TEXTC(I1).NE.'-a35').AND.
C    *        (TEXTC(I1).NE.'+a45').AND.(TEXTC(I1).NE.'-a45').AND.
C    *        (TEXTC(I1).NE.'+a55').AND.(TEXTC(I1).NE.'-a55').AND.
C    *        (TEXTC(I1).NE.'+a16').AND.(TEXTC(I1).NE.'-a16').AND.
C    *        (TEXTC(I1).NE.'+a26').AND.(TEXTC(I1).NE.'-a26').AND.
C    *        (TEXTC(I1).NE.'+a36').AND.(TEXTC(I1).NE.'-a36').AND.
C    *        (TEXTC(I1).NE.'+a46').AND.(TEXTC(I1).NE.'-a46').AND.
C    *        (TEXTC(I1).NE.'+a56').AND.(TEXTC(I1).NE.'-a56').AND.
C    *        (TEXTC(I1).NE.'+a66').AND.(TEXTC(I1).NE.'-a66').AND.
C    *        (TEXTC(I1).NE.'+q11').AND.(TEXTC(I1).NE.'-q11').AND.
C    *        (TEXTC(I1).NE.'+q12').AND.(TEXTC(I1).NE.'-q12').AND.
C    *        (TEXTC(I1).NE.'+q22').AND.(TEXTC(I1).NE.'-q22').AND.
C    *        (TEXTC(I1).NE.'+q13').AND.(TEXTC(I1).NE.'-q13').AND.
C    *        (TEXTC(I1).NE.'+q23').AND.(TEXTC(I1).NE.'-q23').AND.
C    *        (TEXTC(I1).NE.'+q33').AND.(TEXTC(I1).NE.'-q33').AND.
C    *        (TEXTC(I1).NE.'+q14').AND.(TEXTC(I1).NE.'-q14').AND.
C    *        (TEXTC(I1).NE.'+q24').AND.(TEXTC(I1).NE.'-q24').AND.
C    *        (TEXTC(I1).NE.'+q34').AND.(TEXTC(I1).NE.'-q34').AND.
C    *        (TEXTC(I1).NE.'+q44').AND.(TEXTC(I1).NE.'-q44').AND.
C    *        (TEXTC(I1).NE.'+q15').AND.(TEXTC(I1).NE.'-q15').AND.
C    *        (TEXTC(I1).NE.'+q25').AND.(TEXTC(I1).NE.'-q25').AND.
C    *        (TEXTC(I1).NE.'+q35').AND.(TEXTC(I1).NE.'-q35').AND.
C    *        (TEXTC(I1).NE.'+q45').AND.(TEXTC(I1).NE.'-q45').AND.
C    *        (TEXTC(I1).NE.'+q55').AND.(TEXTC(I1).NE.'-q55').AND.
C    *        (TEXTC(I1).NE.'+q16').AND.(TEXTC(I1).NE.'-q16').AND.
C    *        (TEXTC(I1).NE.'+q26').AND.(TEXTC(I1).NE.'-q26').AND.
C    *        (TEXTC(I1).NE.'+q36').AND.(TEXTC(I1).NE.'-q36').AND.
C    *        (TEXTC(I1).NE.'+q46').AND.(TEXTC(I1).NE.'-q46').AND.
C    *        (TEXTC(I1).NE.'+q56').AND.(TEXTC(I1).NE.'-q56').AND.
C    *        (TEXTC(I1).NE.'+q66').AND.(TEXTC(I1).NE.'-q66').AND.
C    *        (TEXTC(I1).NE.'srf' )
     *        ) THEN
C           WFSRF-07
            CALL ERROR('WFSRF-07: Wrong value of COLUMNii.')
C           See allowed values of COLUMNii in the
C           description of file SEP.
          ENDIF
          I1=I1+1
          IF (I1.GT.69) THEN
            CALL RSEP3T('COLUMN70',TEXT,' ')
            IF (TEXT.NE.' ') THEN
C             WFSRF-08
              CALL ERROR('WFSRF-08: More than 69 COLUMNs.')
C             Currently up to 69 values of COLUMNii may be specified.
C             See allowed values of COLUMNii in the
C             description of file SEP.
            ENDIF
          ELSE
            GOTO 5
          ENDIF
        ENDIF
c       NQPOI=I1-1+1
        NQPOI=7
C     End of the loop over future columns of the output file.
C
C
C     Reading filenames of the output files, recording
C     which quantities are to be written into the files:
      CHOUT(1)='MTT'
      CHOUT(2)='MP1'
      CHOUT(3)='MP2'
      CHOUT(4)='MP3'
      NOUT=4
C     CALL RSEP3T('NUM',FOUT(0),'mtt-num.out')
C     CALL RSEP3T('MTT',FOUT(1),'mtt-tt.out')
C     IF (FOUT(1).NE.' ') THEN
C       CHOUT(1)='MTT'
C       CHOUT(2)=' '
C       CHOUT(3)=' '
C       CHOUT(4)=' '
C       NOUT=4
C       CALL RSEP3T('MP1',FOUT(2),' ')
C       IF (FOUT(2).NE.' ') CHOUT(2)='MP1'
C       CALL RSEP3T('MP2',FOUT(3),' ')
C       IF (FOUT(3).NE.' ') CHOUT(3)='MP2'
C       CALL RSEP3T('MP3',FOUT(4),' ')
C       IF (FOUT(4).NE.' ') CHOUT(4)='MP3'
C     ELSE
C       NOUT=0
C       CALL RSEP3T('MP1',FOUT(NOUT+1),' ')
C       IF (FOUT(NOUT+1).NE.' ') THEN
C         CHOUT(NOUT+1)='MP1'
C         NOUT=NOUT+1
C       ENDIF
C       CALL RSEP3T('MP2',FOUT(NOUT+1),' ')
C       IF (FOUT(NOUT+1).NE.' ') THEN
C         CHOUT(NOUT+1)='MP2'
C         NOUT=NOUT+1
C       ENDIF
C       CALL RSEP3T('MP3',FOUT(NOUT+1),' ')
C       IF (FOUT(NOUT+1).NE.' ') THEN
C         CHOUT(NOUT+1)='MP3'
C         NOUT=NOUT+1
C       ENDIF
C     ENDIF
C     CALL RSEP3T('MHIST',FOUT(NOUT+1),' ')
C     IF (FOUT(NOUT+1).NE.' ') THEN
C       CHOUT(NOUT+1)=  'MHIST'
C       NOUT=NOUT+1
C     ENDIF
C
C     CALL MTTQ
C
      NQ=5+NOUT
      RETURN
C
C=======================================================================
C
      ENTRY WFQRI
C
C.......................................................................
C
C Entry designed to read the output of the CRT at an initial point of
C a ray.
C
C     Travel time:
      RAM(NRAMP+(5+1))=YI(1)
C     Three components of the slowness vector:
      RAM(NRAMP+(5+1+1))=YI(6)
      RAM(NRAMP+(5+1+2))=YI(7)
      RAM(NRAMP+(5+1+3))=YI(8)
C     DO 10, I1=1,NOUT
C       IF (CHOUT(I1).EQ.'MTT') THEN
C         Travel time:
C         RAM(NRAMP+(5+I1))=YI(1)
C         Three components of the slowness vector:
C         RAM(NRAMP+(5+I1+1))=YI(6)
C         RAM(NRAMP+(5+I1+2))=YI(7)
C         RAM(NRAMP+(5+I1+3))=YI(8)
C       ELSEIF (CHOUT(I1).EQ.'MP1') THEN
C         First component of the slowness vector:
C         RAM(NRAMP+(5+I1))=YI(6)
C       ELSEIF (CHOUT(I1).EQ.'MP2') THEN
C         Second component of the slowness vector:
C         RAM(NRAMP+(5+I1))=YI(7)
C       ELSEIF (CHOUT(I1).EQ.'MP3') THEN
C         Third component of the slowness vector:
C         RAM(NRAMP+(5+I1))=YI(8)
C       ELSEIF (CHOUT(I1).EQ.'MHIST') THEN
C         Ray history:
C         IRAM(NRAMP+(5+I1))=ISHEET
C       ELSEIF (CHOUT(I1).EQ.' ') THEN
C         Nothing to do:
C         CONTINUE
C       ELSE
C         User-defined quantity:
C         RAM(NRAMP+(5+I1))=0.
C         CALL MTTQRI(CHOUT(I1),RAM(NRAMP+(5+I1)))
C       ENDIF
C 10  CONTINUE
      NRAMP=NRAMP+NQ
      RETURN
C=======================================================================
C
      ENTRY WFQRA
C
C.......................................................................
C
C Entry designed to read the output of the CRT at other than initial
C point of a ray.
C
C     Travel time:
      RAM(NRAMP+(5+1))=Y(1)
C     Three components of the slowness vector:
      RAM(NRAMP+(5+1+1))=Y(6)
      RAM(NRAMP+(5+1+2))=Y(7)
      RAM(NRAMP+(5+1+3))=Y(8)
C     DO 20, I1=1,NOUT
C       IF (CHOUT(I1).EQ.'MTT') THEN
C         Travel time:
C         RAM(NRAMP+(5+I1))=Y(1)
C         Three components of the slowness vector:
C         RAM(NRAMP+(5+I1+1))=Y(6)
C         RAM(NRAMP+(5+I1+2))=Y(7)
C         RAM(NRAMP+(5+I1+3))=Y(8)
C       ELSEIF (CHOUT(I1).EQ.'MHIST') THEN
C         Ray history:
C         IRAM(NRAMP+(5+I1))=ISHEET
C       ELSEIF (CHOUT(I1).EQ.'MP1') THEN
C         First component of the slowness vector:
C         RAM(NRAMP+(5+I1))=Y(6)
C       ELSEIF (CHOUT(I1).EQ.'MP2') THEN
C         Second component of the slowness vector:
C         RAM(NRAMP+(5+I1))=Y(7)
C       ELSEIF (CHOUT(I1).EQ.'MP3') THEN
C         Third component of the slowness vector:
C         RAM(NRAMP+(5+I1))=Y(8)
C       ELSEIF (CHOUT(I1).EQ.' ') THEN
C         Nothing to do:
C         CONTINUE
C       ELSE
C         User-defined quantity:
C         RAM(NRAMP+(5+I1))=0.
C         CALL MTTQRA(CHOUT(I1),RAM(NRAMP+(5+I1)))
C       ENDIF
C 20  CONTINUE
      NRAMP=NRAMP+NQ
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WFINTP(IPA,IPB,TI,X1I,X2I,X3I,P1I,P2I,P3I)
C
C----------------------------------------------------------------------
      INTEGER IPA,IPB
      REAL TI,X1I,X2I,X3I,P1I,P2I,P3I
C
C Subroutine for linear interpolation between two points A and B.
C
C----------------------------------------------------------------------
C
C Input:
C     IPA,IPB ... Indices of positions of points A and B in array RAM.
C             The points must be in the same complex block.
C     TI ...  Time for interpolation. Time TI is supposed to be between
C             the travel times TA and TB in points A and B.
C
C Output:
C     X1I,X2I,X3I ... Three coordinates of interpolated point.
C     P1I,P2I,P3I ... Three components of slowness vector
C             in the interpolated point.
C
C.......................................................................
C
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C For the description of the parameters stored in the common block
C WFSRFC refer to the file wfsrf.inc.
C ...........................
C Auxiliary storage locations:
      REAL WA,WB,AUX
      REAL TA,TB,X1A,X2A,X3A,P1A,P2A,P3A,X1B,X2B,X3B,P1B,P2B,P3B
C-----------------------------------------------------------------------
C
      IF (.NOT.(IRAM(IPA+5).EQ.IRAM(IPB+5))) THEN
C       WFSRF-09
        CALL ERROR('WFSRF-09: Different indices of complex blocks.')
C       This error should not appear.
C       Please contact the author.
      ENDIF
C
      TA =RAM(IPA+6)
      X1A=RAM(IPA+1)
      X2A=RAM(IPA+2)
      X3A=RAM(IPA+3)
      P1A=RAM(IPA+7)
      P2A=RAM(IPA+8)
      P3A=RAM(IPA+9)
      TB =RAM(IPB+6)
      X1B=RAM(IPB+1)
      X2B=RAM(IPB+2)
      X3B=RAM(IPB+3)
      P1B=RAM(IPB+7)
      P2B=RAM(IPB+8)
      P3B=RAM(IPB+9)
C
      AUX=TB-TA
      IF ((TB.EQ.TI).OR.(AUX.EQ.0.)) THEN
        X1I=X1B
        X2I=X2B
        X3I=X3B
        P1I=P1B
        P2I=P2B
        P3I=P3B
      ELSEIF (TA.EQ.TI) THEN
        X1I=X1A
        X2I=X2A
        X3I=X3A
        P1I=P1A
        P2I=P2A
        P3I=P3A
      ELSE
        WB=(TI-TA)/AUX
        WA=1.-WB
        X1I=WA*X1A+WB*X1B
        X2I=WA*X2A+WB*X2B
        X3I=WA*X3A+WB*X3B
        P1I=WA*P1A+WB*P1B
        P2I=WA*P2A+WB*P2B
        P3I=WA*P3A+WB*P3B
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WFSIDE(IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2,
     *                  I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT)
C
C.......................................................................
C
C Subroutine for search for points between rays at a side of a ray tube.
C
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C For the description of the parameters stored in the common block
C WFSRFC refer to the file wfsrf.inc.
C ...........................
      INTEGER IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2,
     *                  I1MI,I1MA,I2MI,I2MA,NEXT
      REAL WFTIM
      REAL TMIN,TMAX
      INTEGER IPA,IPB
      REAL X1I,X2I,X3I,P1I,P2I,P3I
      SAVE    IPA,IPB
C.......................................................................
C
      IF (NEXT.EQ.0) THEN
C       Search for first point:
        IF (KRAY1.EQ.0) THEN
C         The point must be at the top of the tube:
          IPA=IRAM(IRAY1-ORAYA)-NQ
          IPB=IRAM(IRAY2-ORAYA)-NQ
          TMIN=AMIN1(RAM(IPA+6),RAM(IPB+6))
          TMAX=AMAX1(RAM(IPA+6),RAM(IPB+6))
          IF (.NOT.((TMIN.LE.WFTIM).AND.(WFTIM.LE.TMAX))) THEN
C           WFSRF-10
            CALL ERROR('WFSRF-10: Wrong invocation of WFSIDE.')
C           This error should not appear.
C           Please contact the author.
          ENDIF
          CALL WFINTP(IPA,IPB,WFTIM,X1I,X2I,X3I,P1I,P2I,P3I)
          CALL WFWRPT(IRAM(IPA+5),X1I,X2I,X3I,P1I,P2I,P3I)
          IF (IRAM(IPA+5).NE.ICBR2) THEN
C           The new point is in another block than the point on ray 2,
C           there must be some other points:
            NEXT=-1
          ENDIF
        ELSE
C         The point must be within (or at the top of) the tube.
C         Starting from the source:
          IPA=I1MI
          IPB=I2MI
          NEXT=1
          CALL WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT)
          IF (IRAM(IPA+5).NE.ICBR1) THEN
C           Starting from the top of the tube:
            IPA=I1MA
            IPB=I2MA
            NEXT=-1
            CALL WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT)
            IF (IRAM(IPA+5).NE.ICBR1) THEN
C             WFSRF-11
              CALL ERROR('WFSRF-11: Point in propper block not found.')
C             The point in the same complex block were found neither
C             in the upwards search, nor in the downwards search.
C             This error should not appear.
C             Please contact the author.
            ENDIF
          ENDIF
          CALL WFINTP(IPA,IPB,WFTIM,X1I,X2I,X3I,P1I,P2I,P3I)
          CALL WFWRPT(IRAM(IPA+5),X1I,X2I,X3I,P1I,P2I,P3I)
          IF ((IRAM(IPA+5).EQ.ICBR2).OR.
     *        ((ICBR2.EQ.0).AND.(IRAM(IPA+5).EQ.IRAM(I2MA+5)))) THEN
C           The new point is in the same block as the point on ray B:
            NEXT=0
          ENDIF
        ENDIF
      ELSE
C       Search upwards or downwards the tube:
        IPA=IPA+NEXT*NQ
        IPB=IPB+NEXT*NQ
        CALL WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT)
        CALL WFINTP(IPA,IPB,WFTIM,X1I,X2I,X3I,P1I,P2I,P3I)
        CALL WFWRPT(IRAM(IPA+5),X1I,X2I,X3I,P1I,P2I,P3I)
        IF ((IRAM(IPA+5).EQ.ICBR2).OR.
     *      ((ICBR2.EQ.0).AND.(IRAM(IPA+5).EQ.IRAM(I2MA+5)))) THEN
C         The new point is in the same block as the point on ray B:
          NEXT=0
        ENDIF
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WFWRPT(ICB,X1,X2,X3,P1,P2,P3)
C
C.......................................................................
C
C Subroutine to write a point on a wavefront to array (I)RAM.
C
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C For the description of the parameters stored in the common block
C WFSRFC refer to the file wfsrf.inc.
C............................
      INTEGER ICB
      REAL X1,X2,X3,P1,P2,P3
C
      IF (NPOI+NQPOI.GT.NPLGN) CALL WFEROR(2)
      IRAM(NPOI+1)=ICB
      RAM(NPOI+2)=X1
      RAM(NPOI+3)=X2
      RAM(NPOI+4)=X3
      RAM(NPOI+5)=P1
      RAM(NPOI+6)=P2
      RAM(NPOI+7)=P3
      NPOI=NPOI+NQPOI
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WFLINE(IPA,IPB,I1MI,I1MA,I2MI,I2MA,WFTIM,NEXT)
C
C.......................................................................
C
C Subroutine for a search for a line at a side of a ray tube.
C
C Common block /WFSRFC/ to store triangles and rays:
      INCLUDE 'wfsrf.inc'
C For the description of the parameters stored in the common block
C WFSRFC refer to the file wfsrf.inc.
C............................
      INTEGER IPA,IPB,I1MI,I1MA,I2MI,I2MA,NEXT
      REAL WFTIM
C
      INTEGER NQL
      REAL TMIN,TMAX
C
      IF ((NEXT.NE.1).AND.(NEXT.NE.-1)) THEN
C       WFSRF-12
        CALL ERROR('WFSRF-12: Wrong invocation of WFLINE.')
C       This error should not appear.
C       Please contact the author.
      ENDIF
      NQL=NEXT*NQ
C
  10  CONTINUE
        TMIN=AMIN1(RAM(IPA+6),RAM(IPB+6))
        TMAX=AMAX1(RAM(IPA+6),RAM(IPB+6))
        IF ((TMIN.LE.WFTIM).AND.(WFTIM.LE.TMAX)) THEN
C         The point is between IPA and IPB:
          RETURN
        ENDIF
  20    CONTINUE
          IPA=IPA+NQL
          IF ((IPA.LT.I1MI).OR.(IPA.GT.I1MA)) THEN
C           WFSRF-13
            CALL ERROR('WFSRF-13: Point A not reached.')
C           The point on first ray not reached.
C           This error should not appear.
C           Please contact the author.
          ENDIF
        IF (.NOT.((IRAM(IPA+4).NE.0).OR.
     *            (IPA.EQ.I1MI).OR.(IPA.EQ.I1MA))) GOTO 20
  30    CONTINUE
          IPB=IPB+NQL
          IF ((IPB.LT.I2MI).OR.(IPB.GT.I2MA)) THEN
C           WFSRF-14
            CALL ERROR('WFSRF-14: Point B not reached.')
C           The point on second ray not reached.
C           This error should not appear.
C           Please contact the author.
          ENDIF
        IF (.NOT.((IRAM(IPB+4).NE.0).OR.
     *            (IPB.EQ.I2MI).OR.(IPB.EQ.I2MA))) GOTO 30
      GOTO 10
      END
C
C=======================================================================
C
      SUBROUTINE WFEROR(IERR)
C
C.......................................................................
C
C Subroutine to write error messages.
C
      INTEGER IERR
      IF     (IERR.EQ.1) THEN
C       WFSRF-15
        CALL ERROR('WFSRF-15: Small array RAM.')
C       Dimensions of array RAM exceeded when writing points
C       on rays.
C       Try to enlarge the dimension MRAM in the file
C       ram.inc.
      ELSEIF (IERR.EQ.2) THEN
C       WFSRF-16
        CALL ERROR('WFSRF-16: Small array RAM.')
C       Dimensions of array RAM exceeded when writing points
C       at wavefronts.
C       Try to enlarge the dimension MRAM in the file
C       ram.inc.
      ELSEIF (IERR.EQ.3) THEN
C       WFSRF-17
        CALL ERROR('WFSRF-17: Small array RAM.')
C       Dimensions of array RAM exceeded when writing polygons
C       at wavefronts.
C       Try to enlarge the dimension MRAM in the file
C       ram.inc.
      ELSEIF (IERR.EQ.4) THEN
C       WFSRF-18
        CALL ERROR('WFSRF-18: Small array KPOL.')
C       This error may be caused by small dimension MPOL of array
C       KPOL defined in subroutine WFTUBE. Try to enlarge MPOL.
      ELSE
C       WFSRF-19
        CALL ERROR('WFSRF-19: Wrong IERR.')
C       This subroutine was invocated with wrong value of IERR.
C       This error should not appear.
C       Please contact the author.
      ENDIF
      END
C
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'indexx.for'
C     indexx.for
C
C=======================================================================
C