C 
C Subroutine file 'writ.for' to create the output files of complete ray
C tracing
C
C Date: 1998, November 9
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of:
C     WRITB...Block data subroutine defining common blocks /WRITT/,
C             /WRITC/ and /WRITN/ to store the input data on the names
C             of output files.
C             WRITB
C     WRIT1...Subroutine called when starting the complete ray tracing
C             program, and when starting the computation of a new
C             elementary wave.
C             WRIT1
C     WRIT2...Subroutine called when starting the complete tracing of a
C             new ray.
C             WRIT2
C     WRIT31..Subroutine storing the computed quantities along a ray
C             (C.R.T.5.5.1).  It is called with constant step STORE (see
C             input data in the file 'ray.for') of the independent
C             variable along the ray, and at the points of intersection
C             with interfaces either before and after the
C             transformation.
C             WRIT31
C     WRIT32..Subroutine storing the computed quantities at specified
C             surfaces (C.R.T.5.5.2).  It is called at the points of
C             intersection of the ray with the specified surfaces.
C             WRIT32
C     WRIT33..Subroutine storing the computed quantities at the
C             endpoints of the individual elementary waves
C             (C.R.T.5.5.3).  It is called at the endpoints of the
C             elements of rays, situated at structural interfaces.
C             WRIT33
C     WRIT4...Subroutine storing the quantities at the initial point of
C             the ray (C.R.T.6.1).  It is called after termination of
C             tracing the ray.
C             WRIT4
C     CHECK...Auxiliary subroutine to WRIT4.
C             CHECK
C     WRIT5...Subroutine called after termination of the computation of
C             an elementary wave, and when terminating the complete ray
C             tracing program.
C             WRIT5
C     FNAME...Subroutine analyzing the specified strings and composing
C             the filename of them.  Auxiliary subroutine to WRIT1.
C             FNAME
C     WRITA...Subroutine writing the indices of rays forming the
C             homogeneous triangles in the ray-parameter domain.
C             WRITA
C
C.......................................................................
C
C Description of data files:
C                                                    
C Input data WRIT to specify the names of the output files with the
C computed quantities:
C     The data are read in by the list directed input (free format).  In
C     the list of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new READ statement).  All
C     input variables are of the type CHARACTER.  The data are read by
C     the subroutine WRIT1.
C     The filenames are generated using several components (strings) in
C     the following way:  Each component (string) is divided into words,
C     each word being followed by just one space.  Here a word is a
C     substring containing no space, preceded and followed by spaces or
C     by the beginning or end of the string.  An empty word lies between
C     two consecutive spaces.  The filename is composed of
C             the 1-st word of the 1-st component,
C             the 1-st word of the 2-nd component,
C             ...
C             the 1-st word of the last component,
C             the 2-nd word of the 1-st component,
C             the 2-nd word of the 2-nd component,
C             ...
C             the last word of the last component.
C     Example:
C             1-st component = 'd/  .out',
C             2-nd component = 'a b',
C             3-rd component = '1 2',
C             resulting filename = 'd/a1b2.out'.
C     Note that only first 16 characters of each filename component are
C     significant.
C (1) TEXTW
C     String describing the data.  Only the first 80 characters of the
C     string are significant.
C (2) K2P,/
C     K2P...  If non-zero, only two-point rays are written to the files
C             with quantities stored along rays.  Storage of quantities
C             at specified surfaces is not affected by this parameter.
C     Default:  K2P=0.
C (3) FN1A,FN1I,FN1T,/
C     FN1A... String containing the first component of the name of the
C             output file with the computed quantities stored along the
C             rays (C.R.T.5.5.1).
C             Description of output file C.R.T.5.5.1
C     FN1I... String containing the first component of the name of the
C             output file with the quantities at the initial points of
C             rays (C.R.T.6.1), related to the output file with the
C             computed quantities stored along the rays (C.R.T.5.5.1).
C             The other components of the filename are common with the
C             related file containing the computed quantities.
C             Description of output file C.R.T.6.1
C     FN1T... String containing the first component of the name of the
C             formatted output file specifying the homogeneous triangles
C             in the ray-parameter domain.  The other components of the
C             filename are common with the related file containing the
C             computed quantities along rays.
C             Description of file with triangles
C     /...    An obligatory slash.
C     Default: all filenames blank.
C (4) FN2A,FN2I,(FN2B(I),I=1,NENDF),/
C     FN2A... String containing the first component of the name of the
C             output file with the computed quantities stored at the
C             specified surfaces (C.R.T.5.5.2).
C             Description of output file C.R.T.5.5.2
C     FN2I... String containing the first component of the name of the
C             output file with the quantities at the initial points of
C             rays (C.R.T.6.1), related to the output file with the
C             computed quantities stored at the specified surfaces
C             (C.R.T.5.5.2).  The other components of the filename are
C             common with the related file containing the computed
C             quantities.
C             Description of output file C.R.T.6.1
C     FN2B(I)... String containing the second component of the name of
C             the output file with the computed quantities stored at the
C             I-th specified surface (C.R.T.5.5.2).  Here NENDF is the
C             number of surfaces specified for storing the computed
C             quantities.
C     /...    An obligatory slash.
C     Default: all filenames blank.
C (5) Either (5.1) or (5.2):
C (5.1) /
C     A slash is fully sufficient if MODCRT.LE.1, see input data
C     'dcrt.dat' described in source code file 'ray.for'.
C     If MODCRT.LE.1, the filenames (5.2) are ignored if specified.
C     If MODCRT.GE.2, the filenames (5.2) are obligatory.
C (5.2) FN3A,FN3I,(FN3B(I),I=1,NELEM),/
C     FN3A... String containing the first component of the name of the
C             output file with the computed quantities stored at the
C             endpoints of the elements of rays (C.R.T.5.5.3).
C             Description of output file C.R.T.5.5.3
C     FN3I... String containing the first component of the name of the
C             output file with the quantities at the initial points of
C             rays (C.R.T.6.1), related to the output file with the
C             computed quantities stored at the endpoints of the
C             elements of rays (C.R.T.5.5.3).  The other components of
C             the filename are common with the related file containing
C             the computed quantities.
C             Description of output file C.R.T.6.1
C     FN3B(I)... String containing the second component of the name of
C             the output file with the computed quantities stored at the
C             endpoints of the I-th element of rays (C.R.T.5.5.3).
C             These strings need not be specified for MODCRT.LE.1, see
C             input data 'dcrt.dat'.  Otherwise, none of FN3B(i)
C             corresponding to the endpoint of a ray element at which
C             quantities are stored (see MODCRT in the input data
C             'dcrt.dat') may equal spaces.
C     /...    An obligatory slash.
C     Default: all filenames blank.
C (6) For each elementary wave IWAVE the following data (5.1):
C (6.1) FN1C,FN2C,FN3C,(FN2D(I),I=1,NENDF),/
C     FN1C... String containing the second component of the name of the
C             output file with the computed quantities stored along the
C             rays (C.R.T.5.5.1).
C             If both the first and second component of the filename are
C             blank, the file is not created.
C     FN2C... String containing the third component of the name of the
C             output file with the computed quantities stored at the
C             specified surfaces (C.R.T.5.5.2).
C             If the second, third and fourth components of the filename
C             are all blank, the file is not created.
C     FN3C... String containing the third component of the name of the
C             output file with the computed quantities stored at the
C             endpoints of the elements of rays (C.R.T.5.5.3).
C             This string may not equal spaces for MODCRT.GE.2, see
C             input data 'dcrt.dat'.
C             If MODCRT.LE.1, see input data 'dcrt.dat', the file is not
C             created.
C     FN2D(I)... String containing the fourth component of the name of
C             the output file with the computed quantities stored at the
C             I-th specified surface (C.R.T.5.5.2).  Here NENDF is the
C             number of surfaces specified for storing the computed
C             quantities.
C     /...    An obligatory slash.
C     Default: all filenames blank.
C Example of data set WRIT two-point rays stored
C Example of data set WRIT all rays stored
C Example of data set WRIT no rays stored
C
C                                                    
C Output files CRT-R with the computed quantities stored along the rays
C (C.R.T.5.5.1):
C     Unformatted output.  The computed quantities are stored with the
C     step STORE (see the input data in the subroutine file 'ray.for')
C     along all rays, and at all points of intersection of rays with
C     structural interfaces either before and after the transformation.
C     For each point - one record containing the following quantities:
C (1) IWAVE,IRAY,NY,ICB1,ISRF,X,(YL(I),I=1,6),(Y(I),I=1,NY)
C     IWAVE...Index of the elementary wave (output of the subroutine
C             CODE1).  Elementary waves are indexed by 1,2,3,... .
C     IRAY... Index of the ray (output of the subroutine RPAR2).  Rays
C             within each elementary wave are indexed by 1,2,3,... .
C             Note that some of the indexed rays need not exist.
C     NY=IY(1)=27+NAMPL... Number of the basic quantities describing the
C             point of the ray, see the file 'ray.for'.
C     ICB1=IY(5)... Index of the complex block in which the stored point
C             of the ray is situated, supplemented by a sign '+' for P
C             wave and sign '-' for S wave.
C     ISRF=IY(6)... Index of the surface at which the endpoint of the
C             computed element of the ray is situated, supplemented by
C             a sign '+' or '-' for the endpoint situated at the
C             positive or negative side of the surface, respectively.
C             Zero inside the complex block.  Note: the sign of this
C             quantity is the exception to the definition in the
C             original paper on C.R.T.
C     X...    Independent variable along a ray (see YY(1) in
C             C.R.T.5.2.2).
C     YL...   Array containing local quantities at the point of the ray
C             (C.R.T.5.5.4).
C             Description of YL
C     Y...    Array containing basic quantities computed along the ray
C             (C.R.T.5.2.1).
C             Description of Y
C     The detailed description of the quantities YL, YY and IY may be
C     found in the file 'ray.for'.
C
C                                                    
C Output files CRT-S with the computed quantities stored at the
C specified surfaces (C.R.T.5.5.2):
C     Unformatted output.  The computed quantities are stored at the
C     points of intersection of rays with the specified surfaces.
C     For each point - one record containing the following quantities:
C (1) IWAVE,IRAY,NY,ICB1,ISRF,X,(YL(I),I=1,6),(Y(I),I=1,NY)
C     The same quantities as in the data set described above, except
C     that the vectorial reduced amplitudes Y(28)-Y(IY(1)) may (but need
C     not) be replaced by the reduced amplitudes YC(1) to YC(NY-27)
C     involving appropriate conversion coefficients, see C.R.T.5.5.4.
C     YC(1) to YC(NY-27) are described in the subroutine CONV (file
C     'trans.for') and in the subroutine WRIT32 (this file).  Note that
C     NY=33 for P wave and NY=39 for S wave at the initial point of the
C     ray.
C
C                                                    
C Output files CRT-E with the computed quantities stored at the
C endpoints of the elements of rays (C.R.T.5.5.3):
C     Unformatted output.  The files, generated only if KSTORE.GE.2 (see
C     input data in the subroutine file 'ray.for'),  are auxiliary disk
C     storage locations to the program CRT and are not intended to be
C     the output of the complete ray tracing used by the user
C     application programs.  The computed quantities are stored at the
C     endpoints of the elements of all rays before the transformation.
C     For each point - one record containing the following quantities:
C (1) IRAY,(IY(I),I=1,12),(YL(I),I=1,6),(Y(I),I=1,IY(1)),(YY(I),I=1,5)
C     IRAY... Index of the ray (output of the subroutine RPAR2).  Rays
C             within each elementary wave are indexed by 1,2,3,... .
C             Note that some of the indexed rays may not exist.
C     YL...   Array containing local quantities at the point of the ray.
C             (C.R.T.5.5.4).
C             Description of YL
C     Y...    Array containing basic quantities computed along the ray.
C             (C.R.T.5.2.1).
C             Description of Y
C     YY...   Array containing real auxiliary quantities computed along
C             the ray (C.R.T.5.2.2).
C             Description of YY
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray (C.R.T.5.2.2).
C             Description of IY
C     The detailed description of the quantities YL, Y, YY and IY may be
C     found in the file 'ray.for'.
C
C                                                    
C Output files CRT-I with the quantities at the initial points of rays
C (C.R.T.6.1):
C     Unformatted output.  The quantities are stored at the initial
C             points of rays.
C     For whole rays (C.R.T.5.5.1), only initial points of two-point
C             rays are stored if only two-point rays are stored.
C     For intersections with surfaces (C.R.T.5.5.2) and endpoints of the
C             elements of rays (C.R.T.5.5.3) initial points of all rays
C             are stored.
C     For initial point - one record containing following 34 quantities:
C (1) (-IWAVE),IRAY,ICB1I,IEND,ISHEET,IREC,(YLI(I),I=1,6),(YI(I),I=1,29)
C     IWAVE...Index of the elementary wave (output of the subroutine
C             CODE1).  Elementary waves are indexed by 1,2,3,... .
C     IRAY... Index of the ray (output of the subroutine RPAR2).  Rays
C             within each elementary wave are indexed by 1,2,3,... .
C             Note that some of the indexed rays may not exist.
C     ICB1I...Index of the complex block in which the initial point of
C             the ray is situated, supplemented by a sign '+' for P wave
C             and sign '-' for S wave (see C.R.T.6.1).
C     IEND... Reason of the termination of the computation of a ray (see
C             C.R.T.5.4).  For a detailed description see subroutines
C             RAY2 (file 'ray.for') and INIT2 (file 'init.for').
C             Description of IEND in RAY2
C             Description of IEND in INIT2
C     ISHEET..Ray-history index.  The different ray histories are
C             consecutively indexed by positive integers 1,2,3,...
C             According to their appearance during ray tracing.
C             The ray histories are indexed independently within each
C             elementary wave.  The indices are the output of subroutine
C             RPAR4 of the file 'rpar.for'.
C             The ray-history indices are complemented with sign:
C             Positive - successful ray (crossing reference surface),
C             negative - unsuccessful ray (terminating before crossing
C             reference surface).
C     IREC... Index of the receiver for a two-point ray, determined in
C             subroutine RPAR4.
C     YLI...  Array containing the values of the quantities YL(1)-YL(6),
C             see C.R.T.5.5.4, describing the local properties of the
C             model at the initial point of the ray.  See the list of
C             YL(1) to YL(6) in the file 'ray.for'.
C             Description of YL
C     YI...   Array containing the quantities describing the properties
C             of the rays and of the travel-time field at the initial
C             point of the ray, see C.R.T.6.1.  These quantities are
C             listed in the file 'init.for'.
C             Description of YI
C
C                                                    
C Output files CRT-T specifying homogeneous triangles in the ray-parameter
C ray-parameter domain:
C     Formatted output.
C     For each homogeneous triangle - one record containing 3 integers:
C (1) IRAY1,IRAY2,IRAY3
C     IRAY1,IRAY2,IRAY3... Indices of rays forming the vertices of the
C             homogeneous triangle.
C             For one-parametric (2-D) systems of rays, triangles are
C             replaced with intervals and IRAY3=0.
C
C.......................................................................
C
C Storage in the memory:
C     The input data WRIT (2) to (6) describing the names of the output
C     files with the computed quantities are stored in the common block
C     /WRITT/.  Other important variables shared by the subroutines of
C     this file are stored in the common blocks /WRITC/ and /WRITN/.
C     The common blocks are defined in the following subroutine:
C     ------------------------------------------------------------------
C     
C
      BLOCK DATA WRITB
      INCLUDE 'writ.inc'
