C
C Program WFSRF to generate wavefronts composed of polygons for display
C purposes
C
C Version: 5.60
C Preliminary version tested only within history file len-wf.h.
C                    len-wf.h
C Date: 2001, October 24
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.  Only the wavefront
C corresponding to the first elementary wave stored in the CRT output
C file CRT-T is generated.
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.  For general description of file CRTOUT refer
C             to file writ.for.
C             A single set of CRT output files is read from file CRTOUT.
C             Only files CRT-R,
C             CRT-I and
C             CRT-T are read by WFSRF.
C             File CRT-R must contain all rays traced by CRT.
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               'MTT' .. Travel time of the wavefront.
C               'MHIST'. Ray history in the vertex.
C               'MPQ11'. Component 1,1 of the 4x4 ray propagator matrix.
C               'MPQ21'. Component 2,1 of the 4x4 ray propagator matrix.
C               'MPQ31'. Component 3,1 of the 4x4 ray propagator matrix.
C               'MPQ41'. Component 4,1 of the 4x4 ray propagator matrix.
C               'MPQ12'. Component 1,2 of the 4x4 ray propagator matrix.
C               'MPQ22'. Component 2,2 of the 4x4 ray propagator matrix.
C               'MPQ32'. Component 3,2 of the 4x4 ray propagator matrix.
C               'MPQ42'. Component 4,2 of the 4x4 ray propagator matrix.
C               'MPQ13'. Component 1,3 of the 4x4 ray propagator matrix.
C               'MPQ23'. Component 2,3 of the 4x4 ray propagator matrix.
C               'MPQ33'. Component 3,3 of the 4x4 ray propagator matrix.
C               'MPQ43'. Component 4,3 of the 4x4 ray propagator matrix.
C               'MPQ14'. Component 1,4 of the 4x4 ray propagator matrix.
C               'MPQ24'. Component 2,4 of the 4x4 ray propagator matrix.
C               'MPQ34'. Component 3,4 of the 4x4 ray propagator matrix.
C               'MPQ44'. Component 4,4 of the 4x4 ray propagator matrix.
C               'GBW11'. Sum of squares of Gaussian beam widths
C                 corresponding to GBY11. The sum of squares of Gaussian
C                 beam widths is defined here as the trace of the
C                 inverse imaginary part of the matrix of second
C                 derivatives of the travel time of the Gaussian beam.
C                 The sum of squares of Gaussian beam widths depends on
C                 parameters GBEij, GBR11, GBR12, GBR22, GBY11 and GBY22
C                 described below.  It may be expressed as sum W11+W22,
C                 where W11 is dependent on GBY11 but independent
C                 on GBY22, and W22 is dependent on GBY22 but
C                 independent on GBY11. Values of W11 are written to
C                 the column given by COLUMNii='GBW11'.
C                 Note that positive GBY11 must be specified for
C                 this option.
C               'GBW1'.. Values of SQRT(W11). See the description of
C                 COLUMNii='GBW11'. Positive GBY11 must be specified.
C               'GBW22'.. Values of W22. See the description of
C                 COLUMNii='GBW11'. Positive GBY22 must be specified.
C               'GBW2'.. Values of SQRT(W22). See the description of
C                 COLUMNii='GBW11'. Positive GBY22 must be specified.
C               'AMP' .. Amplitude of the Green function in the vertex.
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 Numerical parameters more precisely specifying the above quantities:
C     Refer to the list in program mtt.for.
C Numerical parameters specifying optional ray-parameter grid:
C     Refer to the list in program mtt.for.
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL WFEROR,WFREAD,WFERAS,WFTUBE,WFINTP,
     *         WFQD,WFQRI,WFQRA,WFSIDE,WFLINE,
     *MTTQ,
     *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 ... This file.
