C
C Subroutine file 'writ.for' to create the output files of complete ray
C tracing
C
C Date: 2005, June 20
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             see
C             C.R.T.5.5.1.
C             It is called with constant step STORE (see input data
C             in the file 'ray.for') of the independent variable along
C             the ray, and at the points of intersection with interfaces
C             either before and after the transformation.
C             WRIT31
C     WRIT32..Subroutine storing the computed quantities at specified
C             surfaces, see
C             C.R.T.5.5.2.
C             It is called at the points of intersection of the ray
C             with the specified surfaces.
C             WRIT32
C     WRIT33..Subroutine storing the computed quantities at the
C             endpoints of the individual elementary waves, see
C             C.R.T.5.5.3.
C             It is called at the endpoints of the elements of rays,
C             situated at structural interfaces.
C             WRIT33
C     WRIT4...Subroutine storing the quantities at the initial point of
C             the ray, see
C             C.R.T.6.1.
C             It is called after termination of 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     WRITTR..Subroutine writing the indices of rays forming the
C             homogeneous triangles in the ray-parameter domain.
C             WRITTR
C     WRITBR...Subroutine writing the indices of boundary rays.
C             WRITBR
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).
C     The strings have thus to be enclosed in apostrophes.  In the list
C     of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new READ statement).  All
C     input variables, except K2P, are of the type CHARACTER.
C     The data are read by 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 read.
C     Default: TEXTW=' '.
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     /...    An obligatory slash.
C     Default:  K2P=0.
C (3) 'FN1A','FN1I','FN1T','FN1BR',/
C     'FN1A'... String containing the first component of the name of the
C             output file 'CRT-R' with the computed quantities stored
C             along the rays, see
C             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 'CRT-I' with the quantities at the initial
C             points of rays, see
C             C.R.T.6.1,
C             related to the output file 'CRT-R' with the computed
C             quantities stored along the rays.  The other components
C             of the filename are common with the related file 'CRT-R'
C             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 'CRT-T' specifying the homogeneous
C             triangles in the ray-parameter domain.  The other
C             components of the filename are common with the related
C             file 'CRT-R' with the computed quantities along rays.
C             Description of file with triangles
C     'FN1BR'... String containing the first component of the name of
C             the formatted output file 'CRT-B' specifying the indices
C             of boundary rays.  The other components of the filename
C             are common with the related file 'CRT-R' containing
C             the computed quantities along rays.
C            Description of file with boundary rays
C     /...    An obligatory slash.
C     Default: all strings 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 'CRT-S' with the computed quantities stored at
C             the specified surfaces, see
C             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 'CRT-I' with the quantities at the initial
C             points of rays, see
C             C.R.T.6.1,
C             related to the output file 'CRT-S' with the computed
C             quantities stored at the specified surfaces.  The other
C             components of the filename are common with the related
C             file 'CRT-S' containing the computed 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 'CRT-S' with the computed quantities
C             stored at the I-th specified surface.  Here NENDF is
C             the number of surfaces specified for storing the computed
C             quantities.
C     /...    An obligatory slash.
C     Default: all strings blank.
C (5) Either (5.1) or (5.2):
C (5.1) /
C     A slash is fully sufficient if MODCRT.LE.1,
C     see input data set DCRT.
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 'CRT-E' with the computed quantities stored
C             at the endpoints of the elements of rays, see
C             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 'CRT-I' with the quantities at the initial
C             points of rays, see
C             C.R.T.6.1,
C             related to the output file 'CRT-E' with the computed
C             quantities stored at the endpoints of the elements
C             of rays.  The other components of the filename are common
C             with the related file 'CRT-E' containing the computed
C             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 'CRT-E' with the computed quantities
C             stored at the endpoints of the I-th element of rays.
C             None of the strings FN3B(I), corresponding to the endpoint
C             of a ray element at which quantities are stored, may equal
C             spaces for MODCRT.GE.2.
C     /...    An obligatory slash.
C     Default: all strings blank.
C (6) For each elementary wave IWAVE the following data (6.1):
C (6.1) 'FN1C','FN2C','FN3C',('FN2D(I)',I=1,NENDF),/
C     'FN1C'... String containing the second component of the name of
C             the output file 'CRT-R' with the computed quantities
C             stored along the rays, see
C             C.R.T.5.5.1.
C     'FN2C'... String containing the third component of the name of the
C             output file 'CRT-S' with the computed quantities stored at
C             the specified surfaces, see
C             C.R.T.5.5.2.
C     'FN3C'... String containing the third component of the name of the
C             output file 'CRT-E' with the computed quantities stored at
C             the endpoints of the elements of rays, see
C             C.R.T.5.5.3.
C             This string may not equal spaces for MODCRT.GE.2.
C             If MODCRT.LE.1, the file is not created.
C     'FN2D(I)'... String containing the fourth component of the name of
C             the output file 'CRT-S' with the computed quantities
C             stored at the I-th specified surface, see
C             C.R.T.5.5.2.
C             Here NENDF is the number of surfaces specified for storing
C             the computed quantities.
C     /...    An obligatory slash.
C     Default: all strings blank.
C Following table summarizes the components, of which the individual
C filenames of the output files are composed:
C       Output file:                Components:
C     'CRT-R'                   FN1A +            FN1C(IW)
C     corresponding 'CRT-I'     FN1I +            FN1C(IW)
C     corresponding 'CRT-T'     FN1T +            FN1C(IW)
C     corresponding 'CRT-B'     FN1BR+            FN1C(IW)
C     'CRT-S'                   FN2A + FN2B(IS) + FN2C(IW) + FN2D(IS,IW)
C     corresponding 'CRT-I'     FN2I + FN2B(IS) + FN2C(IW) + FN2D(IS,IW)
C     'CRT-E'                   FN3A + FN3B(IE) + FN3C(IW)
C     corresponding 'CRT-I'     FN3I + FN3B(IE) + FN3C(IW)
C   IW is a total index of elementary wave, IW=NWAVE*(ISRC-1)+IWAVE,
C   where ISRC is index of the source, IWAVE is the index of
C   the elementary wave for source ISRC, and NWAVE is the number of
C   elementary waves to be calculated for each source.
C   IS is an index of a storing surface, IS=1,NENDF.
C   IE is an index of a storing element, IE=1,NELEM.
C Notes:
C   If the first component of the name of any of the CRT output
C     files is blank, the file is not created.
C   File 'CRT-S' and corresponding file 'CRT-I' are not created if
C     FN2B(IS), FN2C(IW) and FN2D(IS,IW) are all blank.
C   The components FN1A, FN1I, FN1T and FN1BR must be mutually
C     different, if they are not blank.
C   The components FN2A and FN2I must differ, if they are not blank.
C   The components FN3A and FN3I must differ, if they are not blank.
C
C Examples of data set WRIT:
C     writ.dat    two-point rays stored.
C     writall.dat all rays stored.
C     writsrf.dat no rays stored.
C In the above examples, the computed quantities at specified surfaces
C are stored.
C
C The names of the output files of program 'crt.for' are specified here
C in the input data WRIT.  On the other hand, the programs processing
C the results of complete ray tracing assume the names of the same
C output files of program 'crt.for' to be specified in file CRTOUT
C formatted as follows:
C                                                  
C General description of formatted file CRTOUT:
C     The file consists of character strings to be read by list directed
C     (free format) input.  The strings have thus to be enclosed in
C     apostrophes.
C (1) For each set of CRT output files data record (1.1):
C (1.1) 'CRT-R','CRT-S','CRT-I','CRT-T','CRT-B',/
C     'CRT-R'.. Name of the unformatted CRT output file with the
C             quantities stored along rays, see
C             C.R.T.5.5.1.
C             Description of file 'CRT-R'
C     'CRT-S'.. Name of the unformatted CRT output file with the
C             quantities stored at the points of intersections of rays
C             with the specified surfaces, see
C             C.R.T.5.5.2.
C             Description of file 'CRT-S'
C     'CRT-I'.. Name of the unformatted CRT output file with the
C             quantities at the initial points of rays, see
C             C.R.T.6.1.
C             Description of file 'CRT-I'
C     'CRT-T'.. Name of the formatted CRT output file with the
C             homogeneous triangles in the ray-parameter domain.
C             Description of file 'CRT-T'
C     'CRT-B'.. Name of the formatted CRT output file with the
C             indices of boundary rays.
C             Description of file 'CRT-B'
C     /...    An obligatory slash.
C     General defaults for the first set of the CRT output files:
C             The defaults may be blank or following:
C             'CRT-R'='r01.out', 'CRT-S'='s01.out', 'CRT-I'='r01i.out'
C             or 'CRT-I'='s01i.out', 'CRT-T'='t01.out', 'CRT-B'=' '
C     General defaults for the subsequent sets:
C             'CRT-R'=' ', 'CRT-S'=' ', 'CRT-I'=' ', 'CRT-T'=' ',
C             'CRT-B'=' '
C (2) / (a slash or the end of file).
C     The end of the data is also assumed if all filenames are blank.
C
C                                                    
C Output files 'CRT-R' with the computed quantities stored along the
C rays, see 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) ISRC,IWAVE,IRAY,NY,ICB1,ISRF,X,(YL(I),I=1,6),(Y(I),I=1,NY)
C     ISRC .. Index of the source.
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             see
C             C.R.T.5.5.4.
C             Description of YL
C     Y...    Array containing basic quantities computed along the ray,
C             see
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, see
C 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) ISRC,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     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, see
C 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             see
C             C.R.T.5.5.4.
C             Description of YL
C     Y...    Array containing basic quantities computed along the ray,
C             see
C             C.R.T.5.2.1.
C             Description of Y
C     YY...   Array containing real auxiliary quantities computed along
C             the ray, see
C             C.R.T.5.2.2.
C             Description of YY
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray, see
C             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
C rays, see C.R.T.6.1.
C     Unformatted output.  The quantities are stored at the initial
C             points of rays.
C     For whole rays ('CRT-R'), only initial points of two-point
C             rays are stored if only two-point rays are stored.
C     For intersections with surfaces ('CRT-S') and endpoints of the
C             elements of rays ('CRT-E') initial points of all rays
C             are stored.
C     For initial point - one record containing following 34 quantities:
C (1) ISRC,-IWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,(YLI(I),I=1,6),
C     (YI(I),I=1,29)
C     ISRC .. Index of the source.
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             C.R.T.6.1.
C     IEND... Reason of the termination of the computation of a ray, see
C             C.R.T.5.4.
C             See subroutine RAY in the subroutine file 'ray.for' for
C             detailed description of IEND.
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             C.R.T.5.5.4,
C             describing the local properties of the model at the
C             initial point of the ray.   See the list of YL(1) to YL(6)
C             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             C.R.T.6.1.
C             These quantities are listed in the file 'init.for'.
C             Description of YI
C
C                                                    
C Output files 'CRT-T' specifying homogeneous triangles in the
C ray-parameter domain. Formatted output, each record containing
C three integers:
C (1) For each elementary wave IWAVE the following data (1.1) and (1.2):
C (1.1) 0,ISRC,IWAVE
C     0  ...  Zero identifying new elementary wave.
C     ISRC .. Index of the source.
C     IWAVE . Index of the elementary wave.
C (1.2) For each homogeneous triangle of the elementary wave
C     data (1.2.1):
C (1.2.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 Output file 'CRT-B' specifying indices of boundary rays.
C Formatted output, each record containing four integers:
C (1) For each elementary wave IWAVE the following data (1.1) and (1.2):
C (1.1) 0,0,ISRC,IWAVE
C     0  ...  Zero identifying new elementary wave.
C     0  ...  Zero.
C     ISRC .. Index of the source.
C     IWAVE . Index of the elementary wave.
C (1.2) For each boundary ray data (1.2.1):
C (1.2.1) IB1,IB2,IB3,IB4
C     IB1 ... Index of the boundary ray.
C     IB2 ... Index of the boundary ray coupled with IB1.
C     IB3,IB4 . Indices of two boundary rays equal in history with IB1,
C             coupled with rays equal in history with IB2, and connected
C             with boundary ray IB1 by sides of homogeneous triangles.
C             IB4=0 if there is only one ray of the described
C             properties.
C             IB3=0 and IB4=0 if there is no ray of the described
C             properties, and for one-parametric shooting in 2-D.
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,
     *         26,27,28,29,30,31,32,33,34,35,36,37,38,39/
      DATA JENDR/10,21,22,23,24,25,26,30,32,33,37,39,
     *           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...One of files 'scro*.for'.
C     FNAME... This file.
C
C Date: 2005, June 20
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       Initiating the source index
        ISRC=1
C
C       Reading the name of the input data
        TEXTW=' '
        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=' '
        FN1BR=' '
        READ(LUN,*) FN1A,FN1I,FN1T,FN1BR
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.1
          CALL ERROR('553.1 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
        IF ( ((FN1A .NE.' ').AND.
     *     ((FN1A .EQ.FN1I).OR.(FN1A .EQ.FN1T).OR.(FN1A .EQ.FN1BR))).OR.
     *       ((FN1I .NE.' ').AND.
     *     ((FN1I .EQ.FN1A).OR.(FN1I .EQ.FN1T).OR.(FN1I .EQ.FN1BR))).OR.
     *       ((FN1T .NE.' ').AND.
     *     ((FN1T .EQ.FN1A).OR.(FN1T .EQ.FN1I).OR.(FN1T .EQ.FN1BR))).OR.
     *       ((FN1BR.NE.' ').AND.
     *     ((FN1BR.EQ.FN1A).OR.(FN1BR.EQ.FN1I).OR.(FN1BR.EQ.FN1T ))).OR.
     *       ((FN2A.NE.' ').AND.(FN2A.EQ.FN2I)).OR.
     *       ((FN3A.NE.' ').AND.(FN3A.EQ.FN3I)) ) THEN
C         554
          CALL ERROR('554 in WRIT1: Equal first filename components.')
C         The components FN1A, FN1I, FN1T and FN1BR must be mutually
C         different, if they are not blank. Similarly FN2A and FN2I must
C         differ, and FN3A and FN3I must differ, if they are not blank.
        ENDIF
C
      ELSE
C
C       Updating the source index
        IF(IWAVE.LE.JWAVE) THEN
          ISRC=ISRC+1
        END IF
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
          NAME2(1)=' '
          IF(FN1A.NE.' ') THEN
            CALL FNAME(2,NAME1,NAME2(1))
            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
          NAME2(2)=' '
          IF(FN1I.NE.' ') THEN
            CALL FNAME(2,NAME1,NAME2(2))
            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
          NAME2(J)=' '
          IF(FN2A.NE.' ') THEN
            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
          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)
          NAME2(J)=' '
          IF(FN2I.NE.' ') THEN
            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
          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             553.2
              CALL ERROR('553.2 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
              NAME2(J)=' '
              IF(FN3A.NE.' ') THEN
                CALL FNAME(3,NAME1,NAME2(J))
                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)
              NAME2(J)=' '
              IF(FN3I.NE.' ') THEN
                CALL FNAME(3,NAME1,NAME2(J))
                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
        IF (J.GT.NUW) THEN
C         553.3
          CALL ERROR('553.3 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.  The number NUW should be increased.
        ENDIF
        NAME2(J)=' '
        IF(FN1T.NE.' ') THEN
          NAME1(1)=FN1T
          NAME1(2)=FN1C
          CALL FNAME(2,NAME1,NAME2(J))
          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
        END IF
C
C       Opening the output file with boundary rays:
        J=2*(NF1+NF2+NF3)+2
        IF (J.GT.NUW) THEN
C         553.4
          CALL ERROR('553.4 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.  The number NUW should be increased.
        ENDIF
        NAME2(J)=' '
        IF(FN1BR.NE.' ') THEN
          NAME1(1)=FN1BR
          NAME1(2)=FN1C
          CALL FNAME(2,NAME1,NAME2(J))
          OPEN(LUW(J),FILE=NAME2(J),FORM='FORMATTED',IOSTAT=IE1)
          IF(IE1.NE.0) THEN
C           551.8
C           Open file error
            IE2=8
            GO TO 91
          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
C       Writing to the file 'CRT-T' with homogeneous triangles:
        CALL WRITTR(0,ISRC,IWAVE)
C
C       Writing to the file 'CRT-B' with boundary rays:
        CALL WRITBR(0,0,ISRC,IWAVE)
C
      END IF
C
C     Index of the last wave:
      JWAVE=IWAVE
C
C     Screen output
      CALL SCRO1(ISRC,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       8 boundary rays.
C         Location of error 551.8 in the code
C       The type 1 to 8 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...One of files '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.
C It is called with constant step STORE of the independent variable
C along the ray, and at the points of intersection with interfaces
C 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...One of files 'scro*.for'.
C
C Date: 2003, November 24
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)) ISRC,JWAVE,JRAY,IY(1),IY(5),IY(6),IY(12),
     *                    YY(1),YY(2),YL,(Y(I),I=1,IY(1))
          ELSE
            NTMP=NTMP+1
            IF(NTMP.LE.MTMP) THEN
              KTMP(1,NTMP)=ISRC
              KTMP(2,NTMP)=JWAVE
              KTMP(3,NTMP)=JRAY
              KTMP(4,NTMP)=IY(1)
              KTMP(5,NTMP)=IY(5)
              KTMP(6,NTMP)=IY(6)
              KTMP(7,NTMP)=IY(12)
              TMP(1,NTMP)=YY(1)
              TMP(2,NTMP)=YY(2)
              DO 11 I=1,6
                TMP(2+I,NTMP)=YL(I)
   11         CONTINUE
              DO 12 I=1,IY(1)
                TMP(8+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.
C It is called at the points of intersection of the ray with specified
C 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             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: 2003, November 23
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))
     *                      ISRC,JWAVE,JRAY,NY,IY(5),IY(6),IY(12),
     *                      YY(1),YY(2),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 C.R.T.5.5.3.
C It is called at the endpoints of the elements of rays, situated
C 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.
C 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.
C             See subroutine RAY in the subroutine file 'ray.for' for
C             detailed description of IEND.
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...One of files 'scro*.for'.
C     CHECK...This file.
C
C Date: 2003, November 24
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,J
      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, see
C     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)) ISRC,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,7),(TMP(I,J),I=1,8+KTMP(4,J))
   30     CONTINUE
        END IF
      END IF