C     writ.inc
      DATA LUW/10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25/
      DATA JENDR/10,21,22,23,24,25,26,30,32,33,40,50,60,61,71,72,73,74/
      END
C     ------------------------------------------------------------------
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT1(LUN,LULOG,IWAVE,IWAVE0,IKODE)
      INTEGER LUN,LULOG,IWAVE,IWAVE0,IKODE
C
C This subroutine is called when starting the complete ray tracing
C program, and when starting the computation of a new elementary wave.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C     LULOG...Logical unit number of the log output device.
C     IWAVE...Zero when starting the complete ray tracing program,
C             otherwise the index of the elementary wave which will be
C             computed (i.e. the output of the subroutine CODE1 from the
C             file 'code.for').
C     IWAVE0..Index of the already computed elementary wave having the
C             most numerous common elements with the current elementary
C             wave.  Need not be defined if IWAVE=0.
C     IKODE...The length of the common part of the codes of the IWAVE-th
C             and IWAVE0-th elementary waves.  Need not be defined if
C             IWAVE=0.
C
C No output.
C
C Common block /DCRT/ (see subroutine file 'ray.for'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common blocks /WRITT/, /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C All the storage locations of the common blocks are defined in this
C subroutine.
C
C Subroutines and external functions required:
      INTEGER LCODE,NELEM
      EXTERNAL WRITB,LCODE,NELEM,SCRO1,FNAME
C     WRITB.. Block data subroutine of this file.
C     LCODE,NELEM... File 'code.for'.
C     SCRO1... File 'scro.for'.
C     FNAME... This file.
C
C Date: 1997, July 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER NKODE,IE1,IE2,I,J
      CHARACTER*16 NAME1(4)
C
      IF(IWAVE.LE.0) THEN
C
C       Reading the name of the input data
        READ(LUN,*) TEXTW
C
        K2P=0
        READ(LUN,*) K2P
C
C       Reading strings composing the names of output files (CRT.5.5.1)
        FN1A=' '
        FN1I=' '
        FN1T=' '
        READ(LUN,*) FN1A,FN1I,FN1T
C       Number of output files to be open
        IF(STORE.EQ.0.) THEN
          NF1=0
        ELSE
          NF1=1
        END IF
C
C       Reading strings composing the names of output files (CRT.5.5.2)
        FN2A=' '
        FN2I=' '
        DO 21 I=1,MF
          FNB(I)=' '
   21   CONTINUE
        READ(LUN,*) FN2A,FN2I,(FNB(I),I=1,MF)
C       Number of output files to be open
        NF2=NSTOR
        IF(2*(NF1+NF2).GT.NUW) THEN
C         553
          CALL ERROR('553 in WRIT1: Few logical units available')
C         There are NUW logical units available, see the common
C         blocks /WRITC/ and /WRITN/ defined in the block data
C         WRITB.  This number should be increased.
        END IF
C
C       Reading strings composing the names of output files (CRT.5.5.3)
        FN3A=' '
        FN3I=' '
        DO 31 I=NF2+1,MF
          FNB(I)=' '
   31   CONTINUE
        READ(LUN,*) FN3A,FN3I,(FNB(I),I=NF2+1,MF)
C       Number of output files to be open
        IF(MODCRT.LE.1) THEN
          MF3=0
        ELSE
          DO 32 I=NF2+1,MF
            IF(FNB(I).EQ.' ') THEN
              MF3=I-1-NF2
              GO TO 33
            END IF
   32     CONTINUE
          MF3=MF-NF2
   33     CONTINUE
        END IF
C
      ELSE
C
C       Numbers of initial ray elements along which nothing is stored
        IF(MODCRT.EQ.0) THEN
          KWRIT1=0
        ELSE
          KWRIT1=NELEM(IKODE)
        END IF
        KWRIT2=KWRIT1
        IF(MODCRT.LE.1) THEN
          KWRIT3=999999
        ELSE
          KWRIT3=NELEM(IKODE)
        END IF
C
C       Reading strings composing the names of output files
        DO 41 I=JWAVE+1,IWAVE-1
          READ(LUN,*,END=98) FN1C,FN2C,FN3C,(FN2D(J),J=1,MF)
   41   CONTINUE
        FN1C=' '
        FN2C=' '
        FN3C=' '
        DO 42 I=1,MF
          FN2D(I)=' '
   42   CONTINUE
        READ(LUN,*,END=98) FN1C,FN2C,FN3C,(FN2D(I),I=1,MF)
C
C       Opening the output files (CRT.5.5.1)
        IF(NF1.NE.0.) THEN
          NAME1(1)=FN1A
          NAME1(2)=FN1C
          CALL FNAME(2,NAME1,NAME2(1))
          IF(NAME2(1).NE.' ') THEN
            OPEN(LUW(1),FILE=NAME2(1),FORM='UNFORMATTED',IOSTAT=IE1)
            IF(IE1.NE.0) THEN
C             551.1
C             Open file error
              IE2=1
              J=1
              GO TO 91
            END IF
          END IF
          NAME1(1)=FN1I
          CALL FNAME(2,NAME1,NAME2(2))
          IF(NAME2(2).NE.' ') THEN
            OPEN(LUW(2),FILE=NAME2(2),FORM='UNFORMATTED',IOSTAT=IE1)
            IF(IE1.NE.0) THEN
C             551.4
C             Open file error
              IE2=4
              J=2
              GO TO 91
            END IF
          END IF
        END IF
C
C       Opening the output files (CRT.5.5.2)
        DO 44 I=1,NF2
          NAME1(1)=FN2A
          NAME1(2)=FNB(I)
          NAME1(3)=FN2C
          NAME1(4)=FN2D(I)
          IF(NAME1(2).EQ.' '.AND.NAME1(3).EQ.' '.AND.NAME1(4).EQ.' ')
     *                                                              THEN
            NAME1(1)=' '
          END IF
          J=2*(NF1+I)-1
          CALL FNAME(4,NAME1,NAME2(J))
          IF(NAME2(J).NE.' ') THEN
            OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1)
            IF(IE1.NE.0) THEN
