C 
C File 'crtin.for' for reading the input data for complete ray tracing
C
C Version: 7.10
C Date: 2014, June 18
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of:
C     CRTIN...Subroutine designed to open the data files for complete
C             ray tracing and to read the input data sets CRT, MODEL,
C             DCRT and INIT.
C             CRTIN
C     UNIT... Subroutine designed to control connecting and
C             disconnecting the data files to logical units, and to
C             determine the logical units from which the given data sets
C             are to be read.  Called e.g. by the subroutine CRTIN.
C             UNIT
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Input data - main data set 'CRT':
C     This main data file for the complete ray tracing program must have
C     the form of
C     SEP (Stanford
C     Exploration Project) parameter or history file.  It contains the
C     names of other input files and the name of the output log file of
C     program CRT.
C Names of the input files:
C     MODEL='string'... String containing the name of the file with the
C             input data for the model.  The data are read by subroutine
C             MODEL1.  See the description of the data file MODEL
C             in subroutine file 'model.for'.
C             Note that it is recommended to define only surfaces
C             covering structural interfaces (Model Surfaces) in data
C             set MODEL, and to define Auxiliary Surfaces in data set
C             DCRT.
C             Default: MODEL='model.dat'
C     SRC='string'... Name of the input file containing the coordinates
C             of the single initial point and the initial value of the
C             travel time for each calculation.
C             If there are several source points in file SRC, ray
C               tracing is repeated for all these source points,
C               repeating for each source point all elementary waves
C               specified in data set
C               CODE.
C             If optional input SEP parameter RPAREW=' ' (default),
C               data in data set
C               RPAR
C               are not repeated, they should be specified for all
C               elementary waves of all source points.
C             If optional input SEP parameter RPAREW is specified and is
C               not blank, data RPAR-(4) for a single source should be
C               located in file
C               RPAREW.
C               If there are several source points, data RPAR-(4) are
C               repeatedly used for all source points.
C             If optional input SEP parameter WRITEW=' ' (default),
C               data in data set
C               WRIT
C               are not repeated, they should be specified for all
C               elementary waves of all source points.
C             If optional input SEP parameter WRITEW is specified and is
C               not blank, data WRIT-(6) for a single source should be
C               located in file
C               WRITEW.
C               If there are several source points, data WRIT-(6) are
C               repeatedly used for all source points.
C             Parameter SRC has no meaning for line or surface initial
C             conditions for rays. Description of
C             file SRC.
C             Default: SRC='src.dat'
C     REC='string'... Name of the file with profile or receiver
C             coordinates.
C             Required just for two-point ray tracing.  Description of
C             file REC.
C             Default: REC='rec.dat'
C     DCRT='string'... String containing the name of the file with the
C             input data for the complete ray tracing.  The data are
C             read by subroutine RAY1.  See the description of
C             data set DCRT
C             in subroutine file 'ray.for'.
C             Default: DCRT='crt.dat'
C     INIT='string'... String containing the name of the file with the
C             input data specifying the line or surface initial
C             conditions for rays.
C             Parameter INIT has no meaning for a point source.
C             If INIT=DCRT, data set INIT is appended to file DCRT
C             (before possible data sets CODE, INIT or WRIT).  See the
C             description of
C             data set INIT
C             in subroutine file 'init.for'.
C             The data are read by subroutine INIT1.
C             Default: INIT=' '
C     CODE='string'... String containing the name of the file with the
C             codes of elementary waves.  It is recommended to write
C             data set CODE into a separate file. See the description of
C             data set CODE
C             in subroutine file 'code.for'.
C             The data are read by subroutine CODE1.
C             Note: for single source, it is possible to specify
C             CODE=DCRT, and to append data set CODE to file DCRT.
C             It is recommended to append at most one of sets CODE,
C             RPAR, WRIT to DCRT.
C             Default: CODE='code.dat'
C                                                    
C     RPAR='string'... String containing the name of the file with the
C             data specifying the take-off parameters of the required
C             rays.
C             If RPAREW=' ' and there are several source points in file
C             SRC, data RPAR-(4) corresponding to individual elementary
C             waves should be repeated for all source points.
C             Otherwise, file RPAR should contain only data RPAR-(1) to
C             RPAR-(3) common to all elementary waves and data RPAR-(4)
C             should be located in optional file RPAREW described below.
C             If RPAR=DCRT, data set RPAR is appended to file DCRT.
C             It is recommended to append at most one of sets CODE,
C             RPAR, WRIT to DCRT.  See the description of
C             data set RPAR
C             in subroutine file 'rpar.for'.
C             Default: RPAR='rpar.dat'
C                                                  
C     RPAREW='string'... String containing the name of the optional file
C             with data RPAR-(4) specifying the take-off parameters of
C             the required rays, corresponding to individual elementary
C             waves.
C             RPAREW=' ':  All data specifying the take-off parameters
C               are located in file RPAR.  In this case file RPAR should
C               contain data RPAR-(1) to RPAR-(3) common to all
C               elementary waves followed by data RPAR-(4) corresponding
C               to individual elementary waves.  If there are several
C               source points in file SRC, the data corresponding to
C               individual elementary waves should be repeated for all
C               source points.
C             RPAREW.NE.' ':  File RPAR contains only data RPAR-(1) to
C               RPAR-(3) common to all elementary waves.
C               Data RPAR-(4) corresponding to individual elementary
C               waves from a single source are located in file RPAREW.
C               If there are several source points in file SRC, the data
C               corresponding to individual elementary waves are
C               repeatedly used for all source points.
C             Default: RPAREW=' '
C                                                    
C     WRIT='string'... String containing the name of the file specifying
C             the names of the output files with the computed
C             quantities.  These data are read by subroutine WRIT1.
C             If WRITEW=' ' and there are several source points in file
C             SRC, data WRIT-(6) corresponding to individual elementary
C             waves should be repeated for all source points.
C             Otherwise, file WRIT should contain only data WRIT-(1) to
C             WRIT-(5) common to all elementary waves and data WRIT-(6)
C             should be located in optional file WRITEW described below.
C             If WRIT=DCRT, data set WRIT is appended to file DCRT.
C             It is recommended to append at most one of sets CODE,
C             RPAR, WRIT to DCRT.  See the description of
C             data set WRIT
C             in subroutine file 'writ.for'.
C             Default: WRIT='writ.dat'
C                                                  
C     WRITEW='string'... String containing the name of the optional file
C             with data WRIT-(6) specifying the take-off parameters of
C             the required rays, corresponding to individual elementary
C             waves.
C             WRITEW=' ':  All data WRIT-(1) to WRIT-(6) are located in
C               file WRIT.  If there are several source points in file
C               SRC, data WRIT-(6) corresponding to individual
C               elementary waves should be repeated for all source
C               points.
C             WRITEW.NE.' ':  File WRIT contains only data WRIT-(1) to
C               WRIT-(5) common to all elementary waves.  Data WRIT-(6)
C               corresponding to individual elementary waves from a
C               single source are located in file WRITEW.  If there are
C               several source points in file SRC, data WRIT-(6)
C               corresponding to individual elementary waves are
C               repeatedly used for all source points.
C             Default: WRITEW=' '
C Names of the output files:
C     The names of the output files with the computed quantities, see
C             C.R.T.5.5,
C             are specified by the file given by input parameter WRIT
C             described above.
C     CRTLOG='string'... The string containing the name of the output
C             log file.  The data are written by subroutines WRIT1,
C             WRIT2, WRIT4 and WRIT5 of file
C             writ.for.
C             Unfortunately, there is no description of the output file,
C             yet.
C             Default: CRTLOG='crtlog.out'
C Optional numerical parameters to control ray tracing:
C     UEBMUL=real... Enables to overwrite the default upper error bound
C             in determining the point of intersection of a ray with
C             a structural interface, surface limiting the computational
C             volume, storing or end surface.
C             UEBMUL=0: The default upper error bound is used.
C               The default upper error bound is the greater one of the
C               following:
C               (a) Current upper error bound for numerical integration
C                   based on (and often equal to) the input value
C                   UEB
C                   of the upper error bound of travel time per one step
C                   of numerical integration (measured in travel time).
C               (b) The largest absolute value of a coordinate divided
C                   by 1000000 (measured in coordinate units).
C             UEBMUL.GT.0: The upper error bound is UEBMUL times the
C               current upper error bound (a) for numerical integration.
C               Note that option UEBMUL=1 corresponds to CRT ver. 5.70
C               or older.
C             Default: UEBMUL=0 (recommended)
C     KMAHI='integer'... Enables to exclude the KMAH index from the
C             ray-history index.  See the description
C             in file rpar.for.
C                                                     
C Numerical parameters to optionally specify anisotropic ray tracing:
C     CRTANI=real... Switch to anisotropic ray tracing:
C             CRTANI=0... Isotropic ray tracing.
C             CRTANI.NE.0... Anisotropic ray tracing.
C               Anisotropic ray tracing is now applicable only under the
C               following restrictions:
C               (1) Coordinates: Only Cartesian coordinates (raycb.for,
C                   init.for).
C               (2) Model: No structural interfaces
C                   (trans.for: reference isotropic model is now used
C                   for the transformation at the interfaces, direction
C                   of the incident slowness vector is preserved but its
C                   norm is adjusted before transformation).
C               (3) Model: No attenuation is considered with anisotropic
C                   ray tracing.
C               (4) Source: Initial point or constant travel time along
C                   an initial line or initial surface (init.for).
C             Default: CRTANI=0.
C     DEGREE=real... Degree of the homogeneous Hamiltonian to be
C             arithmetically averaged for common S-wave ray tracing.
C             Only values -2, -1, 1 and 2 are now allowed.
C             Default: DEGREE=-1.
C     KSWAVE=integer... Selection of the anisotropic-ray-theory S-wave
C             rays instead of common S-wave rays. Affects only S waves.
C             See the description in file hder.for.
C     DSWAVE=real... Minimum relative difference between the eigenvalues
C             of the S1 and S2 waves. Used only for anisotropic-ray-
C             -theory S waves.
C             See the description in file hder.for.
C                                                    
C Numerical parameters defining the kind of initial conditions (the kind
C             of the source), read by subroutine file
C             init.for:
C     INIDIM=integer... Determines the dimensionality of the source:
C             0...Single initial point (point source), its coordinates
C                 are read from file given by parameter SRC.
C             1...Initial curve (line source), see parameter INIT above.
C             2...Initial surface, see parameter INIT above.
C             Default: INIDIM=0 (point source)
C     INIE1=real, INIE2=real, INIE3=real ... Three components of the
C             optional rotation vector describing the rotation of the
C             local Cartesian basis vectors.
C             The local Cartesian basis vectors are used to define the
C             local spherical coordinates if INIDIM=0, or if INIDIM=2
C             and INIPAR is positive.
C             The contravariant (covariant) components of the local
C             Cartesian basis vectors, expressed in the curvilinear
C             global model coordinates, are given by the square root of
C             the contravariant (covariant) metric tensor.
C             The local Cartesian basis vectors thus coincide with the
C             basis vectors of the global model coordinate system if
C             the model coordinates are Cartesian.
C             The local Cartesian basis vectors may then optionally
C             be rotated about the axis given by the rotation vector.
C             The angle of the rotation in radians is given by the
C             Cartesian length of the rotation vector.
C             The rotation vector is expressed with respect to the
C             local Cartesian basis vectors before rotation.
C             Equations:
C               The i-th component of the j-th rotated local Cartesian
C               basis vector, expressed with respect to the local
C               Cartesian basis vectors before rotation, is
C                 E_ij = cos(E) delta_ij + [1-cos(E)] E_i/E E_j/E
C                      - sin(E) epsilon_ijk E_k/E
C               where (E_1,E_2,E_3)=(INIE1,INIE2,INIE3),
C                 E = sqrt(E_k E_k)
C               is the Cartesian length of the rotation vector,
C               delta_ij is the Kronecker delta, and
C               epsilon_ijk is completely skew Levi-Civita's symbol,
C               with epsilon_123=epsilon_231=epsilon_312=1.
C               The contravariant (covariant) components of the rotated
C               local Cartesian basis vectors, expressed in the
C               curvilinear global model coordinates, are obtained by
C               multiplying E_ij from the left by the square root of
C               the contravariant (covariant) metric tensor.
C             Default:
C               For positive INIPAR: INIE1=INIE2=INIE3=0. (no rotation)
C               For negative INIPAR: INIE1=pi, INIE2=INIE3=0. (rotation
C                 by pi radians (180 degrees) about the local Cartesian
C                 X1 coordinate axis.
C     INIDIR=integer... A special option for the above described
C             rotation vector:
C             INIDIR.EQ.0: Rotation vector is given by parameters INIE1,
C               INIE2 and INIE3.  Parameters INIE1, INIE2, INIE3 thus
C               specify a single rotation vector for the whole run of
C               program CRT.
C             INIDIR.NE.0: The information on the rotation is taken from
C               the input  SRC file.
C               Each source point in input file SRC must be
C               supplemented by directional vector E1DIR,E2DIR,E3DIR,
C               expressed with respect to the local Cartesian basis
C               vectors (before rotation).  The rotation vector is not
C               determined by input parameters INIE1, INIE2 and INIE3,
C               but is derived, for each source point, from the
C               directional vector:
C               the third local Cartesian basis vector is rotated
C               towards the directional vector along the rotation vector
C               perpendicular to the local Cartesian X3 axis and to the
C               specified directional vector E1DIR,E2DIR,E3DIR.
C               This option corrresponds to the minimum rotation angle
C               of the local Cartesian X3 axis towards the specified
C               directional vector.
C             Default: INIDIR=0
C     INIPAR=integer... Determines the parametrization of rays:
C             For INIDIM.LE.0:
C               INIPAR.LT.0: Different default values of components
C                 INIE1, INIE2 and INIE3 of the above described rotation
C                 vector.  This option is left just for backward
C                 compatibility.
C               INIPAR.EQ.1: Ray parameters are polar-like spherical
C                 coordinates (colatitude,longitude) connected with the
C                 local Cartesian coordinate system which basis vectors
C                 are given by the square root of the metric tensor at
C                 the initial point, optionally rotated about rotation
C                 vector given by parameters INIE1, INIE2 and INIE3.
C                 Equator plane coincides with the local (X1,X2)-plane.
C                 Zero longitude is determined by the positive local X1
C                 half-axis.  Longitude PI/2 then corresponds to the
C                 positive local X2 half-axis.  The zero colatitude
C                 corresponds to the positive local X3 half-axis.
C                 If SIN(colatitude).LT.0, the ray is reported out of
C                 the ray-parameter domain: IEND=71 in subroutine INIT2.
C               INIPAR.EQ.2: Ray parameters are geographic-like
C                 spherical coordinates (longitude,latitude) connected
C                 with the local Cartesian coordinate system which basis
C                 vectors are given by the square root of the metric
C                 tensor at the initial point.
C                 Equator plane coincides with the local (X1,X2)-plane.
C                 Zero longitude is determined by the positive local X1
C                 half-axis.  Longitude PI/2 then corresponds to the
C                 positive local X2 half-axis.  The latitude is positive
C                 in the direction given by the positive X3 half-axis.
C                 If COS(latitude).LT.0, the ray is reported out of
C                 the ray-parameter domain: IEND=71 in subroutine INIT2.
C               INIPAR.EQ.3: Azimuthal equidistant projection of a unit
C                 sphere is parametrized by 2 Cartesian coordinates
C                 centred at the projection point.  This option is
C                 suitable especially for reflection seismic studies.
C                 The unit vector tangent to the ray,  expressed in the
C                 local Cartesian coordinate system which basis vectors
C                 are given by the square root of the metric tensor at
C                 the initial point, is given by
C                   T1=PAR1*SIN(R)/R
C                   T2=PAR2*SIN(R)/R
C                   T3=     COS(R)
C                 with  R=SQRT(PAR1*PAR1+PAR2*PAR2).
C                 If R.GT.PI, the ray is reported out of the
C                 ray-parameter domain: IEND=71 in subroutine INIT2.
C             For INIDIM=1:
C               INIPAR must be 1 or 2.  The INIPAR-th ray parameter is
C               identical with the parameter parametrizing the initial
C               curve.  The other ray parameter is the angle between the
C               ray take-off plane and the normal to the interpolated
C               surface.  The ray take-off plane is given by the tangent
C               to the initial line and by the slowness vector.
C               For INIPAR=1, the initial line is the line PAR2=0 at the
C                 interpolated surface and is parametrized by PAR1.
C               For INIPAR=2, the initial line is the line PAR1=0 at the
C                 interpolated surface and is parametrized by PAR2.
C             For INIDIM=2:
C               Ray parameters are identical with two parameters
C               parametrizing the initial surface.
C               INIPAR.LE.0: Initial surface is described in terms of
C                 functions specifying the dependence of general
C                 coordinates (X1,X2,X3) on two parameters of the
C                 initial surface.
C               INIPAR.EQ.1: Initial surface is specified in the
C                 polar-like spherical coordinates (colatitude,
C                 longitude, radius) connected with the local Cartesian
C                 coordinate system which basis vectors are given by the
C                 square root of the metric tensor at the given point.
C                 Colatitude and longitude are the parameters, and the
C                 initial surface is determined by a function specifying
C                 the dependence of the radius on these parameters
C                 (colatitude and longitude).
C               INIPAR.GE.2: Initial surface is specified in the
C                 geographic-like spherical coordinates (longitude,
C                 latitude, radius) connected with the local Cartesian
C                 coordinate system which basis vectors are given by the
C                 square root of the metric tensor at the given point.
C                 Longitude and latitude are the parameters, and the
C                 initial surface is determined by a function specifying
C                 the dependence of the radius on these parameters
C                 (longitude and latitude).
C             Default: INIPAR=2
C     ADVANC=real... Initial point of the ray is shifted by distance
C             ADVANC perpendicularly to the initial surface or line,
C             or tangentially to the ray for the single initial point.
C             All initial and other quantities (except for the metric
C             tensor) are then evaluated at the shifted initial point.
C             Finally, the initial travel time is linearly updated
C             using the initial slowness vector.  This option may be
C             useful if the source is situated close to a structural
C             interface.
C             Default: ADVANC=0. (mostly sufficient)
C                                                    
C Parameters to control
C screen output:
C             The parameters are read by subroutine file
C             scro.for.
C                                                   
C   Parameters to control textual screen output:
C     SCRANSI=integer... Identification how the screen handles the ANSI
C             escape sequences.
C             SCRANSI=0: The screen does not support the ANSI escape
C               sequences.  Program 'crt.for' thus does not use the ANSI
C               escape sequences for screen output.
C             SCRANSI=1: The screen supports the ANSI escape sequences.
C               Program 'crt.for' thus uses the ANSI escape sequences
C               for screen output.
C             Default: SCRANSI=1 (good for MS-DOS with ansi.sys, and
C               for most contemporary versions of Linux).
C     SCRPLUS=integer... Identification how the screen handles the first
C             character of a record.
C             SCRPLUS=0: All characters of a record are displayed,
C               including the first character (e.g., Linux Fortran
C               compiler g77).  Program 'crt.for' thus does not use the
C               first character to control the screen output.
C             SCRPLUS.GT.0: The first character of a record is not
C               displayed. Since the ANSI escape sequences are preferred
C               to the first-column control characters, all options
C               SCRPLUS.GT.0 are equal if SCRANSI=1.
C             SCRPLUS=1: The first character of a record has no impact
C               upon the screen display.
C               (Note that program 'crt.for' uses '+' as the
C               first-column control character even in this case.)
C             SCRPLUS=2: If the first character of a record is '+',
C               the record overwrites the preceding record,
C               but '1' has no impact upon the screen display
C               (e.g., MS-DOS Fortran compilers by Lahey).
C               Program 'crt.for' thus uses '+' as the first-column
C               control character.
C             SCRPLUS=3: If the first character of a record is '+',
C               the record overwrites the preceding record.
C               If the first character of a record is '1',
C               the record is written on the first line of the screen.
C               Program 'crt.for' thus uses '+' and '1' as the
C               first-column control characters.
C             Default: SCRPLUS=1 (good for all compilers).
C     SCRWIDTH=integer... Width in characters of the textual part of the
C             screen.  SCRWIDTH cannot be smaller than 20.  Program
C             'crt.for' generates textual output of maximum line length
C             equal to SCRWIDTH.
C             Hint: Specify SCRWIDTH=20 if you use graphical screen
C                   output.
C             Default: SCRWIDTH=79
C     SCRHEIGHT=integer... Height in lines of the textual part of the
C             screen.  SCRHEIGHT influences only the textual output on
C             histories and ray endpoints, the number of lines used for
C             other textual output may exceed SCRHEIGHT.
C             SCRHEIGHT may be lower than the actual number of lines but
C             should not exceed that.  SCRHEIGHT=25 is a good choice for
C             IBM-PC's, but SCRHEIGHT=24 fits also VAX computers.
C             Used only if SCRANSI=1.
C             Default: SCRHEIGHT=25
C     CRTSCRO='string'... Verbose option of the CRT program, enabling to
C             select various levels of screen output.
C             The 'string' is composed of zero to several characters.
C             Each character may be entered either in lowercase or
C             uppercase.  The characters may be entered in any order,
C             the order of characters does not influence the output.
C             Each character enables a particular output:
C             'S' (source): To display the index of the source.
C             'W' (wave): To display the index of the elementary wave.
C             'R' (ray): To display the index of the ray being traced.
C               This option is not recommended if simultaneously
C               SCRANSI=0 and SCRPLUS.LE.1.
C             'A' (activity): To inform about the current activity of
C               the program, e.g., 'Aiming' or 'Tracing'.
C               This option is not recommended if simultaneously
C               SCRANSI=0 and SCRPLUS.LE.1.
C             'P' (parameters): To display two ray take-off parameters.
C               This option is not recommended if simultaneously
C               SCRANSI=0 and SCRPLUS.LE.1.
C             'H' (history): To display the indices of simple blocks,
C               complex blocks and interfaces passed through by a ray.
C               This option considerably slows down the calculation and
C               should be used just for debugging.
C               This option should be avoided if SCRANSI=0.
C             'E' (end of ray): To display reason IEND of the
C               termination of the computation of a ray.  See the
C               description of IEND in
C               RAY2
C               and INIT2.
C               This option should be avoided if SCRANSI=0.
C             'G' (graphics): Option to switch on and off the graphical
C               screen 
C               output.  This option affects the output only if the
C               CalComp
C               graphics in program
C               crt.for
C               is enabled.
C             A blank string means no screen output on the progress of
C               ray tracing.  Only the messages related to the start and
C               termination of the program are displayed.
C             Notes:
C               The default setting is relatively conservative and
C                 innocuous.
C               If the ANSI escape sequences are not supported, switch
C                 to SCRANSI=0.
C               If the ANSI escape sequences are supported, options
C                 CRTSCRO='SWR' and CRTSCRO='SWRA' may also be
C                 reasonably fast.
C                 These options may also be reasonably fast if the
C                 Fortran compiler supports option SCRPLUS=2.
C               If SCRANSI=0 but SCRPLUS=2, it is recommended to
C                 accumulate the screen output within a single line.
C                 In this case, only options 'S', 'W', 'R', 'A', 'P'
C                 are recommended.  The minimum numbers of output
C                 characters corresponding to individual options are:
C                 'S' 11, 'W' 9, 'R' 12, 'A' 7, 'P' 35, plus the number
C                 of used options (with 'RA' counted as a single option)
C                 minus 1 for spaces separating the output strings.
C                 The resulting number of characters should not exceed
C                 the value of SCRWIDTH.
C                 If SCRWIDTH is greater, up 8 spaces are used to
C                 separate the output strings.
C                 Without counting, the default value of SCRWIDTH=79
C                 is sufficient to accomodate any combination of options
C                 'S', 'W', 'R', 'A', 'P' within a single line.
C                 Options 'H' or 'E' always generate multiline screen
C                 output, which is awfully slow if SCRANSI=0.
C               SCRANSI=0, SCRPLUS=1 or 2, CRTSCRO=' ' correspond to
C                 optional screen output subroutine file 'scronul.for'
C                 in CRT version 5.60 or older.
C               SCRANSI=0, SCRPLUS=2, SCRWIDTH=79, CRTSCRO='SWRA'
C                 correspond to basic screen output subroutine file
C                 'scronum.for' in CRT version 5.60 or older.
C               SCRANSI=1, SCRPLUS=1 or more, SCRWIDTH=20,
C                 CRTSCRO='WRAPHEG' correspond to optional screen output
C                 subroutine file 'scropc.for' in CRT version 5.60 or
C                 older.
C             Default: CRTSCRO='SW'
C                                                   
C   Parameters to control graphical screen output:
C     The graphical
C     screen output
C     is available only if the
C     CalComp
C     graphics in program
C     crt.for is enabled,
C     and character 'G' in the input parameter CRTSCRO is specified.
C     SCRBBOX1=real, SCRBBOX2=real, SCRBBOX3=real, SCRBBOX4=real...
C             Bounding box (plotting area) of the graphical
C             screen
C             output.  Used only if the
C             CalComp
C             graphics in program
C             crt.for
C             is enabled.
C             SCRBBOX1, SCRBBOX2: CalComp coordinates of the left bottom
C               corner of the plotting area in the CalComp units.
C             SCRBBOX3, SCRBBOX4: CalComp coordinates of the right top
C               corner of the plotting area in the CalComp units.
C             Default: SCRBBOX1=2.77, SCRBBOX2=0.02, SCRBBOX3=10.98,
C               SCRBBOX4=8.48 (The rightmost 3/4 of the Landscape US
C               Letter form 11.0in * 8.5in mapped onto the screen area.
C               The leftmost 1/4 of the screen area is reserved for the
C               text output, which is controlled by SCRWIDTH.
C               Note that the analogous values for the metric A4 form
C               29.7cm * 21.0cm are SCRBBOX1=7.45, SCRBBOX2=0.00,
C               SCRBBOX3=29.66, SCRBBOX4=21.00)
C     SCRLINE=real... Estimated thickness of a line plotted on the
C             screen (in the CalComp units).  Program 'crt.for' uses
C             this value to set the thickness of spacing between lines.
C             Used only if the CalComp graphics is enabled, see program
C             crt.for.
C             Default: SCRLINE=0.017 (for low resolution and US Letter,
C               note that the analogous value for the metric A4 is
C               SCRLINE=0.045)
C     CRTPAUSE=integer... Positive value enables waiting for ENTER from
C             the standard input after each elementary wave.
C             This option is useful just when running program CRT
C             manually, with optional 'scro.for' and with CalComp
C             graphics erasing the screen when opening a new plot.
C             The program then waits to confirm erasing of the ray
C             diagram of each elementary wave.
C             Attention:
C               Do not enable this feature when executing the CRT
C               program from a history file.
C             Default: CRTPAUSE=0 (no pause, mostly sufficient)
C Note (very detailed):
C     Filenames MODEL, DCRT, INIT, CODE, RPAR, WRIT and CRTOUT need not
C     be mutually different, several data sets may be read from (or
C     written to) the same data file.  Each data file is closed when
C     read over, and its logical unit number may be connected to another
C     file to be opened.  When more than one elementary wave is
C     computed, it is not recommended to write the LOG output data set
C     to the file containing the CODE, RPAR or WRIT data set.
C 
C Example of data CRT all data sets separated
C 
C Example of data CRT with DCRT and INIT
C
C=======================================================================
C
C     
C
      SUBROUTINE CRTIN(FILE1,LUCODE,LURPAR,LUWRIT,LULOG)
      CHARACTER*(*) FILE1
      INTEGER LUCODE,LURPAR,LUWRIT,LULOG
C
C Subroutine CRTIN is designed to open the data files for complete ray
C tracing and to read the input data sets CRT, MODEL, DCRT and INIT.
C
C Input:
C     FILE1...The name of the main input data file CRT.  The file is
C             opened with the logical unit number LU(1)=5 defined in
C             this subroutine.  The name may be blank to use
C             preconnected input device.  Note that also logical units
C             LU(2)=4, LU(3)=3 and LU(4)=2 may be used to connect other
C             input data files always having non-blank filenames.
C
C Output:
C     LUCODE..The logical unit connected to the file with the CODE data.
C     LURPAR..The logical unit connected to the file with the RPAR data.
C     LUWRIT..The logical unit connected to the file with the WRIT data.
C     LULOG...The logical unit connected to the output LOG file.
C
C Subroutines and external functions required:
      EXTERNAL MODEL1,RAY1,INIT1,UNIT
C     MODEL1..File 'model.for' of the package 'model'.
C     RAY1... File 'ray.for'.
C     INIT1...File 'init.for'.
C     UNIT... This file.
C Note that the above subroutines reference many other external
C procedures from various subroutine files.  These indirectly
C referenced procedures are not named here, but are listed in the
C particular subroutine files.
C
C Date: 2014, May 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Quantities describing data files and logical units:
      INTEGER NFILE,IU,NU
      PARAMETER (NFILE=8)
      CHARACTER*80 FILE(NFILE)
      PARAMETER (NU=4)
      INTEGER LU(NU),KU(NU)
      DATA LU/5,4,3,2/
C     NFILE,FILE,IU,NU,LU,KU... For the description of these quantities
C             see subroutine UNIT below.
C
C.......................................................................
C
C     The name of the main input file.  This file contains the names of
C     other input files
      IF(NU.LT.4) THEN
C       103
        CALL ERROR('103 in CRTIN: Less than 4 logical units')
C       Four logical units must be available to read the input data and
C       to write the output log file.
      END IF
C
C     (1) Main data file - contains the names of other input files
      CALL RSEP1(LU(1),FILE1)
      CALL RSEP3T('MODEL' ,FILE(1),'model.dat' )
      CALL RSEP3T('DCRT'  ,FILE(2),'crt.dat'   )
      CALL RSEP3T('INIT'  ,FILE(3),' '         )
      CALL RSEP3T('CODE'  ,FILE(4),'code.dat'  )
      CALL RSEP3T('RPAR'  ,FILE(5),'rpar.dat'  )
      CALL RSEP3T('WRIT'  ,FILE(6),'writ.dat'  )
      CALL RSEP3T('CRTLOG',FILE(7),'crtlog.out')
      IF(FILE(1).EQ.' ') THEN
C       105
        CALL ERROR('105 in CRTIN: Name of input file MODEL is blank')
      END IF
      IF(FILE(2).EQ.' ') THEN
C       106
        CALL ERROR('106 in CRTIN: Name of input file CRT is blank')
      END IF
      IF(FILE(4).EQ.' ') THEN
C       107
        CALL ERROR('107 in CRTIN: Name of input data DCRT is blank')
      END IF
      IF(FILE(5).EQ.' ') THEN
C       108
        CALL ERROR('108 in CRTIN: Name of input data CODE is blank')
      END IF
      IF(FILE(6).EQ.' ') THEN
C       109
        CALL ERROR('109 in CRTIN: Name of input data RPAR is blank')
      END IF
      IF(FILE(7).EQ.' ') THEN
C       110
        CALL ERROR('110 in CRTIN: Name of output file CRTLOG is blank')
      END IF
C
C     (2) Data for model
      CALL UNIT(1,NFILE,FILE,IU,NU,LU,KU,'OLD')
      CALL MODEL1(LU(IU))
C
C     (3) Data for complete ray tracing
      CALL UNIT(2,NFILE,FILE,IU,NU,LU,KU,'OLD')
      CALL RAY1(LU(IU))
C
C     (4) Data for initial points of rays
      CALL UNIT(3,NFILE,FILE,IU,NU,LU,KU,'OLD')
      IF(IU.EQ.0) THEN
        CALL INIT1(0)
      ELSE
        CALL INIT1(LU(IU))
      END IF
C
C     (5) File containing the codes of elementary waves
      CALL UNIT(4,NFILE,FILE,IU,NU,LU,KU,'OLD')
C     The data file for the subroutine CODE1 remains open
      LUCODE=LU(IU)
      IU=0
C
C     (6) File controlling the take-off parameters of rays
      CALL UNIT(5,NFILE,FILE,IU,NU,LU,KU,'OLD')
C     The data file for the subroutine RPAR1 remains open
      LURPAR=LU(IU)
      IU=0
C
C     (7) File specifying the names of files with the computed
C     quantities
      CALL UNIT(6,NFILE,FILE,IU,NU,LU,KU,'OLD')
C     The data file for the subroutine WRIT1 remains open
      LUWRIT=LU(IU)
      IU=0
C
C     (8) The output LOG file
      CALL UNIT(7,NFILE,FILE,IU,NU,LU,KU,'UNKNOWN')
C     The output file for the subroutines WRIT1, WRIT2, WRIT4 and WRIT5
C     remains open
      LULOG=LU(IU)
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE UNIT(IFILE,NFILE,FILE,IU,NU,LU,KU,STATUS)
      INTEGER IFILE,NFILE,IU,NU,LU(NU),KU(NU)
      CHARACTER*(*) FILE(NFILE),STATUS
C
C Subroutine UNIT is designed to control connecting and disconnecting
C the data files to logical units, and to determine the logical units
C from which the given data sets are to be read.
C
C Input:
C     IFILE...Index of the data set to be read in.  The data sets are
C             indexed by integers from 1 to NFILE.
C     NFILE...The total number of data sets.
C     FILE... Character array containing the names of the files
C             containing individual data sets.  Different data sets may
C             be stored in the same file.  If IFILE=1, only FILE(1) has
C             to be defined.
C     IU...   Index of the logical unit connected to the file containing
C             the last read data set (i.e. the last data set was read
C             from the logical unit LU(IU) connected to the file
C             FILE(KU(IU)).  Zero if no data are read in, or if there is
C             no data file to close.  Need not be defined if IFILE=1.
C     NU...   The maximum number of available logical units.
C     LU...   Array containing reference numbers of logical units.
C     KU...   Auxiliary array of the dimension at least NU.
C             Its elements KU(1) to KU(NU) must not be modified between
C             two invocations of this subroutine.  Its values need not
C             be defined if IFILE=1.
C             KU(I)...Zero if the logical unit LU(I) is closed,
C               otherwise the sequential number of the next data set to
C               be read from this unit.
C
C Output:
C     IU...   Index of the logical unit connected to file with the data
C             set to be read in (i.e. the next data set will be read
C             from the logical unit LU(IU) connected to the file
C             FILE(IFILE)).
C             Zero if the filename is blank or if no logical unit is
C             available.
C     KU...   Updated input values.
C
C Date: 1999, May 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IERR,JU,I
C
      IF(IFILE.EQ.1) THEN
C       No logical units are connected when opening the first data file
        DO 10 JU=1,NU
          KU(JU)=0
   10   CONTINUE
      ELSE
C       Updating indices of next data sets to be read from open files:
        DO 13 JU=1,NU
          IF(0.LT.KU(JU).AND.KU(JU).LT.IFILE) THEN
C           The data set from the file connected to the logical unit
C           LU(JU) has been read.  Determining the next data set
C           contained in the file:
            DO 11 I=IFILE,NFILE
              IF(FILE(KU(JU)).EQ.FILE(I)) THEN
C               The I-th data set will also be read from the last file
C               connected to the logical unit LU(JU)
                KU(JU)=I
                GO TO 12
              END IF
   11       CONTINUE
C           No other data set will be read from the last file.  The file
C           may be closed and the logical unit LU(IU) disconnected
   12       CONTINUE
          END IF
   13   CONTINUE
C       Closing the data file:
        IF(IU.GT.0) THEN
C         There is a file submitted to be closed
          IF(0.LT.KU(IU).AND.KU(IU).LT.IFILE) THEN
C           No other data set will be read from this file.  The file
C           may be closed and the logical unit LU(IU) disconnected
            CLOSE(LU(IU))
            KU(IU)=0
          END IF
        END IF
      END IF
C
C     Opening new data file:
      IF(IFILE.GT.0) THEN
        DO 21 JU=1,NU
          IF(KU(JU).EQ.IFILE) THEN
C           The data file is already open
            IU=JU
            RETURN
          END IF
   21   CONTINUE
C       The data file has to be opened
        DO 22 JU=1,NU
          IF(KU(JU).EQ.0) THEN
C           The logical unit LU(JU) may be connected to the data file
            IF(FILE(IFILE).EQ.' ') THEN
              IU=0
            ELSE
              IU=JU
              KU(IU)=IFILE
              OPEN(LU(IU),FILE=FILE(IFILE)
     *                           ,STATUS=STATUS,IOSTAT=IERR,ERR=90)
            END IF
            RETURN
          END IF
   22   CONTINUE
C       No logical unit available
        IU=0
      END IF
      RETURN
C
   90 CONTINUE
C       104
        WRITE(*,'('' STATUS'',I5,'': '',A)') IERR,FILE(IFILE)
        CALL ERROR('104 in UNIT: Open file error')
C       Error encountered during execution of the OPEN statement.
      END
C
C=======================================================================
C