C 
C File 'crtin.for' for reading the input data for complete ray tracing
C
C Date: 1999, June 12
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 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 will be read in by the
C             subroutine MODEL1.  See
C             data file MODEL.
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 a single initial point and the initial value of the
C             travel time.
C             No meaning for line or surface initial conditions for
C             rays.
C             Description of file SRC
C             Default: SRC='src.dat'
C     REC='string'... Name of the input file containing the coordinates
C             Name of the file with profile or receiver coordinates.
C             Description of file REC
C     Default: 'REC'='rec.dat'.
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 will be
C             read in by the subroutine RAY1.
C             Description of data set DCRT
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             The data will be read in by the subroutine INIT1.
C             If INIT=DCRT, data set INIT is appended to file DCRT
C             (before possible data sets CODE, INIT or WRIT).
C             Description of data set INIT
C             Default: INIT=' '
C     CODE='string'... String containing the name of the file with the
C             codes of elementary waves.  The data will be read in by
C             the subroutine CODE1.
C             If CODE=DCRT, data set CODE is appended to file DCRT.
C             It is recommended to append at most one of sets CODE,
C             RPAR, WRIT to DCRT.
C             Description of data set CODE
C             Default: CODE='code.dat'
C     RPAR='string'... String containing the name of the file with the
C             data specifying the take-off parameters of the required
C             rays.  The data will be read in by the subroutine RPAR1.
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.
C             Description of data set RPAR
C             Default: RPAR='rpar.dat'
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 will be read by the subroutine
C             WRIT1.
C             If RPAR=WRIT, 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.
C             Description of data set WRIT
C             Default: WRIT='writ.dat'
C Names of the output files:
C     The names of the output files with the computed quantities
C             (C.R.T.5.5) are specified by the file given by input
C             parameter WRIT described above.
C     CRTLOG='string'... The string containing the name of the output
C             log file.  The data will be 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                                                    
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     INIPAR=integer... Determines the parametrization of rays:
C             For INIDIM.LE.0:
C               INIPAR.LT.0: The same as for IABS(INIPAR), but with
C                 the unit vector (T1,T2,T3) tangent to the ray
C                 changed to (T1,-T2,-T3).
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.
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.2*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 Note (too 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 Example of data CRT all data sets separated
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: 1999, May 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
C
C     Auxiliary storage locations:
      INTEGER I
C     I...    Auxiliary loop variable
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
C             quantities see the 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       102
        CALL ERROR('102 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')
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       101
        WRITE(*,'('' STATUS'',I5,'': '',A)') IERR,FILE(IFILE)
        CALL ERROR('101 in UNIT: Open file error')
C       Error encountered during execution of the OPEN statement.
      END
C
C======================================================================
C