C             551.2
C             Open file error
              IE2=2
              GO TO 91
            END IF
          END IF
          NAME1(1)=FN2I
          IF(NAME1(2).EQ.' '.AND.NAME1(3).EQ.' '.AND.NAME1(4).EQ.' ')
     *                                                              THEN
            NAME1(1)=' '
          END IF
          J=2*(NF1+I)
          CALL FNAME(4,NAME1,NAME2(J))
          IF(NAME2(J).NE.' ') THEN
            OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1)
            IF(IE1.NE.0) THEN
C             551.5
C             Open file error
              IE2=5
              GO TO 91
            END IF
          END IF
   44   CONTINUE
C
C       Opening the output files (CRT.5.5.3)
        NF3=0
        IF(MODCRT.GE.2) THEN
          IF(FN3C.EQ.' ') THEN
C           557
            CALL ERROR('557 in WRIT1: Few filenames specified')
C           There are more required elementary waves than the
C           corresponding specified filenames.
          ELSE
            NKODE=NELEM(LCODE())
C           NKODE is the number of ray elements
            IF(MODCRT.GE.3) NKODE=MIN0(KWRIT3+1,NKODE)
C           NKODE is the index of the last ray element to be stored
            NF3=NKODE-KWRIT3
            IF(NF3.GT.MF3) THEN