C
C     Writing to the output LOG file
C     Simplecticity of the propagator matrix, see
C     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...One of files '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 WRITTR(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
C     
C
      SUBROUTINE WRITBR(IB1,IB2,IB3,IB4)
      INTEGER IB1,IB2,IB3,IB4
C
C This subroutine is designed to write the indices of boundary rays.
C
C Input:
C     IB1 ... Index of the boundary ray.
C     IB2 ... Index of the boundary ray coupled with IB1.
C     IB3,IB4 . Indices of two boundary rays equal in history with IB1,
C             coupled with rays equal in history with IB2, and connected
C             with boundary ray IB1 by sides of homogeneous triangles.
C             IB4=0 if there is only one ray of the described
C             properties.
C             IB3=0 and IB4=0 if there is no ray of the described
C             properties and for one-parametric shooting 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: 2001, October 29
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IFORM1,IFORM3
      CHARACTER*13 FORMAT
      SAVE IFORM1,IFORM3,FORMAT
      DATA IFORM1/99999/,IFORM3/0/,FORMAT/'(I6,I6,I6,I2)'/
C
C     Writing to the output file with boundary rays:
      IF(NAME2(2*(NF1+NF2+NF3)+2).NE.' ') THEN
        IF(IB1.GT.IFORM1.OR.IB2.GT.IFORM1.OR.IB3.GT.IFORM1.OR.
     *     IB4.GT.IFORM1) THEN
          IFORM1=IFORM1*10+9
          FORMAT(3:3)=CHAR(ICHAR(FORMAT(3:3))+1)
          FORMAT(6:6)=FORMAT(3:3)
          FORMAT(9:9)=FORMAT(3:3)
          IFORM3=0
        END IF
        IF(IB4.GT.IFORM3) THEN
          IFORM3=IFORM1
          FORMAT(12:12)=FORMAT(3:3)
        END IF
        WRITE(LUW(2*(NF1+NF2+NF3)+2),FORMAT) IB1,IB2,IB3,IB4
      END IF
C
      RETURN
      END
C
C=======================================================================
C