C     MTTQ ... File mttq.for.
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,J3,NPOINT,LEN1,LEN2,LENG
      INTEGER MLJPOL
      PARAMETER (MLJPOL=10)
      REAL AUX,Z(MOUT),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
      READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
  10  CONTINUE
        READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
        IF (IRAY1.EQ.0) GOTO 20
        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.MAXRAM) 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.MAXRAM) 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=MAXRAM-KRAMP*(NWFT*NT*5 + NWFT*NR*NQ)
      IF (MRAMP.LE.NRAMP) CALL WFEROR(1)
      NPOI=MRAMP
      NPLGN=MAXRAM+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(MAXRAM/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     Unifying character strings in TEXTC and CHOUT:
      DO 110, I1=1,NQ
        CALL LOWER(CHOUT(I1))
        IF (CHOUT(I1).EQ.'mp1') CHOUT(I1)='norm1'
        IF (CHOUT(I1).EQ.'mp2') CHOUT(I1)='norm2'
        IF (CHOUT(I1).EQ.'mp3') CHOUT(I1)='norm3'
  110 CONTINUE
C
C     Storing points:
      OPEN(LU1,FILE=VRTX,FORM='FORMATTED')
      WRITE(LU1,'(A)') '/'
      NPOINT=(NPOI-MRAMP)/NQ
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)*NQ
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.MOUT)) THEN
            J2=J2+1
            DO 210, I2=1,NQ
              IF (TEXTC(J2).EQ.CHOUT(I2)) THEN
                J3=I2
                GOTO 211
              ENDIF
  210       CONTINUE
C           WFSRF-23
            CALL ERROR ('WFSRF-23: No CHOUT found for TEXTC.')
C           This error should not appear.
C           Please contact the author.
  211       CONTINUE
            IF ((TEXTC(J2).EQ.'icb').OR.(TEXTC(J2).EQ.'mhist')) THEN
              AUX=FLOAT(IRAM(J1+J3))
            ELSE
              AUX=RAM(J1+J3)
            ENDIF
            IF (POWER(J2).EQ.1.) THEN
              Z(J2)=AUX
            ELSE
              Z(J2)=AUX**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=MAXRAM-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)/NQ,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)/NQ
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
        IWAVE0=IWAVE
        IF ((IRAY0.EQ.0).AND.(IWAVE0.EQ.0)) THEN
C         Reading the first point on a ray of the first wave:
          CALL AP00(0,LU1,LU2)
          IF (IWAVE.LT.1) GOTO 20
        ELSEIF (IWAVE.EQ.IWAVE0) THEN
C         Reading the next point on a ray of the actual wave:
          CALL AP00(0,LU1,LU2)
          IF (IWAVE.NE.IWAVE0) GOTO 20
        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:
            CALL WFQRI
          ENDIF
C         Recording the new point:
          CALL WFQRA
        ENDIF
        IRAM(IRAY-ORAYA)=NRAMP
        IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA))
     *     RETURN
      GOTO 10
C
  20  CONTINUE
C     End of rays:
      IF (IRAY0.LT.IRAYMA) THEN
C       Some rays at the end of the CRT output file missing:
        DO 22, I1=IRAY0+1,IRAYMA
          IF (I1.GE.IRAYMI) THEN
            IRAM(I1-ORAYA)=NRAMP
          ENDIF
  22    CONTINUE
      ENDIF
      RETURN
      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 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)
      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(I1))
                KRAY1=NPOI
                ICBR1=IRAM(KRAY1-NQ+5)
                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-NQ+5)
        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(I1))
                KRAY2=NPOI
                ICBR2=IRAM(KRAY2-NQ+5)
                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-NQ+5)
        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(I1))
                  KRAY3=NPOI
                  ICBR3=IRAM(KRAY3-NQ+5)
                  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-NQ+5)
          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)-NQ+5).NE.IRAM(KPOL(I2)-NQ+5))
     *        LSHIFT=.TRUE.
  70      CONTINUE
        ENDIF
        IF (LSHIFT) THEN
  72      CONTINUE
            IF (IRAM(KPOL(1)-NQ+5).EQ.IRAM(KPOL(NPOL)-NQ+5)) 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)-NQ+5).NE.IRAM(KPOL(J2-1)-NQ+5))
     *           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,I2
      CHARACTER*20 FORMA1,FORMA2,FORMA3
      CHARACTER*5 QNAMEL(MOUT),TEXT,TEXTL
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     Reading the list of the quantities which may be interpolated:
      CALL MTTQ
      DO 3, I1=1,MOUT
        QNAMEL(I1)=QNAMES(I1)
        IF (QNAMEL(I1).NE.' ') THEN
          CALL LOWER(QNAMEL(I1))
        ENDIF
  3   CONTINUE