C             556
              CALL ERROR('556 in WRIT1: Few filenames specified')
C             There are more elements of the rays of an elementary wave
C             than the corresponding specified filenames.
            ELSE IF(2*(NF1+NF2+NF3).GT.NUW) THEN
C             554
              CALL ERROR('554 in WRIT1: Few logical units available')
C             There are NUW logical units available, see the common
C             blocks /WRITC/ and /WRITN/ defined in the block data
C             WRITB.  This number should be increased.
            END IF
            DO 45 I=1,NF3
              NAME1(1)=FN3A
              NAME1(2)=FNB(KWRIT3+I)
              NAME1(3)=FN3C
              J=2*(NF1+NF2+I)-1
              CALL FNAME(3,NAME1,NAME2(J))
              IF(NAME2(J).NE.' ') THEN
                OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1)
                IF(IE1.NE.0) THEN
C                 
C                 551.3
C                 Open file error
                  IE2=3
                  GO TO 91
                END IF
              END IF
              NAME1(1)=FN3I
              J=2*(NF1+NF2+I)
              CALL FNAME(3,NAME1,NAME2(J))
              IF(NAME2(J).NE.' ') THEN
                OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1)
                IF(IE1.NE.0) THEN
C                 
C                 551.6
C                 Open file error
                  IE2=6
                  GO TO 91
                END IF
              END IF
   45       CONTINUE
          END IF
        END IF
