C Program 'CRT' for complete ray tracing
C
C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of:
C     CRT...  Main program controling the complete ray tracing.
C
C Input data read from the * external unit:
C     The data consist of a single character string, read by list
C     directed (free format) input.  Thus the string has to be enclosed
C     in apostrophes.  The interactive * external unit may be redirected
C     to the file containing the string.
C (1) 'CRT'
C     the string containing the name of the main input data file for the
C     complete ray tracing program.  The data file 'CRT' will be read in
C     by the subroutine CRTIN, and is described within the FORTRAN77
C     source code file 'crtin.for'.  Only the first 80 characters of the
C     string are significant.
C Default: 'CRT'='crt.dat'.
C
C=======================================================================
C
      PROGRAM CRT
C
C Main program for complete ray tracing.  This program reads the input
C data and then controls the complete ray tracing of the specified
C elementary waves.
C
C Input data: Main input data set containing the names of other input
C files and the name of the output log file is read in by the subroutine
C CRTIN of the file 'crtin.for'.  Thus, the structure of the main input
C data set is described in the file 'crtin.for'.
C The name of the main input data file is given by the first actual
C argument of the subroutine crtin called in the first executive
C statement of the main program.  It is blank in the original version.
C
C Subroutines referenced:
      EXTERNAL CRTIN,RAY2,INIT2,CODE1,RPAR1,RPAR2,RPAR4
      EXTERNAL WRIT1,WRIT2,WRIT4,WRIT5
C     MODEL1...File  'model.for' of the package 'model'.
C     CRTIN,UNIT... File 'crtin.for'.
C     RAY1,RAY2... File 'ray.for'.
C     INIT1,INIT2... File 'init.for'.
C     CODE1... File 'code.for'.
C     RPAR1,RPAR2,RPAR4... File 'rpar.for'.
C     WRIT1,WRIT2,WRIT4, WRIT5... File 'writ.for'.
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: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Storage locations:
C
C     Input data:
      CHARACTER*80 FCRT
C     FCRT... Name of the main input data file for the complete ray
C             tracing program.
C
C     Logical unit numbers:
      INTEGER LUCODE,LURPAR,LUWRIT,LULOG
C
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     Quantities describing elementary waves and their rays:
      REAL PAR1,PAR2,YL(6),Y(35),YY(5)
      INTEGER IWAVE,IWAVE0,IKODE,IRAY,IY(12),IEND,ISHEET,IREC
C
C     IWAVE...Index of the computed elementary wave.
C     IWAVE0..Index of the already computed elementary wave having the
C             most numerous common elements with the current elementary
C             wave.
C     IKODE...The length of the common part of the codes of the IWAVE-th
C             and IWAVE0-th elementary waves.
C     IRAY... Index of the computed ray.
C     PAR1,PAR2... Ray take-off parameters.
C     YL...   Array containing local quantities at a point of a ray (see
C             C.R.T.5.5.4).  The quantities are listed in the subroutine
C             file 'ray.for'.
C     Y...    Array containing basic quantities computed along a ray
C             (see C.R.T.5.2.1).  The quantities are listed in the
C             subroutine file 'ray.for'.
C     YY...   Array containing real auxiliary quantities computed along
C             a ray (see C.R.T.5.2.2).  The quantities are listed in the
C             subroutine file 'ray.for'.
C     IY...   Array containing integer auxiliary quantities computed
C             along a ray (see C.R.T.5.2.2).  The quantities are listed
C             in the subroutine file 'ray.for'.
C     IEND... Reason of the termination of the computation of a ray (see
C             C.R.T.5.4).
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.......................................................................
C
C     Opening data files and reading the input data:
C
      FCRT='crt.dat'
      WRITE(*,'(A)') ' Enter the name of the main input data file:'
      READ(*,*) FCRT
      WRITE(*,'(A)') '+Reading input data.                        '
      CALL CRTIN(FCRT,LUCODE,LURPAR,LUWRIT,LULOG)
      WRITE(*,'(A)') '+Computing.                                 '
C
C.......................................................................
C
C     Complete ray tracing:
C
      CALL WRIT1(LUWRIT,LULOG,0,0,0)
      IWAVE=0
      IRAY=0
C
C     Loop over elementary waves
   30 CONTINUE
C       Computation of a single elementary wave:
C       Reading the input data for the elementary wave
        CALL CODE1(LUCODE,IWAVE,IWAVE0,IKODE)
        IF(IWAVE.EQ.0) THEN
C         All required elementary waves are computed
          GO TO 90
        END IF
        CALL RPAR1(LURPAR,IWAVE)
        CALL WRIT1(LUWRIT,LULOG,IWAVE,IWAVE0,IKODE)
C       Loop over rays
   40   CONTINUE
C         Complete tracing of a single ray:
C         Determination of the take-off parameters
          CALL RPAR2(IRAY,PAR1,PAR2)
          IF(IRAY.EQ.0) THEN
C           All required rays of the elementary wave are computed
            GO TO 80
          END IF
C         Initial conditions for the ray
          CALL INIT2(PAR1,PAR2,YL,Y,YY,IY,IEND,IWAVE0,IKODE)
          CALL WRIT2(LULOG,IRAY)
          IF(IEND.EQ.0) THEN
C           Computation of the ray
            CALL RAY2(YL,Y,YY,IY,IEND)
          END IF
C         The ray is computed
          CALL RPAR4(IRAY,PAR1,PAR2,YL,Y,YY,IY,IEND,ISHEET,IREC)
          CALL WRIT4(LULOG,IRAY,YL,Y,YY,IY,IEND,ISHEET,IREC)
        GO TO 40
   80   CONTINUE
C       The elementary wave is computed
        CALL WRIT5(LULOG,IWAVE)
      GO TO 30
   90 CONTINUE
C
C     End of computation
      CALL WRIT5(LULOG,0)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'model.for'
      INCLUDE 'metric.for'
      INCLUDE 'srfc.for'
      INCLUDE 'parm.for'
      INCLUDE 'val.for'
      INCLUDE 'fit.for'
      INCLUDE 'means.for'
      INCLUDE 'hpcg.for'
      INCLUDE 'crtin.for'
      INCLUDE 'code.for'
      INCLUDE 'ray.for'
      INCLUDE 'raycb.for'
      INCLUDE 'trans.for'
      INCLUDE 'coef.for'
      INCLUDE 'init.for'
      INCLUDE 'rpar.for'
      INCLUDE 'rp2d.for'
      INCLUDE 'rp3d.for'
      INCLUDE 'writ.for'
C And one of the following files 'scro*.for':
*     INCLUDE 'scronul.for'
      INCLUDE 'scronum.for'
*     INCLUDE 'scropc.for'
C Note that 'scropc.for' must be supplemented with CalComp graphics
C subroutines PLOTS, PLOT, and NEWPEN.
C
C=======================================================================
C
