C
C File 'crtin.for' for reading the input data for complete ray tracing C C Date: 2001, July 6 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 are read by subroutine C MODEL1. See the description of the C 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, repeating C for each source point all elementary waves specified in C data set CODE for each source C point. Data in data sets RPAR C and WRIT are not repeated, C they should be specified for all elementary waves of all C source points. C Parameter SRC has no meaning for line or surface initial C conditions for rays. C Description of 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. C Description of 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 C description of 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 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 C description of 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 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 are read by 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. See the C description of data set RPAR C in subroutine file 'rpar.for'. 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 are read by subroutine WRIT1. 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 C description of data set WRIT C in subroutine file 'writ.for'. 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 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 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 C Parameter to control screen output, read by optional subroutine file C scropc.for: 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 'scropc.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. Do not enable this C feature when executing the CRT 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: 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