C
C       Opening the output file with homogeneous triangles:
        J=2*(NF1+NF2+NF3)+1
        NAME2(J)=' '
        IF(FN1T.NE.' ') THEN
          NAME1(1)=FN1T
          NAME1(2)=FN1C
          CALL FNAME(2,NAME1,NAME2(J))
C         IF(NAME2(J).NE.' ') THEN
            OPEN(LUW(J),FILE=NAME2(J),FORM='FORMATTED',IOSTAT=IE1)
            IF(IE1.NE.0) THEN
C             551.7
C             Open file error
              IE2=7
              GO TO 91
            END IF
C         end if
        END IF
C
C       Initial values for the elementary wave
        JRAY=0
        JUEB=0
        JFCT=0
        JOUTP=0
        JTRANS=0
        DO 51 I=1,MENDR+2
          NENDR(I)=0
   51   CONTINUE
        DO 52 I=1,2*(NF1+NF2+NF3)
          JPOINT(I)=0
   52   CONTINUE
C
C       Writing to the output LOG file
        WRITE(LULOG,81) IWAVE
   81   FORMAT(' WAVE',I5,':')
C
      END IF
C
C     Index of the last wave:
      JWAVE=IWAVE
C
C     Screen output
      CALL SCRO1(IWAVE)
      RETURN
C
   91 CONTINUE
C       551
        WRITE(*,'('' STATUS'',I6,''.'',I1,'': '',A)') IE1,IE2,NAME2(J)
        CALL ERROR('551 in WRIT1: Open file error')
C       Open file error when opening the file to store:
C       1 computed quantities along the rays.
C         Location of error 551.1 in the code
C       2 computed quantities at the specified surfaces.
C         Location of error 551.2 in the code
C       3 computed quantities at the endpoints of the elements of
C         the rays of the elementary wave.
C         Location of error 551.3 in the code
C       4 quantities at the initial points of rays, corresponding
C         to the above file 1.
C         Location of error 551.4 in the code
C       5 quantities at the initial points of rays, corresponding
C         to the above file 2.
C         Location of error 551.5 in the code
C       6 quantities at the initial points of rays, corresponding
C         to the above file 3.
C         Location of error 551.6 in the code
C       7 homogeneous triangles in the ray-parameter domain.
C         Location of error 551.7 in the code
C       The type 1 to 7 of the file is given by the first decimal
C       place of the status.
   98 CONTINUE
C       558
        CALL ERROR('558 in WRIT1: End of input file')
C       End of input data file WRIT specifying the names of the
C       output files encountered before all required elementary
C       waves are computed.  The number of elementary waves
C       exceeds the number of specified output filenames.
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT2(LULOG,IRAY)
      INTEGER LULOG,IRAY
C
C This subroutine is called when starting the complete tracing of a new
C ray.
C
C Input:
C     LULOG...Logical unit number of the log output device.
C     IRAY... The index of the ray which will be computed (i.e. the
C             output of the subroutine RPAR2 from the file 'rpar.for').
C
C No output.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C JRAY is redefined in this subroutine.  No other of the storage
C locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL SCRO2
C     SCRO2... File 'scro.for'.
C
C Date: 1996, June 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
C
      NTMP=0
      JRAY=IRAY
      DO 1 I=1,NF2
        JSTOR(I)=0
    1 CONTINUE
C
C     Screen output
      CALL SCRO2(IRAY)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT31(YL,Y,YY,IY)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12)
C
C This subroutine stores the computed quantities along a ray (see
C C.R.T.5.5.1).  It is called with constant step STORE of the
C independent variable along the ray, and at the points of intersection
C with interfaces either before and after the transformation.
C
C Input:
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.
C None of the input parameters are altered.
C
C No output.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C Array JPOINT is modified in this subroutine.
C
C Subroutines and external functions required:
      EXTERNAL NELEM,SCRO3
      INTEGER NELEM
C     NELEM.. File 'code.for'.
C     SCRO3.. File 'scro.for'.
C
C Date: 1996, June 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
C
      IF(NELEM(IY(2)).GT.KWRIT1) THEN
        IF(NAME2(1).NE.' ') THEN
          IF(K2P.EQ.0) THEN
            JPOINT(1)=JPOINT(1)+1
            WRITE(LUW(1))
     *            JWAVE,JRAY,IY(1),IY(5),IY(6),YY(1),YL,(Y(I),I=1,IY(1))
          ELSE
            NTMP=NTMP+1
            IF(NTMP.LE.MTMP) THEN
              KTMP(1,NTMP)=JWAVE
              KTMP(2,NTMP)=JRAY
              KTMP(3,NTMP)=IY(1)
              KTMP(4,NTMP)=IY(5)
              KTMP(5,NTMP)=IY(6)
              TMP(1,NTMP)=YY(1)
              DO 11 I=1,6
                TMP(1+I,NTMP)=YL(I)
   11         CONTINUE
              DO 12 I=1,IY(1)
                TMP(7+I,NTMP)=Y(I)
   12         CONTINUE
            ELSE
C             552
              CALL ERROR('552 in WRIT31: Too many points along a ray')