C
C     Recalling the data specifying the quantities to be written into
C     the individual columns of the output file VRTX with points
C     wavefronts:
      FORMA1='COLUMN01'
      FORMA2='POWER01'
      FORMA3='IVALUE01'
      I1=1
   5  CONTINUE
        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').AND.
     *        (TEXTC(I1).NE.'mtt').AND.
     *        (TEXTC(I1).NE.'mhist')) THEN
            DO 7, I2=1,MOUT
              IF (TEXTC(I1).EQ.QNAMEL(I2)) GOTO 8
  7         CONTINUE
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
  8       CONTINUE
          I1=I1+1
          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.GT.MOUT) THEN
            CALL RSEP3T(FORMA1,TEXT,' ')
            IF (TEXT.NE.' ') THEN
C             WFSRF-08
              CALL ERROR('WFSRF-08: More than MOUT COLUMNs.')
C             Currently up to MOUT 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     End of the loop over future columns of the output file.
C
C
C     Recording which quantities are to be computed and stored
C     in the points on the rays and in the points on the wavefronts:
      CHOUT(1)='X1'
      CHOUT(2)='X2'
      CHOUT(3)='X3'
      CHOUT(4)='ISRF'
      CHOUT(5)='ICB'
      CHOUT(6)='MTT'
      CHOUT(7)='MP1'
      CHOUT(8)='MP2'
      CHOUT(9)='MP3'
      NQ=9
      DO 10, I1=1,MOUT
        IF (TEXTC(I1).EQ.'mhist') THEN
            NQ=NQ+1
            IF (NQ.GT.MOUT) CALL WFEROR(5)
            CHOUT(NQ)='MHIST'
        ENDIF
  10  CONTINUE
      DO 20, I2=1,MOUT
        IF (QNAMES(I2).EQ.' ') GOTO 21
        TEXT=QNAMES(I2)
        TEXTL=QNAMEL(I2)
        DO 15, I1=1,MOUT
          IF (TEXTC(I1).EQ.TEXTL) THEN
            NQ=NQ+1
            IF (NQ.GT.MOUT) CALL WFEROR(5)
            CHOUT(NQ)=TEXT
            CALL MTTQD(TEXT)
          ENDIF
  15    CONTINUE
  20  CONTINUE
  21  CONTINUE
C
      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     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:
      RAM(NRAMP+6)=YI(1)
C     Three components of the slowness vector:
      RAM(NRAMP+7)=YI(6)
      RAM(NRAMP+8)=YI(7)
      RAM(NRAMP+9)=YI(8)
C     Other quantities to be interpolated:
      DO 30, I1=10,NQ
        IF (CHOUT(I1).EQ.'MHIST') THEN
C         Ray history:
          IRAM(NRAMP+I1)=ISHEET
        ELSE
C         User-defined quantity:
          RAM(NRAMP+I1)=0.
          CALL MTTQRI(CHOUT(I1),RAM(NRAMP+I1))
        ENDIF
  30  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     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:
      RAM(NRAMP+6)=Y(1)
C     Three components of the slowness vector:
      RAM(NRAMP+7)=Y(6)
      RAM(NRAMP+8)=Y(7)
      RAM(NRAMP+9)=Y(8)
