C
C Program CRT for complete ray tracing
C
C Version: 6.20
C Date: 2008, February 8
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C             File 'SEP' will be read by subroutine CRTIN, and the
C             parameters used by program CRT are described within
C             source code file
C             'crtin.for'.
C     No default, 'SEP' must be specified and cannot be blank.
C
C.......................................................................
C
C This file consists of:
C     CRT...  Main program controlling the complete ray tracing.
C             CRT
C
C=======================================================================
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 ERROR,CRTIN,RAY2,INIT2,CODE1,RPAR1,RPAR2,RPAR4
      EXTERNAL LUWARN,WRIT1,WRIT2,WRIT4,WRIT5
      INTEGER  LUWARN
C     ERROR,LUWARN...File  'error.for' of package FORMS.
C     MODEL1...File  'model.for' of 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: 2008, February 8
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Storage locations:
C
C     Input data:
      CHARACTER*80 FCRT,FRPAR,FWRIT
C
C     FCRT... Name of the main input data file for the complete ray
C             tracing program.
C     FRPAR.. Name of optional file
C             RPAREW
C             containing data RPAR-(4) corresponding to individual
C             elementary waves.
C     FWRIT.. Name of optional file
C             WRITEW
C             containing data WRIT-(6) corresponding to individual
C             elementary waves.
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.
C             
C             Description of YL
C     Y...    Array containing basic quantities computed along a ray,
C             see
C             C.R.T.5.2.1.
C             Description of Y
C     YY...   Array containing real auxiliary quantities computed along
C             a ray, see
C             C.R.T.5.2.2.
C             
C             Description of YY
C     IY...   Array containing integer auxiliary quantities computed
C             along a ray, see
C             C.R.T.5.2.2.
C             
C             Description of IY
C     IEND... Reason of the termination of the computation of a ray, see
C             C.R.T.5.4.
C             
C             Description of IEND in RAY2
C             
C             Description of IEND in INIT2
C     ISHEET..Ray-history index.  The different ray histories are
C             consecutively indexed by positive integers 1,2,3,...
C             According to their appearance during ray tracing.
C             The ray histories are indexed independently within each
C             elementary wave.
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
C             RPAR4.
C
C.......................................................................
C
C     Opening data files and reading the input data:
C
      WRITE(*,'(A)') '+CRT: Enter input filename: '
      FCRT=' '
      READ(*,*) FCRT
      WRITE(*,'(A)') '+CRT: Reading input data.   '
      IF(FCRT.EQ.' ') THEN
C       100
        CALL ERROR('100 in CRT: No input file specified')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      END IF
      CALL CRTIN(FCRT,LUCODE,LURPAR,LUWRIT,LULOG)
      LULOG=LUWARN(LULOG)
      WRITE(*,'(A)') '+CRT: Computing.            '
C
C.......................................................................
C
C     Complete ray tracing:
C
      CALL RPAR1(LURPAR,0)
      CALL WRIT1(LUWRIT,LULOG,0,0,0)
      CALL RSEP3T('RPAREW',FRPAR,' ')
      IF (FRPAR.NE.' ') THEN
        CLOSE(LURPAR)
        OPEN(LURPAR,FILE=FRPAR,STATUS='OLD',ERR=91)
      END IF
      CALL RSEP3T('WRITEW',FWRIT,' ')
      IF (FWRIT.NE.' ') THEN
        CLOSE(LUWRIT)
        OPEN(LUWRIT,FILE=FWRIT,STATUS='OLD',ERR=92)
      END IF
      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
          CALL INIT5(IRAY)
          IF(IRAY.EQ.0) THEN
C           One more source point - repeating loop over elementary waves
            REWIND(LUCODE)
            IF (FRPAR.NE.' ') THEN
              REWIND(LURPAR)
            END IF
            IF (FWRIT.NE.' ') THEN
              REWIND(LUWRIT)
            END IF
            GO TO 30
          ELSE
C           End of ray tracing
            GO TO 90
          END IF
        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)
      WRITE(*,'(A)') '+CRT: Done.                 '
      STOP
C
C     Open file errors:
   91 CONTINUE
C       101
        CALL ERROR('101 in CRT: Open file RPAREW error')
C       Error encountered when opening the file specified by input SEP
C       parameter
C       RPAREW.
      STOP
   92 CONTINUE
C       102
        CALL ERROR('102 in CRT: Open file WRITEW error')
C       Error encountered when opening the file specified by input SEP
C       parameter
C       WRITEW.
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'eigen.for'
C     eigen.for
      INCLUDE 'model.for'
C     model.for
      INCLUDE 'hder.for'
C     hder.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'hpcg.for'
C     hpcg.for
      INCLUDE 'crtin.for'
C     crtin.for
      INCLUDE 'code.for'
C     code.for
      INCLUDE 'ray.for'
C     ray.for
      INCLUDE 'raycb.for'
C     raycb.for
      INCLUDE 'trans.for'
C     trans.for
      INCLUDE 'coef.for'
C     coef.for
      INCLUDE 'init.for'
C     init.for
      INCLUDE 'rpar.for'
C     rpar.for
      INCLUDE 'rp2d.for'
C     rp2d.for
      INCLUDE 'rp3d.for'
C     rp3d.for
      INCLUDE 'writ.for'
C     writ.for
      INCLUDE 'scro.for'
C     scro.for
C
C Screen graphics:
C     To enable screen graphics, link this program with CalComp graphics
C     subroutines PLOTS, PLOT and NEWPEN designed for your operating
C     system and compiler, and comment the following line:
      INCLUDE 'plotnul.for'
C     plotnul.for
C
C=======================================================================
C