C             Dimension MTMP of array to temporarily store the traced
C             ray in the memory is smaller than the number of points
C             along the ray.
            END IF
          END IF
        END IF
      END IF
C
C     Screen output
      CALL SCRO3(YL,Y,YY,IY)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT32(ISTOR,YL,Y,YY,IY,NAMPL,YC)
      INTEGER ISTOR,IY(12),NAMPL
      REAL YL(6),Y(27),YY(5),YC(NAMPL)
C
C This subroutine stores the computed quantities at specified surfaces
C (see C.R.T.5.5.2).  It is called at the points of intersection of the
C ray with specified surfaces.
C
C Input:
C     ISTOR...The sequential number in the input data of the specified
C             surface.
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C             Quantities Y(28) to Y(IY(1)) representing reduced
C             amplitudes are ignored.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.
C     NAMPL...Number of real quantities representing complex-valued
C             vectorial reduced amplitudes.  If no conversion
C             coefficients are applied NAMPL=IY(1)-27,  otherwise
C             NAMPL=6 or 12 (see C.R.T.5.5.4).
C     YC...   Array containing real quantities representing complex-
C             -valued vectorial reduced amplitudes.  If no conversion
C             coefficients are applied, YC is a copy of Y(28) to
C             Y(IY(1)).  Otherwise,  YC represents the vectorial reduced
C             amplitudes involving appropriate conversion coefficients,
C             expressed in ray-centred coordinate system (see
C             C.R.T.5.5.4):
C             P wave at the initial point of the ray:
C               NAMPL=6,
C               YC(1)=REAL(A13),  YC(2)=AIMAG(A13),
C               YC(3)=REAL(A23),  YC(4)=AIMAG(A23),
C               YC(5)=REAL(A33),  YC(6)=AIMAG(A33).
C             S wave at the initial point of the ray:
C               NAMPL=12,
C               YC(1)=REAL(A11),  YC(2)=AIMAG(A11),
C               YC(3)=REAL(A21),  YC(4)=AIMAG(A21),
C               YC(5)=REAL(A31),  YC(6)=AIMAG(A31),
C               YC(7)=REAL(A12),  YC(8)=AIMAG(A12),
C               YC(9)=REAL(A22),  YC(10)=AIMAG(A22),
C               YC(11)=REAL(A32), YC(12)=AIMAG(A32).
C None of the input parameters are altered.
C
C No output.
C
C Common block /DCRT/ (see subroutine file 'ray.for'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C Array JPOINT is modified in this subroutine.
C
C Subroutines and external functions required:
      EXTERNAL NELEM
      INTEGER NELEM
C     NELEM...File 'code.for'.
C
C Date: 1996, June 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER NY
C
      JSTOR(ISTOR)=JSTOR(ISTOR)+1
      IF(KSTOR(NSTOR+ISTOR).LE.JSTOR(ISTOR)
     *    .AND.JSTOR(ISTOR).LE.KSTOR(2*NSTOR+ISTOR)) THEN
        IF(NELEM(IY(2)).GT.KWRIT2) THEN
          IF(NAME2(2*(NF1+ISTOR)-1).NE.' ') THEN
            JPOINT(2*(NF1+ISTOR)-1)=JPOINT(2*(NF1+ISTOR)-1)+1
            NY=27+NAMPL
            WRITE(LUW(2*(NF1+ISTOR)-1))
     *                           JWAVE,JRAY,NY,IY(5),IY(6),YY(1),YL,Y,YC
          END IF
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT33(YL,Y,YY,IY)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12)
C
C This subroutine stores the computed quantities at the endpoints of the
C individual elementary waves (see C.R.T.5.5.3).  It is called at the
C endpoints of the elements of rays, situated at structural interfaces.
C
C Input:
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.
C None of the input parameters are altered.
C
C No output.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL NELEM
      INTEGER NELEM
C     NELEM...File 'code.for'.
C
C Date: 1996, June 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage location:
      INTEGER I
C
      I=NELEM(IY(2))-KWRIT3
      IF(1.LE.I.AND.I.LE.NF3) THEN
        IF(NAME2(2*(NF1+NF2+I)-1).NE.' ') THEN
          JPOINT(2*(NF1+NF2+I)-1)=JPOINT(2*(NF1+NF2+I)-1)+1
          WRITE(LUW(2*(NF1+NF2+I)-1)) JRAY,IY,YL,(Y(I),I=1,IY(1)),YY
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT4(LULOG,IRAY,YL,Y,YY,IY,IEND,ISHEET,IREC)
C
      INTEGER LULOG,IRAY,IY(12),IEND,ISHEET,IREC
      REAL YL(6),Y(35),YY(5)