C     Other quantities to be interpolated:
      DO 40, I1=10,NQ
        IF (CHOUT(I1).EQ.'MHIST') THEN
C         Ray history:
          IRAM(NRAMP+I1)=ISHEET
        ELSE
C         User-defined quantity:
          RAM(NRAMP+I1)=0.
          CALL MTTQRA(CHOUT(I1),RAM(NRAMP+I1))
        ENDIF
  40  CONTINUE
      NRAMP=NRAMP+NQ
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WFINTP(IPA,IPB,TI)
C
C----------------------------------------------------------------------
      INTEGER IPA,IPB
      REAL TI
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 No output.
C     The quantities in the interpolated point are written into the
C     array (I)RAM.
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:
      INTEGER I1
      REAL WA,WB,AUX
      REAL TA,TB
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
      IF (NPOI+NQ.GT.NPLGN) CALL WFEROR(2)
C
      TA =RAM(IPA+6)
      TB =RAM(IPB+6)
C
      AUX=TB-TA
      IF ((TB.EQ.TI).OR.(AUX.EQ.0.)) THEN
        WA=0.
        WB=1.
      ELSEIF (TA.EQ.TI) THEN
        WA=1.
        WB=0.
      ELSE
        WB=(TI-TA)/AUX
        WA=1.-WB
      ENDIF
C
      DO 10, I1=1,NQ
        IF     (CHOUT(I1).EQ.'ISRF') THEN
          IF (IRAM(IPA+I1).EQ.IRAM(IPB+I1)) THEN
            IRAM(NPOI+I1)=IRAM(IPA+I1)
          ELSE
            IRAM(NPOI+I1)=0
          ENDIF
        ELSEIF (CHOUT(I1).EQ.'ICB') THEN
          IRAM(NPOI+I1)=IRAM(IPA+I1)
        ELSEIF (CHOUT(I1).EQ.'MHIST') THEN
          IRAM(NPOI+I1)=IRAM(IPA+I1)
        ELSE
          RAM(NPOI+I1)=WA*RAM(IPA+I1)+WB*RAM(IPB+I1)
        ENDIF
  10  CONTINUE
C
      NPOI=NPOI+NQ
      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----------------------------------------------------------------------
C
C Input:
C     IRAY1,IRAY2 ... Indices of the two rays.
C     KRAY1,KRAY2 ... Addresses of the interpolated points on the two
C             rays, or zeros if the points do not exist.
C     ICBR1,ICBR2 ... Indices of complex blocks of the above two points.
C     I1MI,I1MA,I2MI,I2MA ... Minimum and maximum addresses of the points
C             of the two rays IRAY1,IRAY2.
C     WFTIM ... Wawefront travel time.
C     NEXT ... Indicator of the direction of the search for the
C             interpolated points at the side of the tube.
C             NEXT=0  ... The search begins.
C             NEXT=1  ... Upwards search.
C             NEXT=-1 ... Downwards search.
C Output:
C     NEXT ... Indicator of the direction of the search for the
C             interpolated points at the side of the tube.
C             NEXT=0  ... The search is terminated.
C             NEXT=1  ... Upwards search.
C             NEXT=-1 ... Downwards search.
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 ...........................
      INTEGER IRAY1,KRAY1,ICBR1,IRAY2,KRAY2,ICBR2,
     *                  I1MI,I1MA,I2MI,I2MA,NEXT
      REAL WFTIM
      REAL TMIN,TMAX
      INTEGER IPA,IPB
      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)
          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)
          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)
        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 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.
      ELSEIF (IERR.EQ.5) THEN
C       WFSRF-20
        CALL ERROR('WFSRF-20: Small array CHOUT.')
C       This error may be caused by small dimension MOUT of array CHOUT.
C       Try to enlarge MOUT in the include file
C       wfsrf.inc.
      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 'mttq.for'
C     mttq.for
      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