C
C This subroutine stores the quantities at the initial point of the ray
C (SEE C.R.T.6.1).  It is called after termination of tracing the ray.
C
C Input:
C     LULOG...Logical unit number of the log output device.
C     IRAY... The index of the ray which has been computed (i.e. the
C             output of the subroutine RPAR2 from the file 'rpar.for').
C     YL...   Array containing local quantities at the point of the ray.
C     Y...    Array containing basic quantities computed along the ray.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray.
C     IEND... Reason of the termination of the computation of a ray (see
C             C.R.T.5.4).  For a detailed description see subroutines
C             RAY2 (file 'ray.for') and INIT2 (file 'init.for').
C     ISHEET..Ray-history index.  The different ray histories are
C             consecutively indexed by positive integers 1,2,3,...
C             According to their appearance during ray tracing.
C             The ray histories are indexed independently within each
C             elementary wave.
C             The ray-history indices are complemented with sign:
C             positive - successful ray (crossing reference surface),
C             negative - unsuccessful ray (terminating before crossing
C             reference surface).
C     IREC... Index of the receiver for a two-point ray,
C             otherwise zero.
C
C No output.
C
C Common block /DCRT/ (see subroutine file 'ray.for'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common block /INITC/ (see subroutine file 'init.for'):
      INCLUDE 'initc.inc'
C     initc.inc
C None of the storage locations are altered.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C JUEB, JFCT, JOUTP, JTRANS, arrays JPOINT and NENDR are modified in
C this subroutine.
C
C Subroutines and external functions required:
      EXTERNAL SCRO4,CHECK
C     SCRO4... File 'scro.for'.
C     CHECK... This file.
C
C Date: 1997, September 27
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,I0,I1,I2,I3,I4,I5,I6,I7,I8,I9
      REAL S12,S13,S23,S14,S24,S34
C
C.......................................................................
C
C     Screen output
      CALL SCRO4(IRAY,YL,Y,YY,IY,IEND,ISHEET)
C
C     Statistics
      DO 11 I=1,MENDR
        IF(IEND.EQ.JENDR(I)) THEN
          NENDR(I)=NENDR(I)+1
          GO TO 12
        END IF
   11 CONTINUE
      NENDR(MENDR+1)=NENDR(MENDR+1)+1
   12 CONTINUE
      NENDR(MENDR+2)=NENDR(MENDR+2)+1
C
C     Check for non-existing rays:
      IF(IEND.GE.70) THEN
        RETURN
      END IF
C
C     Statistics
      JFCT=JFCT+IY(9)
      JOUTP=JOUTP+IY(10)
      JTRANS=JTRANS+IY(11)
C
C     Writing the quantities at the initial point of the ray (C.R.T.6.1)
      I1=-JWAVE
      DO 20 I=1,NF1+NF2+NF3
        IF(NAME2(2*I).NE.' ') THEN
C         For whole rays, only initial points of two-point rays are
C         stored if only two-point rays are stored.
C         For intersections with surfaces and endpoints of elementary
C         waves, initial points of all rays are stored.
          IF(K2P.EQ.0.OR.IREC.GT.0.OR.I.GT.NF1) THEN
            JPOINT(2*I)=JPOINT(2*I)+1
            WRITE(LUW(2*I)) I1,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI
          END IF
        END IF
   20 CONTINUE
C
C     Writing the quantities along the two-point ray:
      IF(NAME2(1).NE.' ') THEN
        IF(K2P.NE.0.AND.IREC.GT.0) THEN
          JPOINT(1)=JPOINT(1)+NTMP
          DO 30 J=1,NTMP
            WRITE(LUW(1)) (KTMP(I,J),I=1,5),(TMP(I,J),I=1,7+KTMP(3,J))
   30     CONTINUE
        END IF
      END IF
C
C     Writing to the output LOG file
C     Simplecticity of the propagator matrix (see C.R.T.5.6l)
      IF(UEBDRT.GT.0.) THEN
        S12=Y(12)*Y(18)+Y(13)*Y(19)-Y(14)*Y(16)-Y(15)*Y(17)
        S13=Y(12)*Y(22)+Y(13)*Y(23)-Y(14)*Y(20)-Y(15)*Y(21)-1.
        S23=Y(16)*Y(22)+Y(17)*Y(23)-Y(18)*Y(20)-Y(19)*Y(21)
        S14=Y(12)*Y(26)+Y(13)*Y(27)-Y(14)*Y(24)-Y(15)*Y(25)
        S24=Y(16)*Y(26)+Y(17)*Y(27)-Y(18)*Y(24)-Y(19)*Y(25)-1.
        S34=Y(20)*Y(26)+Y(21)*Y(27)-Y(22)*Y(24)-Y(23)*Y(25)
      END IF
C     Checking exceeding the specified limits
      I=0
      CALL CHECK(YY(2),UEB,I0,I)
      CALL CHECK(YY(3),UEBPP,I1,I)
      CALL CHECK(YY(4),UEBPH,I2,I)
      CALL CHECK(YY(5),UEBHH,I3,I)
      CALL CHECK(S12,UEBDRT,I4,I)
      CALL CHECK(S13,UEBDRT,I5,I)
      CALL CHECK(S23,UEBDRT,I6,I)
      CALL CHECK(S14,UEBDRT,I7,I)
      CALL CHECK(S24,UEBDRT,I8,I)
      CALL CHECK(S34,UEBDRT,I9,I)
C     I is the number of checked quantities exceeding their specified
C     Upper limits
      IF(I.GT.0) THEN
C       Writing report on this ray
        IF(JUEB.EQ.0) THEN
          WRITE(LULOG,'(2A)') '       RAY:     ',
     *        'CHECKED QUANTITIES IN PERCENTS OF THEIR SPECIFIED LIMITS'
        END IF
        JUEB=JUEB+1
        WRITE(LULOG,'(I10,A,10I6)')
     *                           IRAY,': ',I0,I1,I2,I3,I4,I5,I6,I7,I8,I9
      END IF
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C
C     
C
      SUBROUTINE CHECK(Q,QUEB,IRATE,I)
      REAL Q,QUEB
      INTEGER IRATE,I
C
C     Auxiliary subroutine to WRIT4
C
      IF(QUEB.LE.0.) THEN
        IRATE=0
      ELSE
        IF(ABS(Q).GT.QUEB) THEN
          I=I+1
        END IF
        IF(ABS(Q).GT.10000.*QUEB) THEN
          IRATE=1000000.
        ELSE
          IRATE=INT(100.*Q/QUEB+0.5)
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRIT5(LULOG,IWAVE)
      INTEGER LULOG,IWAVE
C
C This subroutine is called after termination of the computation of an
C elementary wave, and when terminating the complete ray tracing
C program.
C
C Input:
C     LULOG...Logical unit number of the log output device.
C     IWAVE...Zero when terminating the complete ray tracing program,
C             otherwise the index of the elementary wave which has been
C             computed (i.e. the output of the subroutine CODE1 from the
C             file 'code.for').
C
C No output.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL SCRO5
C     SCRO5... File 'scro.for'.
C
C Date: 1992, December 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,J,N,K(4)
C
C     Writing to the output LOG file:
      IF(IWAVE.NE.0) THEN
        IF(JUEB.GT.0) THEN
          WRITE(LULOG,10)
   10     FORMAT(11X,'ABOVE RAYS ARE COMPUTED WITH DECREASED ACCURACY')
        END IF
C       Reasons of termination of tracing rays
        WRITE(LULOG,21) NENDR(MENDR+2),JFCT,JOUTP,JTRANS
        N=0
        DO 12 J=1,MENDR
          IF(NENDR(J).GT.0) THEN
            N=N+1
            K(N)=J
          END IF
          IF(N.EQ.4.OR.(J.EQ.MENDR.AND.N.GT.0)) THEN
            WRITE(LULOG,22) (NENDR(K(I)),JENDR(K(I)),I=1,N)
            N=0
          END IF
   12   CONTINUE
        IF(NENDR(MENDR+1).NE.0) THEN
          WRITE(LULOG,23) NENDR(MENDR+1)
        END IF
C       List of output files
        DO 14 I=1,2*(NF1+NF2+NF3)
          IF(NAME2(I).NE.' ') THEN
            WRITE(LULOG,24) JPOINT(I),NAME2(I)
          END IF
   14   CONTINUE
C       formats
   21   FORMAT(  I10,' RAYS ',I12,' FCT  ',I12,' STEPS',I12,' TRANS')
   22   FORMAT(4(I10,' IEND=',I2))
   23   FORMAT(  I10,' IEND=OTHER')
   24   FORMAT(  I10,' POINTS IN FILE: ',A)
      END IF
C
C     Screen output:
      CALL SCRO5(IWAVE)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FNAME(NUM1,NAME1,NAME2)
      INTEGER NUM1
      CHARACTER*(*) NAME1(NUM1)
      CHARACTER*(*) NAME2
C
C This subroutine is designed to analyze the specified strings and to
C compose the filename of them.
C
C Input:
C     NAME1...Character array of character strings to be analyzed.
C     NUM1... Number of filenames to be analyzed.
C
C Output:
C     NAME2...Filename composed of the analyzed input components
C             (strings).
C
C No subroutines and external functions required.
C
C Date: 1992, December 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER M1
      PARAMETER (M1=4)
      INTEGER J1(M1),L1(M1),I1,J2,L2,I2,N,I
C
C     M1...   Maximum number of analyzed strings.
C     J1(N)...Number of the currently analyzed characters of NAME1(N),
C             i.e. the position of the space after the last analyzed
C             word of NAME1(N).  A word is a substring between two
C             consecutive spaces.
C     L1(N)...Length of NAME1(N).
C     I1...   Beginning of the last analyzed word of NAME1(N), i.e. the
C             previous value of J1(N) increased by 1.
C     J2...   Number of the currently defined characters of NAME2.
C     L2...   Length of NAME2.
C     I2...   Auxiliary variable - previous value of J2 increased by 1.
C     N...    Loop variable - index of input filename.
C     I...    Loop variable - position in a string.
C
C.......................................................................
C
      IF(M1.LT.NUM1) THEN
C       559
        CALL ERROR('559 in FNAME: Too many strings')
C       The number NUM1 of input strings exceeds the dimension M1
C       of arrays J1 and L1.  Parameter M1 must be increased.
      END IF
      DO 11 N=1,NUM1
        J1(N)=0
        L1(N)=LEN(NAME1(N))
   11 CONTINUE
      J2=0
      L2=LEN(NAME2)
C
   20 CONTINUE
      DO 29 N=1,NUM1
C       (a) Analyzing NAME1(N):
        I1=J1(N)+1
        DO 21 I=I1,L1(N)
          IF(NAME1(N)(I:I).EQ.' ') THEN
            J1(N)=I
            GO TO 22
          END IF
   21   CONTINUE
        J1(N)=L1(N)+1
   22   CONTINUE
C       J1(N)-th character of NAME1(N) is found to be blank.
C       (b) Copying the word from NAME1(N) to NAME2:
        IF(J1(N).GT.I1) THEN
          I2=J2+1
          J2=MIN0(J2+J1(N)-I1,L2)
          NAME2(I2:J2)=NAME1(N)(I1:J1(N)-1)
        END IF
   29 CONTINUE
C
      DO 31 N=1,NUM1
        IF(J1(N).LT.L1(N)) GO TO 20
   31 CONTINUE
      DO 32 I=J2+1,L2
        NAME2(I:I)=' '
   32 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WRITA(IRAY1,IRAY2,IRAY3)
      INTEGER IRAY1,IRAY2,IRAY3
C
C This subroutine is designed to write the indices of rays forming
C the homogeneous triangles in the ray-parameter domain.
C
C Input:
C     IRAY1,IRAY2,IRAY3... Indices of rays forming the homogeneous
C             triangle.  IRAY3=0 in 2-D.
C
C No output.
C
C Common blocks /WRITC/ and /WRITN/:
      INCLUDE 'writ.inc'
C     writ.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1997, October 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IFORM1,IFORM3
      CHARACTER*10 FORMAT
      SAVE IFORM1,IFORM3,FORMAT
      DATA IFORM1/99999/,IFORM3/0/,FORMAT/'(I6,I6,I2)'/
C
C     Writing to the output file with homogeneous triangles:
      IF(NAME2(2*(NF1+NF2+NF3)+1).NE.' ') THEN
        IF(IRAY1.GT.IFORM1.OR.IRAY2.GT.IFORM1.OR.IRAY3.GT.IFORM1) THEN
          IFORM1=IFORM1*10+9
          FORMAT(3:3)=CHAR(ICHAR(FORMAT(3:3))+1)
          FORMAT(6:6)=FORMAT(3:3)
          IFORM3=0
        END IF
        IF(IRAY3.GT.IFORM3) THEN
          IFORM3=IFORM1
          FORMAT(9:9)=FORMAT(3:3)
        END IF
        WRITE(LUW(2*(NF1+NF2+NF3)+1),FORMAT) IRAY1,IRAY2,IRAY3
      END IF
C
      RETURN
      END
C
C=======================================================================
C