C
C Program NET for a network shortest path calculation of the first C arrival travel time and the corresponding rays, on a rectangular grid C of points. C C Version: 3.10 C Date: 1998, September 26 C C Authors: C Ludek Klimes C Department of Geophysics, Charles University Prague C Ke Karlovu 3 C 121 16 Praha 2, Czech Republic C E-mail: klimes@seis.karlov.mff.cuni.cz C Michal Kvasnicka C Department of Geophysics, Charles University Prague C Ke Karlovu 3 C 121 16 Praha 2, Czech Republic C E-mail: qasnicka@seis.karlov.mff.cuni.cz C C....................................................................... C C Program 'NET' is able to: C (A) Perform both 2-D and 3-D ray tracing in a rectangular grid C discretization of an arbitrarily complex model according to C the paper by Klimes and Kvasnicka (1994). C (B) Compute first arrivals of both direct and multiply reflected C waves. C (C) Compute first arrivals of both P-waves, S-waves and converted C waves. C (D) Employ forward stars of any kind (template forward stars are C stored in a disk file, and may be modified and optimized in C various ways). C (E) Trace rays both in the whole rectangular model volume and in a C Fresnel volume corresponding to a two-point ray, or in a volume C of any other shape. This may increase accuracy of two-point ray C tracing by one order, or to decrease computational complexity C by about two orders. C (F) Adjust the size of a forward star at a structural interface in C order to maximize the accuracy. C (G) Adjust the size of a forward star according to a local degree of C heterogeneity in order to maximize the accuracy. C (H) Estimate maximum errors of all computed arrival times. C (I) In a special 2-D mode referred hereinafter as TTT perform 2-D C grid travel-time tracing of the second order according to the C paper by Klimes (1996). C C References: C Klimes L. and Kvasnicka M. (1994): 3-D network ray tracing. C Geophys.J.int., 116, 726-738. C Klimes L. (1996): Grid travel-time tracing: second-order method C for the first arrivals in smooth media. C PAGEOPH, 148, 539-563. C C....................................................................... C C Notes concerning the method: C C Although the program deals with arbitrarily complex models, the C accuracy of the shortest path approximation of the travel time and C rays, of course, depends on the model discretization step and on the C model complexity. Fortunately, the relative travel-time error can C be estimated. C C Because of the local optimization of the sizes of forward stars, the C network becomes oriented (i.e. backward stars do not equal forward C stars) and the resulting shortest path travel time need not be C reciprocal. C C This program employs two-point (trapezoidal) travel time approximation C along a network edge. C C Hereinafter the way of building the rectangular grid and the network C is briefly described. C C Big bricks, small bricks, gridpoints: C The rectangular volume bounded by coordinate limits X1MIN,X1MAX, C X2MIN,X2MAX, and X3MIN,X3MAX is divided into N1*N2*N3 big bricks. C If the numbers L1,L2,L3 are specified in addition to N1,N2,N3, C each big brick is subdivided into L1*L2*L3 small bricks. C If the numbers L1,L2,L3 are not specified, each big brick contains C just one small brick, as large as big one. C Gridpoints are defined as centres of the small bricks. C C Gridpoint indexing (positions of network nodes): C Gridpoints are indexed by integers 1,2,...,N1*L1*N2*L2*N3*L3. C Indices of gridpoints describe the their locations within the C model volume. The gridpoint indexing is similar to 1-D indexing C of 3-D Fortran array. A gridpoint index define the location of C the network node situated at the gridpoint. C C Network node indexing (addresses of network nodes): C If L1,L2,L3 are not specified, gridpoint indices equal to the C corresponding gridpoint indices. C Description for experienced users: C Network nodes are indexed just within the computational volume C composed of big bricks. Nodes are indexed consecutively within C each big brick. The node indexing within a big brick is similar C to 1-D indexing of 3-D Fortran array. The number of indices C within each big brick is L1*L2*L3. If the number of big bricks C within the computational volume is L4, then the node indices take C the values 1,2,...,L1*L2*L3*L4. C C Index file: C Index file enables to specify the computational volume of a C complex shape. The network ray tracing is then performed only C in the specified computational volume instead of the whole C rectangular model volume. This enables, e.g., to perform the C network ray tracing just in the estimated Fresnel volume. C The option of the index file is designed for experienced users. C The index file to specifies mapping of gridpoints onto the network C nodes. If it exists, it contains N1*N2*N3 items (integer C indices), each corresponding to one big brick. Each item contains C the index of the first node of the big brick divided by L1*L2*L3. C C Network nodes: C The set of network nodes if composed of: C (1) Given source nodes. C (2) The subset of gridpoints: C A gridpoint belongs to the network if contained in the big C brick belonging to the network. The relation of big C bricks with respect to the network is defined by means of C the index file: C If the index file is not specified, all big bricks belong C to the network. C If the index file is specified, each its item (index) C corresponds to one big brick: C If the index is zero, the big brick does not belong to C the network. C If the index is not zero, the big brick belongs to the C network and the index denotes the set of nodes within C the brick. C Big bricks belonging to the network are assumed to be C indexed by positive integers. C (3) Given receiver nodes. C C Generation of network edges: C Template forward star: C The concept of a template forward star is based on the C idea of a set of nodes (gridpoints) periodic in the space: C Template forward star is of the same geometry for all C network nodes. The template forward star is assumed C symmetric with respect to axis reflections and C interchangings. The edges corresponding to the template C forward star are read from the file 'net.fs*'. C Template forward stars are considered to contain the edge C edge connecting the central node with itself, although C this edge is not explicitly referred in the file C 'net.fs*'. C Size of a template forward star roughly describes the C extent of the forward star in grid intervals. C Many different kinds of template forward stars may be C considered: C Full 3-D cubic forward star, size NFS: C Contains all edges connecting the central node with all C nodes within the cube of the linear size of 2*NFS+1 C nodes, centred at the central node. C Full 2-D square forward star, size NFS: C Contains all edges connecting the central node with all C nodes within the square of the linear size of 2*NFS+1 C nodes, centred at the central node. C Full 3-D spherical forward star, size NFS: C Contains all edges connecting the central node with all C nodes within the sphere of the radius of SQRT(NFS*NFS+2) C grid intervals, centred at the central node. C Full 3-D circular forward star, size NFS: C Contains all edges connecting the central node with all C nodes within the circle of the radius of SQRT(NFS*NFS+1) C grid intervals, centred at the central node. C Moser-Saito 3-D cubic forward star, size NFS: C Full 3-D cubic forward star, size NFS, without the edges C that would be parallel in a homogeneous medium. C Moser-Saito 2-D square forward star, size NFS: C Full 2-D square forward star, size NFS, without the C edges that would be parallel in a homogeneous medium. C Optimized 3-D spherical forward star, size NFS: C Full 3-D spherical forward star, size NFS, without the C edges that can be removed while keeping, C in a homogeneous medium, each angular cone of the C radius grater than 1/(SQRT(2)*NFS) radians covered by C edges. C Optimized 2-D circular forward star, size NFS: C Full 2-D circular forward star, size NFS, without the C edges that can be removed while keeping, C in a homogeneous medium, each angle of the size greater C than 1/NFS radians covered by edges. C Note that this forward star is not a subset of the C optimized 3-D spherical forward star. C C Actual forward stars (without receiver nodes): C All forward stars corresponding to nodes situated within C a same small brick are the same, equal to the forward C star corresponding to its central gridpoint. C A forward star corresponding to a gridpoint is the C intersection of the template forward star with the set of C network nodes. C At source nodes, full spherical or circular template C forward stars are considered. C At gridpoint nodes, optimized spherical or circular C template forward stars from the files 'net.fs3' and C 'net.fs2' are considered. C Receiver backward stars: C The edges generated by the backward stars of the receiver C nodes are added to the edges generated by the above C forward stars. The backward stars are the same as forward C stars from the central nodes of the corresponding small C bricks. C At receivers, full spherical or circular template C backward stars are considered. C C....................................................................... C C Files 'net.fs2' and 'net.fs3' with template forward stars: C The files describing forward stars of various sizes must have the C names 'net.fs2' and 'net.fs3'. Forward stars submitted with this C version: C 'net.fs2'... Optimized spherical 2-D forward stars, sizes 1-100. C 'net.fs3'... Optimized circular 3-D forward stars, sizes 1-27. C C Note: C Files 'net.fs2' and 'net.fs3' may be reduced to approximately 2/3 C of their size if replacing all multiple spaces by a single space. C This may slightly save disk space and speed up program starting. C C TTT mode: C Template forward stars 'net.fs2' and 'net.fs3' are not required. C C....................................................................... C C C Data files: 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) 'SEP',/ C 'SEP'...String in apostrophes containing the name of the input SEP C parameter file containing the data specifying grid C dimensions, references to other data files and optionally C some numerical parameters. C Description of file SEP C Default: 'SEP'='net.h' C C C Data file SEP has the form of the SEP (Stanford Exploration Project) C parameter file: C All the data are specified in the form of PARAMETER=VALUE, e.g. C N1=50, with PARAMETER directly preceding = without intervening C spaces and with VALUE directly following = without intervening C spaces. The PARAMETER=VALUE couple must be delimited by a space C or comma from both sides. C The PARAMETER string is not case-sensitive. C PARAMETER= followed by a space resets the default parameter value. C All other text in the input files is ignored. The file thus may C contain unused data or comments without leading comment character. C Everything between comment character # and the end of the C respective line is ignored, too. C The PARAMETER=VALUE couples may be specified in any order. C The last appearance takes precedence. C Other data filenames: C NET=string... String with the name of the input data file C containing the references to other input and output files. C Description of input data NET C Default: NET='net.dat' C FSTAB=string... String with the name of the optional output file C containing the information on memory required to store C forward stars. Useful if array dimension MFS declared in C net.inc has to be updated. C The file is not created by default. C Default: FSTAB=' ' C Data specifying grid dimensions: C Boundaries of the model volume: C X1MIN=O1-0.5*D1, X1MAX=X1MIN+FLOAT(N1)*D1, C X2MIN=O2-0.5*D2, X2MAX=X2MIN+FLOAT(N2)*D2, C X3MIN=O3-0.5*D3, X3MAX=X3MIN+FLOAT(N3)*D3. C Because the model volume is divided into the small bricks C and the gridpoints are situated in the centres of the C small bricks, the left-hand model boundary is situated C half a grid interval to the left from the leftmost C gridpoints. Similarly for other model boundaries. C In a 2-D model, boundaries parallel with the model plane C are ignored. C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C N3 need not be specified for a 2-D model. C Default: N3=1 C D1=positive real... Grid interval in the direction of the first C coordinate axis. C Default: D1=1. C D2=positive real... Grid interval in the direction of the second C coordinate axis. C Default: D2=1. C D3=positive real... Grid interval in the direction of the third C coordinate axis. C Default: D3=1. C O1=real... First coordinate of the grid origin (first point of the C grid). C Default: O1=0. C O2=real... Second coordinate of the grid origin. C Default: O2=0. C O3=real... Third coordinate of the grid origin. C Default: O3=0. C Additional parameters for network ray tracing in Fresnel volumes: C L1,L2,L3 may be left out and should only be specified by C experienced users. C (If the numbers L1,L2,L3 are not specified, each big brick C contains just one small brick, as large as big one.) C L1=positive integer... Number of small bricks in one big brick in C the direction of axis X1. If specified, must be positive. C Default: L1=1 C L2=positive integer... Number of small bricks in one big brick in C the direction of axis X2. If specified, must be positive. C Default: L2=1 C L3=positive integer... Number of small bricks in one big brick in C the direction of axis X3. If specified, must be positive. C Default: L3=1 C Note: Only one of N1*L1, N2*L2, N3*L3 may equal 1. C In such a case, two-dimensional ray tracing is performed: C Two-dimensional forward stars are considered, C model boundaries in the proper direction are ignored, C source and receiver coordinates perpendicular to the plane C of computation have to equal. C Numerical parameters: C NFSMAX=integer... Options: NFSMAX=positive integer, 0, -1. C Default: NFSMAX=0. C NFSMAX=positive integer: C NFSMAX is the maximum size of a forward star. C Simultaneously, the default size of a forward star if the C size is not adjusted. C The size of a forward star is measured in grid intervals, C and describes the maximum length of the edge between two C consecutive points (nodes) of a ray. C The maximum size of a forward star has to be estimated by C the user according the model, in order to prevent the ray C to skip over some details in the model. These details may C be, e.g., small low or high velocity regions or bumps on C the structural interfaces. C NFSMAX has to be positive and not exceeding dimension C MFSMAX in the common block /FS/. C NFSMAX=0: The maximum size of forward stars is determined C automatically: New NFSMAX**(-2) is average of NFS**(-2), C where NFS is optimum size at each node. This is default. C TTT mode (specified by NFSMAX=-1): C N1,N2,N3... N3=1 is obligatory in the input data. C L1,L2,L3... L1=L2=L3=1 is strongly recommended. C NFSMAX=-1: 2-D 'Second-order' grid travel-time tracing by Klimes C (1996) is used instead of network ray tracing. C Default: NFSMAX=0. C RIDGE1=positive real, RIDGE2=positive real... C Numerical parameters controlling the check for C the ridges of the first-arrival travel time. At the C ridges, the rays going from different directions meet. C The interpolation across the ridges should be avoided C because it would introduce the first-order errors. C Instead, paraxial approximations from gridpoints at each C side of the ridge should be applied. For the paraxial C approximations, the second travel-time derivative along C the gridline crossing the ridge is taken zero. C The serious problem is to indicate the ridge. C The slowness-vector difference across the ridge is C negative and lower than would correspond to the C slowness-vector difference along the neighbouring gridline C segment, i.e. to the second travel-time derivative at the C corresponding side of the ridge. C The following algorithm is now applied: C The slowness-vector difference along the neighbouring C gridline is multiplied by RIDGE1 if negative or by 0 if C positive, and decreased by RIDGE2 times the norm of the C slowness gradient. If the result is greater than the C slowness-vector difference along the gridline being C checked, the ridge is indicated. C RIDGE1 should thus aways be greater than 1, whereas C RIDGE2 should be positive, small with respect to 1. C Values like RIDGE1=1.1, 1.2, 1.3, 1.4, 1.5, 3.0, 6.0, C and RIDGE2=0.05, 0.10, 0.20, 1.00 have been tested with C very unstable results. The questionable defaults of C RIDGE1=1.5 and RIDGE2=0.10 have finally been chosen. C Very large values of RIDGE1 and RIDGE2 should disable C the test for ridges and thus correspond to program TTT-2D C version 0.17 (slightly fixed version 0.07). C Default: RIDGE1=1.5, RIDGE2=0.1. C VER1=real... If the distance of the centre of wavefront curvature C from the gridline segment is greater than VER1 grid C intervals: C Cubic interpolation of travel time along the gridline C segment, C else: cubic interpolation of travel time residual along C the gridline segment. The travel time residual is taken C with respect to spherical wavefronts having the common C centre at the point of intersection of two given C slowness vectors. C For more details refer to Klimes (1996), Sec.3.1.4, where C VER1 is denoted by N. C The results are not very sensitive to this parameter and C there has never arisen a need to change the default value. C Default: VER1=9.999. C VER2=real... Smoothing and stabilizing slowness vectors: C Residual travel times with respect to quadratic or C spherical interpolation are approximated by C TTRES=(1-VER2)*TTCUB+VER2*TTLIN, C where TTCUB is cubic interpolation, and TTLIN is linear C approximation adjusting the resulting slowness vector C according to travel time residuals. C VER1=0 corresponds to interpolation described in Sections C 3.1.3.1 and 3.1.4.1 and is default in TTT-2D ver. 0.06. C VER1=1 corresponds to interpolation described in Sections C 3.1.3.2 and 3.1.4.2 and is default since TTT-2D version C 0.07. C Default: VER2=1.0. C Example of data set SEP C C C Main input file 'NET': C Sequential file, read by list directed (free format) input, C containing model parameters, source/receiver coordinates, and C names of other input and output files. C In the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new read statement). C 'ITEMS' in the list of input variables enclosed in apostrophes C represent character strings enclosed in apostrophes. Otherwise, C if the first letter of the symbolic name in the list of input C variables is I-N, the corresponding value in input data is C integer, otherwise, the input parameter is of the type real. C / in the list of input variables indicates an obligatory slash. C The slash may also be used instead of default values. C (1) 'SRC','REC','RAYS','END',/ C 'SRC'...Name of the input file with source coordinates. C Description of input data SRC C 'REC'...Name of the input file with receiver coordinates. C If blank, no rays are stored within the file 'RAYS'. C Description of input data REC C 'RAYS'..Name of the output file with rays. C If blank, no rays are stored within the file 'RAYS'. C Description of output data RAYS C 'END'...Name of the output file with endpoints of rays (receiver C coordinates, receiver arrival times, and estimates of the C corresponding maximum travel-time errors. C If blank, no file 'END' is generated. C Description of output data END C Default: 'REC'=' ', 'RAYS'=' ', 'END'=' '. C TTT mode: C 'REC'...File read if specified but its data ignored. No receivers C are considered in the present version. C 'RAYS'..Ignored - output file not generated. C 'END'...Ignored - output file not generated. C (2) NREFL,/ C NREFL...Number of reflections. C Attention concerning this version: C If, in the case of reflections (NREFL.GT.0), some C 'VEL(I)' differs from 'VEL(0)' at lines (3.2) below, the C errors cannot be evaluated correctly and the filename C 'END' above should be left blank. C Note: To calculate a reflected wave, it is recommended to C keep NREFL=0 and to submit the reflector points as a set C 'REC' of the receiver points. Then to submit the output C 'END' file as the source file 'SRC' for the subsequent C calculation of the reflected wave. In this way, the C arrival time errors connected with the discretization of C the reflector can be removed. On the other hand, the C whole rays are split into several files 'RAYS' and have C to be put together after finishing the network ray C tracing. C Default: NREFL=0. C (3) Once (3.1), then NREFL-times (3.2) and (3.1): C TTT mode: NREFL=0. C (3.1) 'IND(I)','VEL(I)','ICB(I)','TT(I)','ERR(I)','PRED(I)','NFS(I)',/ C 'IND(I)'... Name of the input index file, specifying for each C big brick if its gridpoints belong to the network. C May be blank - then the default indexing is assumed. C Must not be blank if (L1.GT.1.OR.L2.GT.1.OR.L3.GT.1) at C input (1). C This file should only be used by experienced users, others C should leave it blank. C Description of input data IND(I) C 'VEL(I)'... Name of the input file containing velocities at all C network nodes, for I-times reflected wave. C Has to be specified. C Description of input data VEL(I) C 'ICB(I)'... Name of the input file containing indices of C (geological) blocks. Only network edges corresponding to C a forward star of the size NFS=1 are allowed to cross an C interface between two different blocks (i.e. blocks with C different indices). This considerably increases the C accuracy in presence of structural interfaces (often 5 C times). Note that the limitation to size 1 is considered C to be the optimum one, at least, for velocity contrasts of C 1.18/1 or greater. C 'ICB(I)' may be blank - especially in the case of a single C block (no structural interfaces). C In the presence of structural interfaces, this file is C recommended to be submitted in order to increase the C accuracy. C Description of input data ICB(I) C 'TT(I)'... Name of the output file containing travel-times at all C network nodes after I reflections. C If blank, the file is not generated. C Description of output data TT(I) C 'ERR(I)'... Name of the output file containing estimated upper C bounds for the errors of the computed travel-times at all C network nodes after I reflections. C If blank, the file is not generated. C Attentions concerning this version: C (a) in the presence of structural interfaces, the errors C are evaluated correctly only if the file 'ICB' is C submitted. C (b) if, in the case of reflections, 'VEL(I)' differs from C 'VEL(I-1)', the errors cannot be evaluated correctly and C the corresponding 'ERR(I)' should be left blank. C Description of output data ERR(I) C 'PRED(I)'... Name of the output file containing predecessors of C all network nodes after I reflections. C If blank, the file is not generated. C May be blank for most applications. C Description of output data PRED(I) C 'NFS(I)'... Just for debugging. The parameter should be left out C or blank in the input data. C String 'NFS(I)' may be blank, '*' or the name of the file C containing optimum sizes of forward stars, preceded by + C for input or by - for output. C 'NFS(I)'=' ': optimum sizes of forward stars at the C network nodes are estimated and used for calculation. C 'NFS(I)'='+...' (input): sizes of forward stars at the C network nodes are read from the file. C 'NFS(I)'='-...' (output): optimum sizes of forward stars C at the network nodes are estimated, used for calculation C and written to the file. C 'NFS(I)'='*': sizes of all forward stars equal NFSMAX from C the input (2). C Description of data NFS(I) C Default: 'IND(I)'=' ', 'VEL(I)'=' ', 'ICB(I)'=' ', C 'TT(I)'=' ','ERR(I)'=' ', 'PRED(I)'=' ', 'NFS(I)'=' '. C TTT mode: C 'IND(I)'... Should be blank. C No index file may be used in the presend version. C 'ICB(I)'... Indices of geological blocks are read but C not used. Should be blank. C 'ERR(I)'='P1(I)'... Name of the output file containing the first C slowness vector component. C If blank, the file is not generated. C 'PRED(I)'='P2(I)'... Name of the output file containing the second C slowness vector component. C If blank, the file is not generated. C 'NFS(I)'='P3(I)'... Name of the output file reserved for the third C slowness vector component. Presently should be blank. C No file is generated. C (3.2) 'INTF(I)',/ C 'INTF(I)'... Name of the input file containing refractor points. C Must not be blank if the reflection is considered. C Description of input data INTF(I) C Example of data set NET C C C Input file 'SRC' containing source coordinates: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) Several times (2.1): C (2.1) 'NAMSRC',X1S,X2S,X3S,TTS,TTSERR/ C 'NAMSRC'... Name of the source point. Truncated to the first six C characters. C X1S,X2S,X3S... Coordinates of a point of the source. C In a 2-D model, coordinates perpendicular to the model C plane have to be the same for all source points and all C receivers. C TTS... Initial arrival time at the source point. Must not be C negative. Zero initial time is o.k., but is changed to C a very small positive value. C TTSERR..Initial value of the arrival time error at the source C point. It is likely zero at the actual source and may C thus be omitted. It is introduced especially for the C purposes of tracing reflected waves. In such a case, C TTS is the arrival time of incident wave at the reflector C and TTSERR is its error resulting from network ray tracing C between the actual source and the reflector. C Default: TTS=0, TTSERR=0. C TTT mode: C Only a point source is considered in the present version. C (3) / (a slash). C Example of data set SRC C C C Input file 'REC' containing receiver coordinates: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) Several times (2.1): C (2.1) 'NAMREC(IR)',X1R(IR),X2R(IR),X3R(IR),/ C 'NAMREC(IR)'... Name of the receiver point. Truncated to the C first six characters. C X1R(IR),X2R(IR),X3R(IR)... Coordinates of the IR-th receiver. C (3) / (a slash). C Example of data set REC C C C Output file 'RAYS': C (1) / (a slash). C (2) For each receiver (2.1), (2.2), and (2.3): C (2.1) 'RAYnnnnn FROM NAMSRC TO NAMREC',/ C 'RAYnnnnn FROM NAMSRC TO NAMREC'... String in apostrophes, C composed of: C (a1) substring 'RAY', C (a2) the sequential index of the receiver written using C format (i5), C (b1) substring ' FROM ', C (b2) name of the source point truncated to the first 6 C characters, C (c1) substring ' TO ', C (c2) name of the receiver point truncated to the first 6 C characters. C If the name of the source point is blank, (b1) is replaced C by 6 spaces. C If the name of the receiver point is blank, (c1) and (c2) C are left out. C If both the names of the source and receiver point are C blank, (b1), (b2), (c1) and (c2) are left out. C Examples of the string: C 'RAY 112' C 'RAY 112 TO NAMREC' C 'RAY 112 FROM NAMSRC' C 'RAY 112 FROM NAMSRC TO NAMREC' C (2.2) For each point of the ray (2.1.1): C (2.2.1) X1,X2,X3,TT,/ C X1,X2,X3... Coordinates of the point of the ray. C TT... Arrival time at the point. C Note: Rays are stored backwards, starting with the receiver, C ending with the source point. C (2.3) / (a slash). C (3) / (a slash). C C C Output file 'END': C (1) / (a slash). C (2) For each receiver (2.1): C (2.1) 'NAMREC',X1R,X2R,X3R,TT,TTERR,/ C 'NAMREC'... Name of the receiver point truncated to the first six C characters. C X1R,X2R,X3R... Coordinates of the receiver. C TT... Arrival time at the receiver point. C TTERR...Estimated maximum error of the arrival time. C Attentions concerning this version: C (a) in the presence of structural interfaces, the errors C are evaluated correctly only if the file 'ICB' is C submitted. C (b) if, in the case of reflections, 'VEL(I)' differs from C 'VEL(I-1)', the errors cannot be evaluated correctly and C the corresponding 'ERR(I)' should be left blank. C (3) / (a slash). C C C Input index files 'IND(I)': C (1) (IND(I),I=1,N1*N2*N3) C IND(I)..Zero if the I-th big brick does not belong to the C network, C otherwise the index the big brick. The nodes within the C big brick are indexed by integers from C L1*L2*L3*(IND(I)-1)+1 to L1*L2*L3*IND(I). C Default: IND(I)=I. C C C Input files 'VEL(I)' containing velocities at network nodes: C (1) (V(I),I=1,L1234) C V(I)... Velocity at network node no.I. C L1234.. Maximum index of a network node. L1234=L1*L2*L3*L4, where C L4 is the maximum index of a big brick belonging to the C network, see the file 'IND' below. If the file 'IND' C is not specified, L4=N1*N2*N3 by the default. C Free space is indicated by a zero velocity V(I)=0. C Free space should not be indicated by a small positive C velocity. A small velocity considerably increases time of C computation. C Default: V(I)=1. C TTT mode: C At structural interfaces, it is recommended to specify average C slownesses over bricks. However, the method is unstable in block C models. C Example of data set VEL C C C Input index files 'ICB(I)': C (1) (ICB(I),I=1,L1234) C ICB(I)..Index of (geological) block in which the network node no.I C is situated. C Default: ICB(I)=1. C C C Output files 'TT(I)' (time field): C (1) (TT(I),I=1,L1234) C TT(I)... Travel time at network node no.I. C L1234.. Maximum index of a network node. C The output file is designed to be read by the list-directed input C (free format). The null values are generated in place of undefined C travel times in a free space. The null values are treated as default C values when read by list-directed input (free format). C Example: 124 null values are written as ' 124*'. C C C Output files 'ERR(I)' (time error field): C (1) (ERR(I),I=1,L1234) C ERR(I)..Estimated maximum absolute value of the travel-time error C at network node no.I. C Please, notice that the absolute (i.e. not relative) C travel-time errors are estimated and written. C L1234.. Maximum index of a network node. C The output file is designed to be read by the list-directed input C (free format). The null values are generated in place of undefined C errors in a free space. C C C Output files 'PRED(I)' (predecessors): C (1) (IPRED(I),I=1,L1234) C IPRED(I)... Predecessor to network node no.I. C IPRED(I)=I: node has no predecessor (is in a free space). C IPRED(I)=0: node belongs to a reflector, preceding node C has the same position. C IPRED(I).LT.0: preceding node is (-IPRED(I))-th source C point. C Otherwise, IPRED(I).GT.0.AND.IPRED(I).NE.I.: preceding C node is a gridpoint, IPRED(I) being its position index. C L1234.. Maximum index of a network node. C C C Output files 'NFS(I)': C (1) (NFS(I),I=1,L1234) C NFS(I)..Optimum forward star size at the network node no.I. C If NFS.GT.NFSMAX, it is automatically decreased to NFSMAX, C where NFSMAX takes its positive input value or is C determined automatically if NFSMAX=0 on input. C C C Input files 'INTF(I)' containing indices of gridpoints forming the C reflectors: C (1) (INTF(I),I=1,N),/ C INTF(I)... Index of the I-th gridpoint forming the reflector. C N... Number of points forming the reflector. C C....................................................................... C C List of routines: C C High-level: C NET... Main program calling subroutines: C INPUT, SOURCE, GENER, RESFIL, ERRORS, RECVRS, TRACER, and C OUTPUT. C NET C INPUT...Reads the whole main input data file 'NET', source and C receiver coordinates. C INPUT C SOURCE..Reads input data for each iteration, and initializes the C arrival-time field and other arrays. Under iteration we C understand here the computation between two reflections. C If reflectors are not specified, there is just one C iteration. C SOURCE C GENER...Performs shortest path ray tracing (one iteration). C GENER C RESFIL..In a case of reflections, stores the results of an C iteration in an unformatted direct-access scratch file C (one iteration). C RESFIL C ERRORS..Generates the field of arrival-time errors by means of C backward ray tracing - writes 'ERR(I)' (one iteration). C ERRORS C RECVRS..Updates arrival times at the receivers. C RECVRS C TRACER..Performs backward ray tracing from the receivers, C including the estimation of arrival-time errors. C TRACER C OUTPUT..Rewrites arrival-time fields and predecessors from memory C or scratch files to formatted output files 'TT(I)' and C 'IPRED(I)'. C OUTPUT C Updating: C LOOP,FREVOL,UPDATE,SRCREC... Auxiliary routines to SOURCE, GENER, C and RECVRS, updating a node of a network. C Loop with FREVOL or UPDATE: employed by GENER. C Loop with SRCREC: employed by SOURCE and RECVRS. C LOOP C FREVOL C UPDATE C SRCREC C Template forward stars: C READFS..Called by SOURCE, to read and construct optimized template C forward stars. C READFS C MAKEFS..Called by SOURCE and RECVRS, to construct full spherical C template forward stars. C MAKEFS C SETFS...Auxiliary subroutine to READFS and MAKEFS. C SETFS C STORFS..Auxiliary subroutine to SETFS. C STORFS C Backward ray tracing: C ONERAY..Subroutine called by ERRORS and TRACER, to perform C backward ray tracing. C ONERAY C Accuracy estimation: C SETERR..Subroutine called by ONERAY, to estimate and accumulate C the arrival-time errors during backward ray tracing. C SETERR C OPTNFS..Subroutine called by SOURCE and SETERR, to set up C optimum sizes of forward stars. C OPTNFS C OPTMAT..Auxiliary subroutine to SETERR and OPTNFS, to calculate C the matrix composed of the first and second slowness C derivatives. C OPTMAT C TRYMAT..Auxiliary subroutine to OPTMAT, to try the calculation of C the matrix. C TRYMAT C MIXDER..Auxiliary subroutine to TRYMAT, to calculate mixed second C partial slowness derivatives. C MIXDER C Slowness interpolation: C SLOW... Subroutine called by SOURCE, RECVRS, TRACER, and ONERAY, C to find a close gridpoint, and to interpolate the slowness C within the same geological block. C SLOW C POS... Auxiliary subroutine to SLOW, to find the nearest C gridpoint. C POS C Indexing of nodes: C INDX... Integer function, called by many subroutines, to return C the index of the node at the given gridpoint. C INDX C C....................................................................... C C External subroutines required and not contained within this file: C C Formatted output: C WARRAY..Auxiliary subroutine to ERRORS and OUTPUT, designed to C write the given array into the file. C Source code file 'forms.for'. C Eigenvalues: C EIGEN...Subroutine from the IBM scientific subroutine package, C called by OPTNFS. C Source code file 'eigen.for'. C C----------------------------------------------------------------------- C C C PROGRAM NET C C----------------------------------------------------------------------- C C Common blocks: INCLUDE 'net.inc' C net.inc C All common blocks are declared here. C C----------------------------------------------------------------------- C INTEGER IREFL * INTEGER IT1,IT2 * REAL TIME C C IREFL.. Number of reflections. IREFL=0: direct wave. C IT1,IT2,TIME... Auxiliary variables for time watching. C C....................................................................... C C Reading the input parameters CALL INPUT() C C Loop over reflections: DO 10 IREFL=0,NREFL C C Generating time field of the IREFL-times reflected wave: CALL SOURCE(IREFL) C C Example timer for Lahey Fortran77 compiler: * IT1=0 * CALL TIMER(IT1) C C Calculation of shortest paths and arrival times: CALL GENER(IREFL) C C Example timer for Lahey Fortran77 compiler: * IT2=0 * CALL TIMER(IT2) * TIME=0.01*(FLOAT(IT2)-FLOAT(IT1))/60. * WRITE(*,'(A,F10.3,A)') ' Time=',TIME,' min.' C C Post-processing: CALL RESFIL(IREFL) CALL ERRORS(IREFL) C 10 CONTINUE C End of loop over reflections. C C End of ray computation and writing the results: CALL RECVRS() CALL TRACER() CALL OUTPUT() STOP END C C======================================================================= C C C SUBROUTINE INPUT() C C Procedure to read the whole main input data file 'net', source and C receiver coordinates. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ /TTTPAR/ are C required here: INCLUDE 'net.inc' C net.inc INCLUDE 'ttt.inc' C ttt.inc C C----------------------------------------------------------------------- C INTEGER LU1,LU2 REAL DWARF,UNDEF PARAMETER(LU1=1) PARAMETER (DWARF=1.E-10,UNDEF=-999999.) C CHARACTER*80 TEXT INTEGER ISRC,IREC,I,L REAL D1,D2,D3,O1,O2,O3 C REAL A1111,A1122,A2222,A1133,A2233,A3333,A1123 C REAL A2223,A3323,A2323,A1113,A2213,A3313,A2313 C REAL A1313,A1112,A2212,A3312,A2312,A1312,A1212 C C TEXT...Names of input files NET and SEP, then a dummy storage C location for a text. C ID(II),JD(II),KD(II)... Vector specifying II-th node of the C forward star, in dimensions of a small brick. C A1111,A1122,A2222,A1133,A2233,A3333,A1123,A2223,A3323,A2323,A1113, C A2213,A3313,A2313,A1313,A1112,A2212,A3312,A2312,A1312, C A1212... Intended for future extension (anisotropy). C C....................................................................... C C Name of the main input data file is read from the * external unit: TEXT='net.h' WRITE(*,'(A)') * '+NET: Enter name of main input data file [''net.h'']: ' READ(*,*) TEXT C End of reading from the * unit. C C Reading SEP parameter file: CALL RSEP1(LU1,TEXT) C C Reading input data file NET: CALL RSEP3T('NET',TEXT,'net.dat') OPEN(LU1,FILE=TEXT,STATUS='OLD') WRITE(*,'(2A)') '+NET: Reading input data file: ',TEXT(1:36) C C Numbers of gridpoints: CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3I('L1',L1,1) CALL RSEP3I('L2',L2,1) CALL RSEP3I('L3',L3,1) IF(N1.LE.0) THEN LN1=.TRUE. N1=1 L1=MAX0(1,L1) ELSE LN1=.FALSE. END IF IF(N1.LT.1.OR.N2.LT.1.OR.N3.LT.1.OR. * L1.LT.1.OR.L2.LT.1.OR.L3.LT.1) THEN C NET-10 CALL ERROR('NET-10: Number of gridpoints is not positive') C NET-10: Number of gridpoints is not positive: C N1,N2,N3, and, if specified, also L1,L2,L3, in input C file SEP must be positive. END IF NL1=N1*L1 NL2=N2*L2 NL3=N3*L3 C C Boundaries of the model volume: CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) IF(D1.LE.0..OR.D2.LE.0..OR.D3.LE.0.) THEN C NET-56 CALL ERROR('NET-56: Grid interval D1, D2 or D3 is not positive') C NET-56: Grid interval D1, D2 or D3 is not positive: C Grid interval D1, D2 and D3 must be positive in this C version of the NET program, see the C description of input file SEP. END IF X1MIN=O1-0.5*D1 X2MIN=O2-0.5*D2 X3MIN=O3-0.5*D3 X1MAX=X1MIN+FLOAT(N1)*D1 X2MAX=X2MIN+FLOAT(N2)*D2 X3MAX=X3MIN+FLOAT(N3)*D3 C Space steps in the rectangular grid IF(LN1) THEN DSX1=(X1MAX-X1MIN)/FLOAT(NL2) ELSE IF(NL1.EQ.1) THEN DSX1=0. ELSE DSX1=(X1MAX-X1MIN)/FLOAT(NL1) END IF IF(NL2.EQ.1) THEN DSX2=0. ELSE DSX2=(X2MAX-X2MIN)/FLOAT(NL2) END IF IF(NL3.EQ.1) THEN DSX3=0. ELSE DSX3=(X3MAX-X3MIN)/FLOAT(NL3) END IF IF(LN1) THEN ASX1=0. ASX2=SQRT(DSX1**2+DSX2**2) ELSE ASX1=ABS(DSX1) ASX2=ABS(DSX2) ENDIF ASX3=ABS(DSX3) C C Size of the forward star and the interpolation method for TTT: C Size of the forward star: CALL RSEP3I('NFSMAX',NFSMAX,0) IF(NFSMAX.GT.MFSMAX) THEN C NET-07 CALL ERROR('NET-07: Maximum size of a forward star too great') C NET-07: Maximum size of a forward star too great: C NFSMAX in input data file SEP should C be decreased if positive, adjusted to MFSMAX if zero, C or the dimension MFSMAX in the common block /FS/ C in net.inc should be increased. END IF C Interpolation method: CALL RSEP3R('RIDGE1',RIDGE1,1.5 ) CALL RSEP3R('RIDGE2',RIDGE2,0.10 ) CALL RSEP3R('VER1' ,VER1 ,9.999) CALL RSEP3R('VER2' ,VER2 ,1.000) C C (1) Names of the files with source, receivers, rays, and errors: FREC=' ' FRAYS=' ' FEND=' ' READ(LU1,*) FSRC,FREC,FRAYS,FEND IF(FREC.EQ.' '.AND.(FRAYS.NE.' '.OR.FEND.NE.' ')) THEN C NET-12 CALL ERROR * ('NET-12: Name of input file with receivers not specified') C NET-12: Name of input file with receivers not specified: C If you want to generate output file with rays, or with C arrival times and their errors at the receivers, C you have to specify the name of input file with C receivers. C See the main input file NET, (2). END IF IF(FREC.NE.' '.AND.FRAYS.EQ.' '.AND.FEND.EQ.' ') THEN C NET-13 CALL ERROR * ('NET-13: Name of output file with rays not specified') C NET-13: Name of output file with rays not specified: C If you specify receiver positions at the receiver C file, you should specify the name of the output file C with rays, or with arrival times and their errors at C the receivers. C See the main input file NET, (2). END IF C C (2) Number of reflections: NREFL=0 READ(LU1,*) NREFL IF(NREFL.GT.MREFL) THEN C NET-14 CALL ERROR('NET-14: Too many reflections') C NET-14: Too many reflections: C Number NREFL (main input file NET, C (2)) should be decreased to satisfy the inequality C NREFL.LE.MREFL, C or the dimension MREFL in the common block /FILES/ C in net.inc should be increased. END IF C C (3) Names of the output travel-time and predecessor files, C input velocity and index files, and input refractor-point files: DO 32 L=0,NREFL IF(L.GT.0) THEN FINTF(L)=' ' READ(LU1,*) FINTF(L) IF(FINTF(L).EQ.' ') THEN C NET-15 CALL ERROR('NET-15: Reflector file not specified') C NET-15: Reflector file not specified: C If you want to compute reflected time field, you C must specify name of reflector file. C See main input file NET, (4.2). END IF END IF FIND(L)=' ' FVEL(L)=' ' FICB(L)=' ' FTT(L)=' ' FERR(L)=' ' FPRED(L)=' ' FNFS(L)=' ' READ(LU1,*) * FIND(L),FVEL(L),FICB(L),FTT(L),FERR(L),FPRED(L),FNFS(L) IF(FVEL(L).EQ.' ') THEN C NET-16 CALL ERROR('NET-16: Velocity file not specified') C NET-16: Velocity file not specified: C Velocity file has a fundamental importance for time C field computations. You must always specify name of C velocity file. C See main input file NET, (4.1). END IF 32 CONTINUE C C End of reading main input data file 'NET' CLOSE(LU1) C C....................................................................... C C Reading source coordinates: C WRITE(*,'(A)') * '+NET: Reading source coordinates ' OPEN(LU1,FILE=FSRC,STATUS='OLD') READ(LU1,*) (TEXT,I=1,20) DO 48 ISRC=1,MSRC X1S(ISRC)=UNDEF X2S(ISRC)=UNDEF X3S(ISRC)=UNDEF TTS(ISRC)=0. TTSERR(ISRC)=0. READ(LU1,*) * NAMES,X1S(ISRC),X2S(ISRC),X3S(ISRC),TTS(ISRC),TTSERR(ISRC) IF(X1S(ISRC).EQ.UNDEF.AND.X2S(ISRC).EQ.UNDEF.AND. * X3S(ISRC).EQ.UNDEF) THEN NSRC=ISRC-1 GO TO 49 END IF IF(X1S(ISRC).EQ.UNDEF) X1S(ISRC)=0. IF(X2S(ISRC).EQ.UNDEF) X2S(ISRC)=0. IF(X3S(ISRC).EQ.UNDEF) X3S(ISRC)=0. C Checking the initial arrival time at source points IF(TTS(ISRC).LT.0.) THEN C NET-17 CALL ERROR * ('NET-17: Negative time at the source not acceptable') C NET-17: Negative time at the source not acceptable: C Initial time at the source points must not be negative C for a correct function of this code. ELSE IF(TTS(ISRC).EQ.0.) THEN TTS(ISRC)=DWARF END IF C Checking source positions IF(ISRC.EQ.1) THEN IF(NL1.EQ.1.AND..NOT.LN1) THEN X1MIN=X1S(1) X1MAX=X1S(1) END IF IF(NL2.EQ.1) THEN X2MIN=X2S(1) X2MAX=X2S(1) END IF IF(NL3.EQ.1) THEN X3MIN=X3S(1) X3MAX=X3S(1) END IF ELSE IF(NL1.EQ.1.AND..NOT.LN1) THEN IF(X1S(L).NE.X1MIN) THEN C NET-18 CALL ERROR('NET-18: Different 1-st source coordinates') C NET-18: Different 1-st source coordinates: C In a 2-D model, coordinates perpendicular to the C model plane have to be the same for all source C points and all receivers. END IF END IF IF(NL2.EQ.1) THEN IF(X2S(L).NE.X2MIN) THEN C NET-19 CALL ERROR('NET-19: Different 2-nd source coordinates') C NET-19: Different 2-nd source coordinates: C In a 2-D model, coordinates perpendicular to the C model plane have to be the same for all source C points and all receivers. END IF END IF IF(NL3.EQ.1) THEN IF(X3S(L).NE.X3MIN) THEN C NET-20 CALL ERROR('NET-20: Different 3-rd source coordinates') C NET-20: Different 3-rd source coordinates: C In a 2-D model, coordinates perpendicular to the C model plane have to be the same for all source C points and all receivers. END IF END IF END IF IF(X1S(ISRC).LT.X1MIN.OR.X1MAX.LT.X1S(ISRC).OR. * X2S(ISRC).LT.X2MIN.OR.X2MAX.LT.X2S(ISRC).OR. * X3S(ISRC).LT.X3MIN.OR.X3MAX.LT.X3S(ISRC)) THEN C NET-21 CALL ERROR('NET-21: Source is not inside the model') C NET-21: Source is not inside the model: C Source coordinates X1S,X2S,X3S must satisfy conditions C X1MIN.LE.X1S.AND.X1S.LE.X1MAX, C X2MIN.LE.X2S.AND.X2S.LE.X2MAX, C X3MIN.LE.X3S.AND.X3S.LE.X3MAX. C See the main input file NET, (2), C and input file SRC. ENDIF 48 CONTINUE C NET-22 CALL ERROR('NET-22: Too many source points') C NET-22: Too many source points: C Number of receivers in input file REC C should be less than MREC, C or the dimension MREC in the common block /POINTS/ C should be increased. 49 CONTINUE CLOSE(LU1) IF(NSRC.LE.0) THEN C NET-23 CALL ERROR('NET-23: No source point given') C NET-23: No source point given: C You must specify at least one source point at the C input source file SRC. END IF C C....................................................................... C C Reading receiver coordinates: C WRITE(*,'(A)') * '+NET: Reading receiver coordinates ' IF(FREC.EQ.' '.OR.FRAYS.EQ.' ') THEN NREC=0 ELSE OPEN(LU1,FILE=FREC,STATUS='OLD') READ(LU1,*) (TEXT,I=1,20) DO 58 IREC=1,MREC X1R(IREC)=UNDEF X2R(IREC)=UNDEF X3R(IREC)=UNDEF READ(LU1,*) NAMER(IREC),X1R(IREC),X2R(IREC),X3R(IREC) IF(X1R(IREC).EQ.UNDEF.AND.X2R(IREC).EQ.UNDEF.AND. * X3R(IREC).EQ.UNDEF) THEN NREC=IREC-1 GO TO 59 END IF IF(X1R(IREC).EQ.UNDEF) X1R(IREC)=0. IF(X2R(IREC).EQ.UNDEF) X2R(IREC)=0. IF(X3R(IREC).EQ.UNDEF) X3R(IREC)=0. C Checking the receiver positions IF(NL1.EQ.1.AND..NOT.LN1) THEN IF(X1R(IREC).NE.X1MIN) THEN C NET-24 CALL ERROR * ('NET-24: Different 1-st source and receiver coordinate') C NET-24: Different 1-st source and receiver coordinate: C In a 2-D model, coordinates perpendicular to the C model plane have to be the same for all source C points and all receivers. END IF END IF IF(NL2.EQ.1) THEN IF(X2R(IREC).NE.X2MIN) THEN C NET-25 CALL ERROR * ('NET-25: Different 2-nd source and receiver coordinate') C NET-25: Different 2-nd source and receiver coordinate: C In a 2-D model, coordinates perpendicular to the C model plane have to be the same for all source C points and all receivers. END IF END IF IF(NL3.EQ.1) THEN IF(X3R(IREC).NE.X3MIN) THEN C NET-26 CALL ERROR * ('NET-26: Different 3-rd source and receiver coordinate') C NET-26: Different 3-rd source and receiver coordinate: C In a 2-D model, coordinates perpendicular to the C model plane have to be the same for all source C points and all receivers. END IF END IF IF(X1R(IREC).LT.X1MIN.OR.X1MAX.LT.X1R(IREC).OR. * X2R(IREC).LT.X2MIN.OR.X2MAX.LT.X2R(IREC).OR. * X3R(IREC).LT.X3MIN.OR.X3MAX.LT.X3R(IREC)) THEN C NET-27 CALL ERROR('NET-27: Receiver is not inside the model') C NET-27: Receiver is not inside the model: C Receiver coordinates X1R,X2R,X3R must satisfy C conditions C X1MIN.LE.X1R.AND.X1R.LE.X1MAX, C X2MIN.LE.X2R.AND.X2R.LE.X2MAX, C X3MIN.LE.X3R.AND.X3R.LE.X3MAX. C See the main input file NET, (2), C and input file REC. ENDIF 58 CONTINUE C NET-28 CALL ERROR('NET-28: Too many receivers') C NET-28: Too many receivers: C Number of receivers in the input file c REC should be less than MREC, C or the dimension MREC in the common block /POINTS/ C in net.inc should be increased. 59 CONTINUE CLOSE(LU1) END IF C RETURN END C C======================================================================= C C C SUBROUTINE SOURCE(IREFL) INTEGER IREFL C C Initialization procedure for starting the program from the source or C from the IREFL-th interface. Reads input data for the iteration, C and initializes the arrival-time field and other arrays. C C Input: C IREFL.. Number of reflections. C C No output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ /NODE/ /SRCC/ C are required here: INCLUDE 'net.inc' C net.inc INCLUDE 'netnode.inc' C netnode.inc C C----------------------------------------------------------------------- C EXTERNAL SRCREC,TTTSRC,INDX INTEGER INDX C INTEGER LU1 PARAMETER (LU1=1) REAL GIANT PARAMETER (GIANT=1.0E+10) C INTEGER MIND,MGRID,MP2,MP3,MPRED,MNFS,MICB INTEGER N123,L123,ISRC,INDQ,NFSMAK,I,L,M INTEGER JPOS,JPOS1,JPOS2,JPOS3,JADR,NFSOPT,JCOUNT,NUMOPT REAL VMAX,V,C1,C2MIN,C2MAX,AUX,SUMOPT C C....................................................................... C WRITE(*,'(A)') '+NET: Reading velocities ' C C Reading the index array: IF(IREFL.GT.0) THEN IF(FIND(IREFL).EQ.FIND(IREFL-1)) THEN C The index array has not been changed GO TO 20 END IF END IF L4=0 IND0=0 C IF(FIND(IREFL).EQ.' ') THEN L1234=N1*N2*N3 IF(L1.GT.1.OR.L2.GT.1.OR.L3.GT.1) THEN C NET-29 CALL ERROR('NET-29: No index file specified') C NET-29: No index file specified: C If L1,L2,L3 (input data file SEP) C are specified, the index file IND has to be C submitted. See the main input file c NET, (4.1). END IF MIND=0 ELSE MIND=MRAM N123=N1*N2*N3 L123=L1*L2*L3 IF(N123.GT.MIND) THEN C NET-30 CALL ERROR('NET-30: Too many big bricks') C NET-30: Too many big bricks: C Numbers N1,N2,N3 (input data SEP) C should be decreased to satisfy the inequality C N1*N2*N3.LE.MRAM, C or the dimension MRAM of array RAM in include file C ram.inc should be C increased. END IF DO 11 L=1,N123 IRAM(IND0+L)=L 11 CONTINUE CALL RARRAI(LU1,FIND(IREFL),'FORMATTED', * .FALSE.,0,N123,IRAM(IND0+1)) DO 12 L=1,N123 L4=MAX0(IRAM(IND0+L),L4) IRAM(IND0+L)=(IRAM(IND0+L)-1)*L123+1 12 CONTINUE L1234=L1*L2*L3*L4 MIND=N123 END IF C C Dynamic array allocation: MGRID=L1234 IF(NFSMAX.GE.0) THEN C Network ray tracing: IF(FICB(IREFL).EQ.' ') THEN MICB=0 ELSE MICB=L1234 END IF IF(FNFS(IREFL).EQ.'*') THEN MNFS=0 ELSE MNFS=L1234 END IF IF(FPRED(IREFL).EQ.' '.AND. * FERR(IREFL).EQ.' '.AND.FRAYS.EQ.' '.AND.FEND.EQ.' ') THEN MPRED=0 ELSE MPRED=L1234 END IF MP2=0 MP3=0 ELSE C Second-order grid travel-time tracing: MICB=0 MNFS=0 MPRED=0 MP2=L1234 MP3=0 END IF ITT0 =IND0 +MIND IPOSQ0=ITT0 +MGRID IP0 =IPOSQ0+MGRID IP20 =IP0 +MGRID IP30 =IP20 +MP2 IPRED0=IP30 +MP3 NFS0 =IPRED0+MPRED ICB0 =NFS0 +MNFS I =ICB0 +MICB IF(I.GT.MRAM) THEN C NET-01 CALL ERROR('NET-01: Small array RAM') C NET-01: Small array RAM: C Dimension MRAM of array RAM allocated in C ram.inc C must be at least C MIND+L1*L2*L3*L4*(4+M1PRED+M1NFS+M1ICB), C where: C MIND=0 if 'IND'=' ', see 'index file' (usually C applicable), C MIND=N1*N2*N3 otherwise. C L4 is the number of nonempty big bricks. C L4=N1*N2*N3 and L1=1,L2=1,L3=1 without volume indexing C ('IND'=' '). C N1,N2,N3,L1,L2,L3 are given by input data file C SEP. C M1NFS=-1 without the optimization of sizes of forward C stars, C M1NFS=0 with the optimization of sizes of forward C stars and also for TTT, C M1PRED=0 if travel-time errors are not calculated and C the predecessor file is not generated and also C for TTT, C M1PRED=1 otherwise. C M1ICB=0 in models without structural interfaces and C also for TTT, C M1ICB=1 in models with structural interfaces. C N1,N2,N3 (input data file SEP) should C be decreased to satisfy the inequality N1*N2*N3.LE.MRAM, C or the dimension MRAM of array RAM in include file C ram.inc should be C increased. END IF C C Reading the velocity field 20 CONTINUE IF(IREFL.GT.0) THEN IF(FVEL(IREFL).EQ.FVEL(IREFL-1)) THEN C The velocity field has not been changed GO TO 30 END IF END IF IF(L1234.GT.MGRID) THEN C NET-31 CALL ERROR('NET-31: Too many network nodes') C NET-31: Too many network nodes: C The number of network nodes, see input file C SEP, should not exceed MGRID. C This error should not appear. Contact the authors. END IF IF(IREFL.GT.0) THEN IF(FERR(IREFL).NE.' ') THEN C NET-43 CALL ERROR * ('NET-43: Time errors cannot be computed after reflection') C NET-43: Time errors cannot be computed after reflection: C If the velocity field of the reflected wave is C different from the velocity field before the C reflection, the errors of the arrival times cannot be C computed correctly by this version. When the rays are C traced, only time fields and predecessors are C alternated in the memory, unlike slowness fields, see C subroutine SETERR. C Thus, if the filenames VEL(I) at the lines (4.1) of C main input data NET differ, the C filenames ERR(I) C should not be specified at the lines (4.1) for the C reflected waves. Simultaneously, if the filename C END C at the line (2) is specified, errors written in the C output file END are incorrect. END IF IF(FEND.NE.' ') THEN C NET-44 CALL ERROR * ('NET-44: Time errors cannot be computed after reflection.') C NET-44: Time errors cannot be computed after reflection: C If, in a case of reflection, filenames VEL(I) at the C lines (4.1) of main input data NET c differ and the filename END at line (2) is specified, C errors written in the output file END are incorrect. C For more details see also error NET-43 above. END IF END IF DO 21 L=1,L1234 RAM(IP0+L)=1. 21 CONTINUE CALL RARRAY(LU1,FVEL(IREFL),'FORMATTED', * .FALSE.,0.,L1234,RAM(IP0+1)) VMAX=0. AUX=1.001/GIANT DO 22 L=1,L1234 V=RAM(IP0+L) IF(ABS(V).LT.AUX) THEN RAM(IP0+L)=GIANT ELSE VMAX=AMAX1(V,VMAX) RAM(IP0+L)=1.0/V ENDIF 22 CONTINUE C C Reading the indices of blocks: 30 CONTINUE IF(IREFL.GT.0) THEN IF(FICB(IREFL).EQ.FICB(IREFL-1)) THEN C The indices of blocks have not been changed GO TO 40 END IF END IF IF(FICB(IREFL).EQ.' ') THEN LICB=.FALSE. ELSE LICB=.TRUE. IF(L1234.GT.MICB) THEN C NET-32 CALL ERROR * ('NET-32: Insufficient memory for indices of blocks') C NET-32: Insufficient memory for indices of blocks: C The number of network nodes, see input file C SEP, should not exceed MICB. C This error should not appear. Contact the authors. END IF DO 31 L=1,L1234 IRAM(ICB0+L)=1 31 CONTINUE CALL RARRAI(LU1,FICB(IREFL),'FORMATTED', * .FALSE.,0,L1234,IRAM(ICB0+1)) END IF C C Switch for storing predecessors: 40 CONTINUE IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing (array for P2): LPRED=.FALSE. IF(L1234.GT.MP2) THEN C NET-52 CALL ERROR('NET-52: Insufficient memory for P2') C NET-52: Insufficient memory for P2: C This error should not appear. Contact the authors. END IF DO 42 L=1,L1234 RAM(IP20+L)=0. 42 CONTINUE ELSE IF(FREC.EQ.' '.AND.FPRED(IREFL).EQ.' ' * .AND.FERR (IREFL).EQ.' ') THEN LPRED=.FALSE. ELSE LPRED=.TRUE. IF(L1234.GT.MPRED) THEN C NET-33 CALL ERROR('NET-33: Insufficient memory for predecessors') C NET-33: Insufficient memory for predecessors: C The number of network nodes, see input file C SEP, should not exceed MPRED. C This error should not appear. Contact the authors. END IF END IF C C Computing the forward star sizes: IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing (array for P3): LNFS=.FALSE. IF(NL3.GT.1.AND.L1234.GT.MP3) THEN C NET-53 CALL ERROR('NET-53: Insufficient memory for P3') C NET-53: Insufficient memory for P3: C This error should not appear. Contact the authors. END IF C DO 43 L=1,L1234 C RAM(IP30+L)=0. C 43 CONTINUE ELSE IF(FNFS(IREFL)(1:1).EQ.'*') THEN LNFS=.FALSE. IF(NFSMAX.LE.0) THEN C NET-08 CALL ERROR('NET-08: Zero maximum size of a forward star') C NET-08: Zero maximum size of a forward star: C For NFSMAX=0 in input data SEP, C option 'NFS(I)'='*' cannot be used in input data C NET in record (4.1). END IF ELSE LNFS=.TRUE. IF(L1234.GT.MNFS) THEN C NET-34 CALL ERROR * ('NET-34: Insufficient memory for forward star sizes') C NET-34: Insufficient memory for forward star sizes: C The number of network nodes, see input file C SEP, should not exceed MNFS. C This error should not appear. Contact the authors. END IF IF(FNFS(IREFL)(1:1).EQ.'+') THEN DO 51 L=1,L1234 IRAM(NFS0+L)=NFSMAX 51 CONTINUE CALL RARRAI(LU1,FNFS(IREFL)(2:LEN(FNFS(IREFL))),'FORMATTED', * .FALSE.,0,L1234,IRAM(NFS0+1)) IF(NFSMAX.EQ.0) THEN DO 52 L=1,L1234 IF(IRAM(NFS0+L).GT.NFSMAX) THEN NFSMAX=IRAM(NFS0+L) ENDIF 52 CONTINUE ELSE DO 53 L=1,L1234 IF(IRAM(NFS0+L).GT.NFSMAX) THEN IRAM(NFS0+L)=NFSMAX ENDIF 53 CONTINUE END IF ELSE JCOUNT=0 DO 58 JPOS3=0,NL3-1 DO 57 JPOS2=0,NL2-1 JPOS=NL1*(JPOS2+NL2*JPOS3) DO 56 JPOS1=0,NL1-1 JPOS=JPOS+1 C JPOS=1+JPOS1+NL1*(JPOS2+NL2*JPOS3) JADR=INDX(JPOS) IF(JADR.GT.0) THEN CALL OPTNFS * (JPOS,JPOS1,JPOS2,JPOS3,JADR,NFSOPT,C1,C2MIN,C2MAX) IRAM(NFS0+JADR)=NFSOPT JCOUNT=JCOUNT+1 IF(JCOUNT/1000*1000.EQ.JCOUNT) THEN WRITE(*,'(A,I7,A,I7)') * '+NET: Generating optimum sizes of forward stars:', * JCOUNT,' of',L1234 END IF END IF 56 CONTINUE 57 CONTINUE 58 CONTINUE IF(NFSMAX.EQ.0) THEN C Automatic estimation of NFSMAX: NUMOPT=0 SUMOPT=0. DO 61 L=1,L1234 I=IRAM(NFS0+L) IF(I.GT.0) THEN NUMOPT=NUMOPT+1 SUMOPT=SUMOPT+1./FLOAT(I)**2 END IF 61 CONTINUE SUMOPT=SQRT(FLOAT(NUMOPT)/SUMOPT) NFSMAX=INT(SUMOPT+0.5) NFSMAX=MIN0(NFSMAX,MAX0(N1*L1,N2*L2,N3*L3)-1) IF(NFSMAX.LE.0) THEN C NET-09 CALL ERROR('NET-09: Zero maximum size of a forward star') C NET-09: Zero maximum size of a forward star: C Contact the authors. END IF IF(NFSMAX.GT.MFSMAX) THEN C NET-11 CALL ERROR * ('NET-11: Maximum size of a forward star too great') C NET-11: Maximum size of a forward star too great: C NFSMAX in input data file SEP C should be decreased if positive, C adjusted to MFSMAX if zero, C or the dimension MFSMAX in the common block /FS/ C in net.inc should be increased. END IF DO 62 L=1,L1234 IF(IRAM(NFS0+L).GT.NFSMAX) THEN IRAM(NFS0+L)=NFSMAX END IF 62 CONTINUE END IF WRITE(*,'(A,I7,A,I7)') * '+NET: Generating optimum sizes of forward stars:', * JCOUNT,' of',L1234 IF(FNFS(IREFL).NE.' ') THEN CALL WARRAI(LU1,FNFS(IREFL)(2:LEN(FNFS(IREFL))),'FORMATTED', * .FALSE.,0,.FALSE.,0,L1234,IRAM(NFS0+1)) END IF END IF END IF C C 2-D or 3-D model: IF(NL3.EQ.1) THEN IF(NL2.EQ.1) THEN IF(NL1.EQ.1) THEN C NET-35 CALL ERROR('NET-35: One-point model') C NET-35: One-point model: C N1*L1=N2*L2=N3*L3=1 in input file C SEP. ELSE C NET-36 CALL ERROR('NET-36: Line model') C NET-36: Line model: C N2*L2=N3*L3=1 in input file SEP. END IF ELSE IF(NL1.EQ.1) THEN C NET-37 CALL ERROR('NET-37: Line model') C NET-37: Line model: C N1*L1=N3*L3=1 in input file SEP. ELSE DMIN=AMIN1(ASX1,ASX2)/VMAX END IF END IF ELSE IF(NL2.EQ.1) THEN IF(NL1.EQ.1) THEN C NET-38 CALL ERROR('NET-38: Line model') C NET-38: Line model: C N1*L1=N2*L2=1 in input file SEP. ELSE DMIN=AMIN1(ASX1,ASX3)/VMAX END IF ELSE IF(NL1.EQ.1) THEN DMIN=AMIN1(ASX2,ASX3)/VMAX ELSE DMIN=AMIN1(ASX1,ASX2,ASX3)/VMAX END IF END IF END IF IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing: AUX=ASX1*ASX1+ASX2*ASX2+ASX3*ASX3 DMIN=SQRT(AUX)-SQRT(AUX-DMIN*DMIN) END IF C C....................................................................... C C Writing parameters on the screen: IF(IREFL.EQ.0) THEN WRITE(*,'(A,I2,A,3(I4,A),I9,A,I2,A,I4,A)') * '+',NREFL,' reflections,',NL1,'*',NL2,'*',NL3,'=',NL1*NL2*NL3, * ' gridpoints, f.s.size:',NFSMAX,',',NREC,' receivers' WRITE(*,*) END IF C C....................................................................... C C Initialization: C WRITE(*,'(A)') * '+NET: Initializing source nodes ' C IF(IREFL.LE.0) THEN C C Given source points: C C Initialization of arrays DO 71 I=1,NL1*NL2*NL3 IADR=INDX(I) IF(IADR.GT.0) THEN RAM(ITT0+IADR)=GIANT IF(LPRED) THEN IRAM(IPRED0+IADR)=I END IF END IF 71 CONTINUE C C MINQ,MAXQ describe the extent of the queue. MINQ=1 MAXQ=0 TTMIN=TTS(1) C C Extent of a dense forward star at source points: IF(LNFS) THEN NFSMAK=0 ELSE IF(NFSMAX.GE.0) THEN NFSMAK=NFSMAX CALL MAKEFS(NFSMAK) END IF C C Loop over source points DO 79 ISRC=1,NSRC C C Minimum first-arrival time: TTMIN=AMIN1(TTS(ISRC),TTMIN) C C Parameters of the source point: CALL SLOW(X1S(ISRC),X2S(ISRC),X3S(ISRC),DPOS1,DPOS2,DPOS3, * IPOS1,IPOS2,IPOS3,IPOS,IADR,PI) TTI=TTS(ISRC) ISRCI=-ISRC C C Updating: IF(NFSMAX.GE.0) THEN C Network ray tracing: C Adjusting extent of dense forward star at the source point IF(LNFS) THEN IF(IRAM(NFS0+IADR).NE.NFSMAK) THEN NFSMAK=IRAM(NFS0+IADR) CALL MAKEFS(NFSMAK) END IF END IF C Updating CALL LOOP(SRCREC) ELSE C Second-order grid travel-time tracing: CALL TTTSRC() END IF C 79 CONTINUE C C Reading and assembling optimized forward star: CALL READFS() C C....................................................................... C ELSE C C Given reflector: C C Reading the interface file: C Creating queue for travel-times on the interface DO 81 INDQ=1,MGRID IRAM(IPOSQ0+INDQ)=0 81 CONTINUE CALL RARRAI(LU1,FINTF(IREFL),'FORMATTED', * .FALSE.,0,L1234,IRAM(IPOSQ0+1)) DO 82 INDQ=1,MGRID IF(IRAM(IPOSQ0+INDQ).LE.0) THEN MAXQ=INDQ-1 GO TO 83 END IF 82 CONTINUE C NET-39 CALL ERROR('NET-39: Too many points of reflector') C NET-39: Too many points of reflector: C Number of reflector points should be less than the C number of network nodes (gridpoints). 83 CONTINUE C C Labeling the time field in the queue (TT(I)=RAM(ITT0+I).LT.0): M=0 TTMIN=GIANT DO 84 I=1,MAXQ IADR=INDX(IRAM(IPOSQ0+I)) C Check for the computational volume and for a free space: IF(IADR.GT.0) THEN IF(RAM(IP0+IADR).LT.GIANT) THEN M=M+1 IRAM(IPOSQ0+M)=IRAM(IPOSQ0+I) TTMIN=AMIN1(RAM(ITT0+IADR),TTMIN) RAM(ITT0+IADR)=-RAM(ITT0+IADR) END IF END IF 84 CONTINUE MINQ=1 MAXQ=M C C Defining the 1st approximation of the s.p.t. And initially C updating the interface: DO 88 I=1,NL1*NL2*NL3 IADR=INDX(I) IF(IADR.GT.0) THEN IF(RAM(ITT0+IADR).LT.0..AND.RAM(IP0+IADR).LT.GIANT) THEN RAM(ITT0+IADR)=-RAM(ITT0+IADR) IF(LPRED) THEN IRAM(IPRED0+IADR)=0 ENDIF ELSE RAM(ITT0+IADR)=GIANT IF(LPRED) THEN IRAM(IPRED0+IADR)=I ENDIF ENDIF ENDIF 88 CONTINUE C ENDIF IF(MAXQ.LT.MINQ) THEN C NET-50 CALL ERROR('NET-50: No source point') C NET-50: No source point: C This error should not appear. Contact the authors. ENDIF RETURN END C C======================================================================= C C C SUBROUTINE GENER(IREFL) INTEGER IREFL C C Procedure generating travel-time field and predecessors by performing C shortest path ray tracing (one iteration). C C Input: C IREFL.. Number of reflections. C C No output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NODE/ are required here: INCLUDE 'net.inc' C net.inc INCLUDE 'netnode.inc' C netnode.inc C C----------------------------------------------------------------------- C REAL GIANT PARAMETER (GIANT=1.E+10) C EXTERNAL INDX,UPDATE,TTTUPD,FREVOL INTEGER INDX C INTEGER NUMNOD,NUMA,NUMINC,NUMOLD,INDQ,I REAL TTINDQ PARAMETER (NUMINC=1000) C C NUMNOD..Number of network nodes not situated in a free space. C NUMA... Counter of finished nodes (nodes moved to set A). C NUMINC..Minimum increment of the number of finished nodes between C the reports on the screen. If the increment is smaller, C screen output is suppressed. C NUMOLD..Counter of finished nodes reported on the screen. C INDQ... Loop variable (index in Q). C I... Loop variable (index of a node). C TTINDQ..Travel time in the INDQ-th Q element. C C....................................................................... C C Check for a free space: WRITE(*,'(A)') * '+NET: Eliminating nodes in a free space. ' NUMNOD=L1234 DO 8 I=1,L1234 IF(RAM(IP0+I).GE.GIANT) THEN RAM(ITT0+I)=-GIANT NUMNOD=NUMNOD-1 ENDIF 8 CONTINUE C C Counter of finished nodes (nodes moved to set A): NUMA=0 NUMOLD=-NUMINC C C Loop for intervals: 10 CONTINUE C New interval: TTMIN=TTMIN+DMIN C C Determination of the first element MINQ of set B in Q: DO 21 INDQ=MINQ,MAXQ TTINDQ=RAM(ITT0+INDX(IRAM(IPOSQ0+INDQ))) IF(0..LT.TTINDQ) THEN MINQ=INDQ GO TO 22 END IF 21 CONTINUE MINQ=MAXQ+1 22 CONTINUE C C Screen output: IF(NUMA.GE.NUMOLD+NUMINC.OR.MINQ.GT.MAXQ) THEN WRITE(*,'(A,I2,4(A,I7))') * '+',IREFL,'-th reflection',NUMA,' nodes of',NUMNOD, * ' finished, QMIN=',MINQ,', QMAX=',MAXQ NUMOLD=NUMA END IF C C (4) Iteration check: testing the end of time field generation C condition to finish the generation of the s.p.t. IF(MINQ.GT.MAXQ) THEN DO 31 IADR=1,L1234 RAM(ITT0+IADR)=ABS(RAM(ITT0+IADR)) 31 CONTINUE RETURN ENDIF C DO 40 INDQ=MINQ,MAXQ C C (2) Selection: TTINDQ=RAM(ITT0+INDX(IRAM(IPOSQ0+INDQ))) IF(0..LT.TTINDQ.AND.TTINDQ.LT.TTMIN) THEN C C New node 'I' (position IPOS, address IADR): IPOS=IRAM(IPOSQ0+INDQ) IPOS1=IPOS-1 IPOS3=IPOS1/(NL1*NL2) IPOS2=IPOS1/NL1-IPOS3*NL2 IPOS1=IPOS1-(IPOS2+IPOS3*NL2)*NL1 IADR=INDX(IPOS) PI=RAM(IP0+IADR) TTI=RAM(ITT0+IADR) C C (3) Updating nodes 'J' of the forward star FS(I), which are C the neighbours to the node 'I': IF(NFSMAX.GE.0) THEN C Network ray tracing: IF(L4.EQ.0) THEN CALL LOOP(UPDATE) ELSE CALL LOOP(FREVOL) END IF ELSE C Second-order grid travel-time tracing: IF(L4.EQ.0) THEN CALL TTTUPD() ELSE C NET-51 CALL ERROR('NET-51: Fresnel volumes disabled') C NET-51: Fresnel volumes disabled: C Second-order grid travel-time tracing cannot be C performed in Fresnel or other volumes specified C by means of the index file. END IF END IF C C Moving 'I' from set B to set A: IADR=INDX(IPOS) RAM(ITT0+IADR)=-RAM(ITT0+IADR) C Increment of the counter (number of nodes in the set A): NUMA=NUMA+1 C END IF 40 CONTINUE C End of loop for nodes 'I', which are secondary sources of the C time field. C GO TO 10 C End of loop for intervals. END C C======================================================================= C C C SUBROUTINE RESFIL(IREFL) INTEGER IREFL C C In a case of reflections, stores the results of an iteration in an C unformatted direct-access scratch file (one iteration). C C Input: C IREFL.. Number of reflections. C C No output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FILES/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER LU0,LU PARAMETER(LU0=10) C INTEGER I C C I... Temporary loop variable. C C....................................................................... C WRITE(*,'(A)') * ' NET: Writing travel time field ' C IF(IREFL.LT.NREFL) THEN LU=LU0+IREFL OPEN(LU,RECL=8,FORM='UNFORMATTED',ACCESS='DIRECT', * STATUS='SCRATCH') IF(LPRED) THEN DO 11 I=1,L1234 WRITE(LU,REC=I) RAM(ITT0+I),IRAM(IPRED0+I) 11 CONTINUE ELSE DO 12 I=1,L1234 WRITE(LU,REC=I) RAM(ITT0+I) 12 CONTINUE END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE ERRORS(IREFL) INTEGER IREFL C C Procedure generating the travel-time errors at all network nodes C coinciding with gridpoints by means of backward ray tracing. C Writes 'ERR(I)' (one iteration). C C Input: C IREFL.. Number of reflections. C C No output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /FILES/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER LUERR PARAMETER (LUERR=5) C EXTERNAL INDX INTEGER INDX LOGICAL LTRACE,LRAYS,LEND INTEGER IPREDE,IEND1,IEND2,IEND3,IEND,IENDA,IAUX,JCOUNT REAL DUMMY C C ERR(I)=RAM(IERR0+I): INTEGER IERR0 EQUIVALENCE (IERR0,IPOSQ0) C C....................................................................... C IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing: Errors are not generated RETURN END IF C WRITE(*,'(A)') * '+NET: Generating travel-time errors. ' C IF(FERR(IREFL).NE.' ') THEN LTRACE=.FALSE. LRAYS=.FALSE. LEND=.FALSE. C DO 10 IEND=1,L1234 RAM(IERR0+IEND)=-1. 10 CONTINUE C JCOUNT=0 DO 23 IEND3=0,NL3-1 DO 22 IEND2=0,NL2-1 IEND=NL1*(IEND2+NL2*IEND3) DO 21 IEND1=0,NL1-1 IEND=IEND+1 C IEND=1+IEND1+NL1*(IEND2+NL2*IEND3) IENDA=INDX(IEND) IF(IENDA.GT.0) THEN IF(RAM(IERR0+IENDA).LT.0.) THEN IPREDE=IRAM(IPRED0+IENDA) IF(IPREDE.NE.IEND) THEN CALL ONERAY(LTRACE,LRAYS,LEND,IAUX,IREFL,IPREDE, * IEND1,IEND2,IEND3,IENDA,RAM(ITT0+IENDA), * RAM(IERR0+1),DUMMY,IAUX) END IF END IF JCOUNT=JCOUNT+1 IF(JCOUNT/1000*1000.EQ.JCOUNT) THEN WRITE(*,'(A,I7,A,I7)') * '+NET: Generating travel-time errors:', * JCOUNT,' of',L1234 END IF END IF 21 CONTINUE 22 CONTINUE 23 CONTINUE C C....................................................................... C C Writing the travel-time error estimates: C WRITE(*,'(A)') * '+NET: Writing travel-time errors. ' C CALL WARRAY(LUERR,FERR(IREFL),'FORMATTED', * .TRUE.,-1.,.FALSE.,0.,L1234,RAM(IERR0+1)) C END IF RETURN END C C======================================================================= C C C SUBROUTINE RECVRS() C C Procedure updating arrival times at the receivers. C C No input, no output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ /NODE/ /SRCC/ C are required here: INCLUDE 'net.inc' C net.inc INCLUDE 'netnode.inc' C netnode.inc C C----------------------------------------------------------------------- C EXTERNAL SRCREC C REAL GIANT PARAMETER (GIANT=1.E+10) C INTEGER IR,ISRC,NFSREC,NFSSRC,NFSMAK,I1,I2,I3,I4,I5 REAL DIM,PS,A1,A2,A3,DIST2,TTSR C C IR... Index of the receiver. C ISRC... Index of the source. C NFSREC..Size of the f.s. at the receiver. C NFSSRC..Size of the f.s. at the source. C NFSMAK..Size of recently generated full forward star at source C points. C DIM... 2 or 3 for 2-D or 3-D calculation, respectively. C PS... Slowness at the source. C A1,A2,A3,I1,I2,I3,I4,I5... Other source parameters. C DIST2.. Square of the source-reiver distance. C TTSR... Arrival time at the receiver from the source under C consideration. C C....................................................................... C IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing: IF(NREC.GT.0) THEN WRITE(*,'(A)') *'+NET: Receivers are not supported by the 2-D second-order method' WRITE(*,'(A)') ' ' END IF RETURN END IF C C Computing approximate arrival time and rays at the receiver nodes: WRITE(*,'(A)') * '+NET: Updating receivers... ' C C Extent of a dense forward star at receiver points: IF(LNFS) THEN NFSMAK=0 ELSE NFSMAK=NFSMAX CALL MAKEFS(NFSMAK) END IF C C Loop over receivers: DO 90 IR=1,NREC C C Initially, receiver is considered not connected: IPREDR(IR)=0 C C Parameters of the receiver point: CALL SLOW(X1R(IR),X2R(IR),X3R(IR),DPOS1,DPOS2,DPOS3, * IPOS1,IPOS2,IPOS3,IPOS,IADR,PI) TTI=GIANT ISRCI=IR C C....................................................................... C C For no reflections, testing the source-receiver position: C IF(NREFL.EQ.0) THEN C Loop over source points: DO 39 ISRC=1,NSRC C C Square of the source-reiver distance: DIST2=0. IF(LN1) THEN DIST2=DIST2+((X1S(ISRC)-X1R(IR))**2 * +(X2S(ISRC)-X2R(IR))**2)/ASX2**2 END IF IF(NL1.GT.1) THEN DIST2=DIST2+((X1S(ISRC)-X1R(IR))/ASX1)**2 END IF IF(NL2.GT.1) THEN DIST2=DIST2+((X2S(ISRC)-X2R(IR))/ASX2)**2 END IF IF(NL3.GT.1) THEN DIST2=DIST2+((X3S(ISRC)-X3R(IR))/ASX3)**2 END IF C C Size of the f.s. at the receiver point: IF(LNFS) THEN NFSREC=IRAM(NFS0+IADR) ELSE NFSREC=NFSMAX END IF C C Dimension of the model and forward stars (2 or 3): IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN DIM=2. ELSE DIM=3. END IF C IF(DIST2.LE.FLOAT(NFSREC)**2+DIM-1.) THEN C Source is situated within the receiver f.s. ellipsoid CALL SLOW(X1S(ISRC),X2S(ISRC),X3S(ISRC),A1,A2,A3, * I1,I2,I3,I4,I5,PS) C C Check for crossing an interface: IF(LICB) THEN IF(IRAM(ICB0+IADR).NE.IRAM(ICB0+I5)) THEN IF(DIST2.GT.DIM) GO TO 38 END IF END IF C C Size of the f.s. at the source point: IF(LNFS) THEN NFSSRC=IRAM(NFS0+I5) ELSE NFSSRC=NFSMAX END IF C C Updating: IF(DIST2.LE.FLOAT(NFSSRC)**2+DIM-1.) THEN C Travel time from source to the receiver TTSR=SQRT((X1R(IR)-X1S(ISRC))**2+ * (X2R(IR)-X2S(ISRC))**2+ * (X3R(IR)-X3S(ISRC))**2)*0.5*(PI+PS)+TTS(ISRC) C Updating IF(IPREDR(IR).LT.0) THEN IF(TTR(IR).GT.TTSR) THEN IPREDR(IR)=-ISRC TTR(IR)=TTSR END IF ELSE IPREDR(IR)=-ISRC TTR(IR)=TTSR END IF END IF C 38 CONTINUE END IF 39 CONTINUE END IF C C....................................................................... C C Connecting the receiver still not connected: C IF(IPREDR(IR).EQ.0) THEN C Adjusting extent of the dense forward star at the receiver C point: IF(LNFS) THEN IF(IRAM(NFS0+IADR).NE.NFSMAK) THEN NFSMAK=IRAM(NFS0+IADR) CALL MAKEFS(NFSMAK) END IF END IF C C Updating: CALL LOOP(SRCREC) TTR(IR)=TTI ENDIF C 90 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE TRACER() C C Procedure performing backward ray tracing from the receivers, C including the estimation of arrival-time errors. C C No input, no output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /FILES/ are required: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER LURAYS,LUEND PARAMETER (LURAYS=4) PARAMETER (LUEND=5) C LOGICAL LTRACE,LRAYS,LEND INTEGER IPREDE,IEND1,IEND2,IEND3,IAUX,IENDA,IREC,ISRC REAL DUMMY1,DUMMY2,DUMMY3,PEND,TERR,DUMMY C C IREC... Loop variable - index of a receiver. C C....................................................................... C IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing: No posterior ray tracing RETURN END IF C LTRACE=.TRUE. WRITE(*,'(A)') * '+NET: Backward ray tracing... ' IF(FRAYS.EQ.' ') THEN LRAYS=.FALSE. ELSE LRAYS=.TRUE. END IF IF(FEND.EQ.' ') THEN LEND=.FALSE. ELSE LEND=.TRUE. END IF C IF(LRAYS) THEN OPEN(LURAYS,FILE=FRAYS) WRITE(LURAYS,'(A)') '/' END IF IF(LEND) THEN OPEN(LUEND,FILE=FEND) WRITE(LUEND,'(A)') '/' END IF DO 80 IREC=1,NREC IPREDE=IPREDR(IREC) IF(IPREDE.NE.0) THEN IF(LRAYS) THEN IF(NAMES.EQ.' ') THEN IF(NAMER(IREC).EQ.' ') THEN WRITE(LURAYS,'(2A,I5,5A)') * '''','RAY',IREC, ''' /' ELSE WRITE(LURAYS,'(2A,I5,5A)') * '''','RAY',IREC,' ',NAMES,' TO ',NAMER(IREC),''' /' END IF ELSE IF(NAMER(IREC).EQ.' ') THEN WRITE(LURAYS,'(2A,I5,5A)') * '''','RAY',IREC,' FROM ',NAMES, ''' /' ELSE WRITE(LURAYS,'(2A,I5,5A)') * '''','RAY',IREC,' FROM ',NAMES,' TO ',NAMER(IREC),''' /' END IF END IF WRITE(LURAYS,'(4F12.6,A)') * X1R(IREC),X2R(IREC),X3R(IREC),TTR(IREC),' /' END IF IF(LEND) THEN CALL SLOW(X1R(IREC),X2R(IREC),X3R(IREC), * DUMMY1,DUMMY2,DUMMY3,IEND1,IEND2,IEND3,IAUX,IENDA,PEND) END IF CALL ONERAY(LTRACE,LRAYS,LEND,LURAYS,NREFL, * IPREDE,IEND1,IEND2,IEND3,IENDA,TTR(IREC),DUMMY,TERR,ISRC) IF(LRAYS) THEN WRITE(LURAYS,'(4F12.6,A)') * X1S(ISRC),X2S(ISRC),X3S(ISRC),TTS(ISRC),' /' WRITE(LURAYS,'(''/'')') END IF IF(LEND) THEN WRITE(LUEND,'(3A,5F12.6,A)') '''',NAMER(IREC),'''', * X1R(IREC),X2R(IREC),X3R(IREC),TTR(IREC),TERR,' /' END IF END IF 80 CONTINUE IF(LRAYS) THEN WRITE(LURAYS,'(''/'')') CLOSE(LURAYS) END IF IF(LEND) THEN WRITE(LUEND,'(''/'')') CLOSE(LUEND) END IF RETURN END C C======================================================================= C C C SUBROUTINE OUTPUT() C C Subroutine to rewrite arrival-time fields and predecessors from memory C or scratch files to formatted output files 'TT(I)' and 'IPRED(I)'. C C No input, no output. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /FILES/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER LU0,LU1,LU PARAMETER(LU0=10) PARAMETER(LU1=1) REAL GIANT PARAMETER (GIANT=1.0E+10) C INTEGER NF,I C NF,I... Loop variables. C C....................................................................... C WRITE(*,'(A)') * '+NET: Generating output files... ' C DO 30 NF=NREFL,0,-1 C C Reading scratch file: IF(NF.NE.NREFL) THEN LU=LU0+NF IF(FTT(NF).NE.' '.OR.FPRED(NF).NE.' ') THEN IF(LPRED) THEN DO 11 I=1,L1234 READ(LU,REC=I) RAM(ITT0+I),IRAM(IPRED0+I) 11 CONTINUE ELSE DO 12 I=1,L1234 READ(LU,REC=I) RAM(ITT0+I) 12 CONTINUE END IF END IF C Closing the scratch file C ********* CLOSE(LU) C ********* END IF C C Writing arrival times: IF(FTT(NF).NE.' ') THEN CALL WARRAY(LU1,FTT(NF),'FORMATTED', * .FALSE.,0.,.TRUE.,GIANT,L1234,RAM(ITT0+1)) END IF C C Writing predecessors: IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing C (writing slowness vector): IF(N1.GT.1.AND.FERR(NF).NE.' ') THEN CALL WARRAY(LU1,FERR(NF),'FORMATTED', * .FALSE.,0.,.FALSE.,GIANT,L1234,RAM(IP0+1)) END IF IF(N2.GT.1.AND.FPRED(NF).NE.' ') THEN CALL WARRAY(LU1,FPRED(NF),'FORMATTED', * .FALSE.,0.,.FALSE.,GIANT,L1234,RAM(IP20+1)) END IF C IF(N3.GT.1.AND.FNFS(NF).NE.' ') THEN C CALL WARRAY(LU1,FNFS(NF),'FORMATTED', C * .FALSE.,0.,.FALSE.,GIANT,L1234,RAM(IP30+1)) C END IF ELSE IF(FPRED(NF).NE.' ') THEN CALL WARRAI(LU1,FPRED(NF),'FORMATTED', * .FALSE.,0,.FALSE.,0,L1234,IRAM(IPRED0+1)) END IF 30 CONTINUE C WRITE(*,'(A)') * '+NET: Done. ' RETURN END C C======================================================================= C C C SUBROUTINE LOOP(UPDATE) EXTERNAL UPDATE C C This subroutine performs the loop over the nodes 'J' of the forward C star FS(I). A clear but exact arrangement of the loop is expressed C inside asterisks. The code following the asterisk rectangular is C equivalent, but much longer and somewhat faster. C Auxiliary subroutine to SOURCE, GENER, and RECVRS. C LOOP(FREVOL) or LOOP(UPDATE): employed by GENER, C LOOP(SRCREC): employed by SOURCE and RECVRS. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NODE/ are required here: INCLUDE 'net.inc' C net.inc INCLUDE 'netnode.inc' C netnode.inc C C----------------------------------------------------------------------- C INTEGER I1MIN,I1MAX,I2MIN,I2MAX,I3MIN,I3MAX,I1,I2,I3 INTEGER NFSOLD,NLEX1,NLEX2,NLEX3,IFSMIN,IFSMAX SAVE NFSOLD,NLEX1,NLEX2,NLEX3,IFSMIN,IFSMAX C Still no forward star used: DATA NFSOLD/0/ C C I1MIN,I1MAX,I2MIN,I2MAX,I3MIN,I3MAX... C I1,I2,I3... C NFSOLD... C NLEX1,NLEX2,NLEX3... C IFSMIN,IFSMAX... C C....................................................................... C C Additional check of input data IF(NFSMAX.LT.0) THEN C Second-order grid travel-time tracing -- inaccessible branch: C NET-55 CALL ERROR('NET-55: LOOP should not be called') C NET-55: LOOP should not be called: C This error should not appear. Contact the authors. END IF C C Index of the block (if applicable): IF(LICB) THEN ICBI=IRAM(ICB0+IADR) END IF C C Adjusting extent of the forward star: IF(NFSOLD.EQ.0) THEN NFSOLD=NFSMAX NLEX1=MIN0(NFSOLD,NL1-1) NLEX2=MIN0(NFSOLD,NL2-1) NLEX3=MIN0(NFSOLD,NL3-1) END IF IF(LNFS) THEN IF(IRAM(NFS0+IADR).NE.NFSOLD) THEN NFSOLD=IRAM(NFS0+IADR) NLEX1=MIN0(NFSOLD,NL1-1) NLEX2=MIN0(NFSOLD,NL2-1) NLEX3=MIN0(NFSOLD,NL3-1) END IF END IF C C Location of the forward star in the memory: IFSMIN=KFS0(NFSOLD-1)+1 IFSMAX=KFS0(NFSOLD) C C Extent of the intersection of the forward star with the grid: I1MIN= -IPOS1 I1MAX=NL1-1-IPOS1 I2MIN= -IPOS2 I2MAX=NL2-1-IPOS2 I3MIN= -IPOS3 I3MAX=NL3-1-IPOS3 C ************************************************************************ * DO 163 IFS=IFSMIN,IFSMAX * * I1=KFS1(IFS) * * I2=KFS2(IFS) * * I3=KFS3(IFS) * * IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * * * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * * * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN * * CALL UPDATE() * * END IF * * 163 CONTINUE * ************************************************************************ C C This is an optimized version of the above loop: IF(I3MIN.LE.-NLEX3) THEN IF(NLEX3.LE.I3MAX) THEN IF(I2MIN.LE.-NLEX2) THEN IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 100 IFS=IFSMIN,IFSMAX CALL UPDATE() 100 CONTINUE ELSE DO 101 IFS=IFSMIN,IFSMAX IF(KFS1(IFS).LE.I1MAX) THEN CALL UPDATE() END IF 101 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 102 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS)) THEN CALL UPDATE() END IF 102 CONTINUE ELSE DO 103 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX) THEN CALL UPDATE() END IF 103 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 104 IFS=IFSMIN,IFSMAX IF(KFS2(IFS).LE.I2MAX) THEN CALL UPDATE() END IF 104 CONTINUE ELSE DO 105 IFS=IFSMIN,IFSMAX IF(KFS1(IFS).LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX) THEN CALL UPDATE() END IF 105 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 106 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * KFS2(IFS).LE.I2MAX) THEN CALL UPDATE() END IF 106 CONTINUE ELSE DO 107 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX) THEN CALL UPDATE() END IF 107 CONTINUE END IF END IF END IF ELSE IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 108 IFS=IFSMIN,IFSMAX IF(I2MIN.LE.KFS2(IFS)) THEN CALL UPDATE() END IF 108 CONTINUE ELSE DO 109 IFS=IFSMIN,IFSMAX IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS)) THEN CALL UPDATE() END IF 109 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 110 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.KFS2(IFS)) THEN CALL UPDATE() END IF 110 CONTINUE ELSE DO 111 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS)) THEN CALL UPDATE() END IF 111 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 112 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF(I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN CALL UPDATE() END IF 112 CONTINUE ELSE DO 113 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN CALL UPDATE() END IF 113 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 114 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN CALL UPDATE() END IF 114 CONTINUE ELSE DO 115 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I2=KFS2(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX) THEN CALL UPDATE() END IF 115 CONTINUE END IF END IF END IF END IF ELSE IF(I2MIN.LE.-NLEX2) THEN IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 116 IFS=IFSMIN,IFSMAX IF(KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 116 CONTINUE ELSE DO 117 IFS=IFSMIN,IFSMAX IF(KFS1(IFS).LE.I1MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 117 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 118 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 118 CONTINUE ELSE DO 119 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 119 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 120 IFS=IFSMIN,IFSMAX IF(KFS2(IFS).LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 120 CONTINUE ELSE DO 121 IFS=IFSMIN,IFSMAX IF(KFS1(IFS).LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 121 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 122 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * KFS2(IFS).LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 122 CONTINUE ELSE DO 123 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 123 CONTINUE END IF END IF END IF ELSE IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 124 IFS=IFSMIN,IFSMAX IF(I2MIN.LE.KFS2(IFS).AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 124 CONTINUE ELSE DO 125 IFS=IFSMIN,IFSMAX IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS).AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 125 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 126 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.KFS2(IFS).AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 126 CONTINUE ELSE DO 127 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS).AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 127 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 128 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF(I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 128 CONTINUE ELSE DO 129 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 129 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 130 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 130 CONTINUE ELSE DO 131 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I2=KFS2(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * KFS3(IFS).LE.I3MAX) THEN CALL UPDATE() END IF 131 CONTINUE END IF END IF END IF END IF END IF ELSE IF(NLEX3.LE.I3MAX) THEN IF(I2MIN.LE.-NLEX2) THEN IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 132 IFS=IFSMIN,IFSMAX IF(I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 132 CONTINUE ELSE DO 133 IFS=IFSMIN,IFSMAX IF( KFS1(IFS).LE.I1MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 133 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 134 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 134 CONTINUE ELSE DO 135 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 135 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 136 IFS=IFSMIN,IFSMAX IF( KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 136 CONTINUE ELSE DO 137 IFS=IFSMIN,IFSMAX IF( KFS1(IFS).LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 137 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 138 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 138 CONTINUE ELSE DO 139 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 139 CONTINUE END IF END IF END IF ELSE IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 140 IFS=IFSMIN,IFSMAX IF(I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 140 CONTINUE ELSE DO 141 IFS=IFSMIN,IFSMAX IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 141 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 142 IFS=IFSMIN,IFSMAX IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 142 CONTINUE ELSE DO 143 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 143 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 144 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF(I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 144 CONTINUE ELSE DO 145 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 145 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 146 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 146 CONTINUE ELSE DO 147 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I2=KFS2(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.KFS3(IFS)) THEN CALL UPDATE() END IF 147 CONTINUE END IF END IF END IF END IF ELSE IF(I2MIN.LE.-NLEX2) THEN IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 148 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF(I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 148 CONTINUE ELSE DO 149 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 149 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 150 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 150 CONTINUE ELSE DO 151 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I3=KFS3(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 151 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 152 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF( KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 152 CONTINUE ELSE DO 153 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 153 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 154 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 154 CONTINUE ELSE DO 155 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I3=KFS3(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * KFS2(IFS).LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 155 CONTINUE END IF END IF END IF ELSE IF(NLEX2.LE.I2MAX) THEN IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 156 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF(I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 156 CONTINUE ELSE DO 157 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 157 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 158 IFS=IFSMIN,IFSMAX I3=KFS3(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 158 CONTINUE ELSE DO 159 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I3=KFS3(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.KFS2(IFS).AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 159 CONTINUE END IF END IF ELSE IF(I1MIN.LE.-NLEX1) THEN IF(NLEX1.LE.I1MAX) THEN DO 160 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) I3=KFS3(IFS) IF(I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 160 CONTINUE ELSE DO 161 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) I3=KFS3(IFS) IF( KFS1(IFS).LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 161 CONTINUE END IF ELSE IF(NLEX1.LE.I1MAX) THEN DO 162 IFS=IFSMIN,IFSMAX I2=KFS2(IFS) I3=KFS3(IFS) IF(I1MIN.LE.KFS1(IFS).AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 162 CONTINUE ELSE DO 163 IFS=IFSMIN,IFSMAX I1=KFS1(IFS) I2=KFS2(IFS) I3=KFS3(IFS) IF(I1MIN.LE.I1.AND.I1.LE.I1MAX.AND. * I2MIN.LE.I2.AND.I2.LE.I2MAX.AND. * I3MIN.LE.I3.AND.I3.LE.I3MAX) THEN CALL UPDATE() END IF 163 CONTINUE END IF END IF END IF END IF END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE FREVOL() C C Procedure determining index of the node 'J', if its offset C respectively to the node 'I' is given. C This procedure has to be called if the indices of nodes differ from C the indices of gridpoints. C Called by LOOP if LOOP is called by GENER. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NODE/ are required here: INCLUDE 'net.inc' C net.inc INCLUDE 'netnode.inc' C netnode.inc C C----------------------------------------------------------------------- C REAL GIANT PARAMETER (GIANT=1.E+10) C INTEGER JPOS1,JPOS2,JPOS3,JN1,JN2,JN3,JL1,JL2,JL3,JADR,J C C JPOS1,JPOS2,JPOS3... Position of the gridpoint 'J' within the C model. JPOSI=0,1,2,...,NLI-1. C JN1,JN2,JN3... Position of the big brick with 'J' within the C model. JNI=0,1,2,...,NI-1. C JL1,JL2,JL3... Position of the gridpoint 'J' within the big brick. C JLI=0,1,2,...,LI-1. C JADR... Index of the node 'J'. C J... Index of the first node of the big brick. C REAL TTJ,TTIJ C C TTJ... Arrival time at the node 'J'. C TTIJ... Arrival time through the node 'I' to the node 'J'. C C....................................................................... C C 22 integer operations JPOS1=IPOS1+KFS1(IFS) JPOS2=IPOS2+KFS2(IFS) JPOS3=IPOS3+KFS3(IFS) JN1=JPOS1/L1 JN2=JPOS2/L2 JN3=JPOS3/L3 J=IRAM(IND0+1+JN1+N1*(JN2+N2*JN3)) IF(J.LE.0) RETURN JL1=JPOS1-JN1*L1 JL2=JPOS2-JN2*L2 JL3=JPOS3-JN3*L3 JADR=J+JL1+L1*(JL2+L2*JL3) GO TO 20 C C======================================================================= C C C ENTRY UPDATE() C C Procedure updating travel time at node 'J' of the network. C Called by LOOP if LOOP is called by GENER. C C----------------------------------------------------------------------- C C Index of the node 'J': JADR=IPOS+KFS4(IFS) C 20 CONTINUE TTJ=RAM(ITT0+JADR) C IF(TTJ.LT.0.) RETURN C Node 'J' is not in the set A, may be updated. C C Check for crossing an interface: IF(LICB) THEN IF(IRAM(ICB0+JADR).NE.ICBI) THEN IF(KFS5(IFS).GT.1) RETURN END IF END IF C C*********************************************************************** C Consolidating forward and backward stars (optional): * IF(LNFS) THEN * IF(KFS5(IFS).GT.IRAM(NFS0+JADR)) RETURN * END IF C To guarantee the coincidence of forward and backward stars, it C is necessary not only to enable the above three statements, but C also to submit such template forward stars that each template C forward star is a subset of the forward stars of greater sizes. C*********************************************************************** C C Arrival time from the node 'I' to the node 'J': TTIJ=TTI+DFS(IFS)*(RAM(IP0+JADR)+PI) C C Test of Bellman's condition for shortest path IF(TTIJ.GE.TTJ) RETURN C C Testing if the node 'J' is in the queue IF(TTJ.GE.GIANT) THEN C Updating queue by adding 'J' as Q(MAXQ) MAXQ=MAXQ+1 IRAM(IPOSQ0+MAXQ)=IPOS+KFS4(IFS) END IF C C Updating RAM(ITT0+JADR)=TTIJ IF(LPRED) THEN IRAM(IPRED0+JADR)=IPOS END IF C C END of IF C END of IF RETURN END C C======================================================================= C C C SUBROUTINE SRCREC() C C Auxiliary routine to LOOP updating a source or receiver node. C Called by LOOP if LOOP is called by SOURCE or RECVRS. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ /NAMESR/ /POINTS/ /NODE/ /SRCC/ are C required here: INCLUDE 'net.inc' C net.inc INCLUDE 'netnode.inc' C netnode.inc C C----------------------------------------------------------------------- C EXTERNAL INDX INTEGER INDX C REAL GIANT PARAMETER (GIANT=1.E+10) C INTEGER JPOS,JADR REAL TTJ,TTIJ C C JPOS... Position of the gridpoint 'J' within the model. C JPOS=1,2,3,...,NL1*NL2*NL3. C JADR... Index of the node 'J'. C TTJ... Arrival time at the node 'J'. C TTIJ... Travel time from the node 'I' to the node 'J'. C C....................................................................... C C Index of the node 'J': JPOS=IPOS+KFS4(IFS) JADR=INDX(JPOS) C C Check for the Fresnel volume: IF(JADR.LE.0) RETURN C C Check for a free space: IF(RAM(IP0+JADR).GE.GIANT) RETURN C TTJ=RAM(ITT0+JADR) C C Check for crossing an interface: IF(LICB) THEN IF(IRAM(ICB0+JADR).NE.ICBI) THEN IF(KFS5(IFS).GT.1) RETURN END IF END IF C C*********************************************************************** C Consolidating forward and backward stars (optional): * IF(LNFS) THEN * IF(KFS5(IFS).GT.IRAM(NFS0+JADR)) RETURN * END IF C*********************************************************************** C C Travel time from the node 'I' to the node 'J': IF(LN1) THEN TTIJ=SQRT((FLOAT(KFS2(IFS))*DSX1-DPOS1)**2 * +(FLOAT(KFS2(IFS))*DSX2-DPOS2)**2 * +(FLOAT(KFS3(IFS))*DSX3-DPOS3)**2) * *0.5*(RAM(IP0+JADR)+PI) ELSE TTIJ=SQRT((FLOAT(KFS1(IFS))*DSX1-DPOS1)**2 * +(FLOAT(KFS2(IFS))*DSX2-DPOS2)**2 * +(FLOAT(KFS3(IFS))*DSX3-DPOS3)**2) * *0.5*(RAM(IP0+JADR)+PI) END IF C IF(ISRCI.LE.0) THEN IF(TTI+TTIJ.GE.TTJ) RETURN IF(TTJ.GE.GIANT) THEN C Updating queue by adding 'J' as Q(MAXQ) MAXQ=MAXQ+1 IRAM(IPOSQ0+MAXQ)=IPOS+KFS4(IFS) END IF RAM(ITT0+JADR)=TTI+TTIJ IF(LPRED) THEN IRAM(IPRED0+JADR)=ISRCI END IF C END of IF ELSE IF(TTJ+TTIJ.GE.TTI) RETURN TTI=TTJ+TTIJ IPREDR(ISRCI)=JPOS C END of IF END IF C C END of IF C END of IF RETURN END C C======================================================================= C C C SUBROUTINE READFS() C C Procedure called by SOURCE, to read and construct optimized template C forward stars. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER LU1 PARAMETER(LU1=1) CHARACTER*80 FSTAB C INTEGER MD PARAMETER (MD=292) C MD=292 is good for optimized spherical 3-D forward stars up to C size MFSMAX=27. For minimum values of MD refer to files C net.fs2 and net.fs3. C INTEGER ID(MD),JD(MD),KD(MD) INTEGER JFS,ND,I1,I2,I3,I,L C C ID(II),JD(II),KD(II)... Vector specifying II-th node of the C forward star, in dimensions of a small brick. C JFS... Size of a template forward star currently being read. C ND... Number of nodes stored in the file for a forward star. C I1,I2,I3,I,L... Temporary and loop variables. C C....................................................................... C C Reading and assembling optimized forward star: WRITE(*,'(A)') * '+NET: Reading and assembling the forward star. ' C C Reading nodes of the forward stars: IF(NFSMAX.GT.0) THEN IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN OPEN(LU1,FILE='net.fs2',STATUS='OLD') ELSE OPEN(LU1,FILE='net.fs3',STATUS='OLD') END IF END IF KFS0(0)=0 DO 79 JFS=1,NFSMAX READ(LU1,*) I,ND IF(I.NE.JFS) THEN C NET-40 CALL ERROR * ('NET-40: Template forward star not found in the file') C NET-40: Template forward star not found in the file: C The file 'net.fs3' or 'net.fs2' does not contain a C template forward star of the size NFSMAX. NFSMAX in C main input data file NET, (3), C should be decreased if positive, C adjusted if zero, C or larger forward stars should be added to file C 'net.fs3' or 'net.fs2'. END IF IF(ND.GT.MD) THEN C NET-41 CALL ERROR('NET-41: Too many input nodes of a forward star') C NET-41: Too many input nodes of a forward star: C NFSMAX in main input data file NET, C (3), should be decreased if positive, C adjusted if zero, C or the dimension MD in the subroutine READFS should be C increased. END IF READ(LU1,*) (ID(L),JD(L),KD(L),L=1,ND) C C Assembling whole optimized forward star from the given nodes: KFS0(JFS)=KFS0(JFS-1) DO 69 L=1,ND I1=ID(L) I2=JD(L) I3=KD(L) CALL SETFS(JFS,I1,I2,I3) 69 CONTINUE C 79 CONTINUE CLOSE(LU1) C C Printing the table of the numbers of nodes corresponding to the C template forward stars of sizes 1,2,...,NFSMAX: CALL RSEP3T('FSTAB',FSTAB,' ') IF(FSTAB.NE.' ') THEN OPEN(LU1,FILE=FSTAB) WRITE(LU1,'(A)') ' Size Number Sum of' WRITE(LU1,'(A)') ' (NFS) of nodes numbers' WRITE(LU1,'(3I10)') * (JFS,KFS0(JFS)-KFS0(JFS-1),KFS0(JFS),JFS=1,NFSMAX) CLOSE(LU1) WRITE(*,'(2A)') '+NET: FS table written to file ',FSTAB(1:48) C The central node of a forward star is not stored in the memory C and is not counted. END IF C RETURN END C C======================================================================= C C C SUBROUTINE MAKEFS(NFSMAK) INTEGER NFSMAK C C Procedure called by SOURCE and RECVRS, to construct full spherical C template forward stars. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER NDIM,NH,IH,I1,I2,I3,I3I3,I REAL DH C C....................................................................... C IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN NDIM=2 ELSE NDIM=3 END IF NH=NFSMAK*NFSMAK+NDIM-1 DO 11 I=1,NFSMAK KFS0(I)=0 11 CONTINUE C DH=1./SQRT(12.*FLOAT(NH)) DH=DH/2. DO 69 IH=0,NH DO 58 I1=INT(SQRT(FLOAT(IH)/3.)-DH+1.),INT(SQRT(FLOAT(IH))+DH) I=IH-I1*I1 DO 57 I2=INT(SQRT(FLOAT(I)/2.)-DH+1.), * MIN0(INT(SQRT(FLOAT(I))+DH),I1) I3I3=I-I2*I2 I3=INT(SQRT(FLOAT(I3I3))+0.500) IF(I3*I3.EQ.I3I3.AND.I3.LE.I2) THEN C Assembling the forward star: CALL SETFS(NFSMAK,I3,I2,I1) END IF 57 CONTINUE 58 CONTINUE 69 CONTINUE C DO 71 I=NFSMAK,NFSMAX KFS0(I)=KFS0(NFSMAK) 71 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE SETFS(JFS,I1,I2,I3) INTEGER JFS,I1,I2,I3 C C Subroutine designed to set up all possible mirrorings of an edge C of the forward star. Auxiliary subroutine to READFS and MAKEFS. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER NLEX1,NLEX2,NLEX3,LFS,I REAL DIM C C NLEX1,NLEX2,NLEX3... C LFS... Number of template forward star nodes already stored. C NODES KFS0(JFS),...,LFS are being mirrored with respect to C the central node, while storing the results of mirroring. C I... Loop variable. C DIM... Dimension of the model (2-D or 3-D). C C....................................................................... C LFS=KFS0(JFS) IF(NL1.EQ.1.OR.NL2.EQ.1.OR.NL3.EQ.1) THEN DIM=2. ELSE DIM=3. END IF NLEX1=MIN0(JFS,NL1-1) NLEX2=MIN0(JFS,NL2-1) NLEX3=MIN0(JFS,NL3-1) C IF(I1.LE.NLEX1.AND.I2.LE.NLEX2.AND.I3.LE.NLEX3) THEN CALL STORFS(LFS,I1,I2,I3) END IF IF(I1.NE.I3) THEN IF(I2.LE.NLEX1.AND.I3.LE.NLEX2.AND.I1.LE.NLEX3) THEN CALL STORFS(LFS,I2,I3,I1) END IF IF(I3.LE.NLEX1.AND.I1.LE.NLEX2.AND.I2.LE.NLEX3) THEN CALL STORFS(LFS,I3,I1,I2) END IF IF(I1.NE.I2.AND.I2.NE.I3) THEN IF(I3.LE.NLEX1.AND.I2.LE.NLEX2.AND.I1.LE.NLEX3) THEN CALL STORFS(LFS,I3,I2,I1) END IF IF(I1.LE.NLEX1.AND.I3.LE.NLEX2.AND.I2.LE.NLEX3) THEN CALL STORFS(LFS,I1,I3,I2) END IF IF(I2.LE.NLEX1.AND.I1.LE.NLEX2.AND.I3.LE.NLEX3) THEN CALL STORFS(LFS,I2,I1,I3) END IF END IF END IF DO 61 I=KFS0(JFS)+1,LFS IF(KFS1(I).NE.0) THEN CALL STORFS(LFS,-KFS1(I), KFS2(I), KFS3(I)) END IF 61 CONTINUE DO 62 I=KFS0(JFS)+1,LFS IF(KFS2(I).NE.0) THEN CALL STORFS(LFS, KFS1(I),-KFS2(I), KFS3(I)) END IF 62 CONTINUE DO 63 I=KFS0(JFS)+1,LFS IF(KFS3(I).NE.0) THEN CALL STORFS(LFS, KFS1(I), KFS2(I),-KFS3(I)) END IF 63 CONTINUE DO 64 I=KFS0(JFS)+1,LFS KFS4(I)=KFS1(I)+NL1*(KFS2(I)+NL2*KFS3(I)) KFS5(I)=INT(SQRT(AMAX1(FLOAT(KFS1(I))**2 * +FLOAT(KFS2(I))**2 * +FLOAT(KFS3(I))**2-DIM+1.,1.))+0.999) DFS(I)=SQRT((FLOAT(KFS1(I))*ASX1)**2 * +(FLOAT(KFS2(I))*ASX2)**2 * +(FLOAT(KFS3(I))*ASX3)**2)/2. 64 CONTINUE C KFS0(JFS)=LFS RETURN END C C======================================================================= C C C SUBROUTINE STORFS(LFS,I1,I2,I3) INTEGER LFS,I1,I2,I3 C C Subroutine designed to store one edge of the forward star in the C memory. Auxiliary subroutine to SETFS. C C----------------------------------------------------------------------- C C Common block /FS/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C LFS=LFS+1 IF(LFS.GT.MFS) THEN C NET-42 CALL ERROR('NET-42: Too big forward star') C NET-42: Too big forward star: C NFSMAX in main input data file NET, C (3), should be decreased if positive, adjusted if zero, C or the dimension MFS in the common block /FS/ C in net.inc should be increased. C The minimum value of the dimension MFS may be determined C as follows: C (a) Choose MFS as large as possible. C (b) Add parameter FSTAB='filename' into the input SEP C parameter file and adjust NFSMAX. C (c) Compile and run the program. C (d) Look at file 'filename' just having been C generated. The last integer in the file is the C minimum dimension MFS corresponding to the given C value of NFSMAX. C (e) Update MFS, restore NFSMAX and remove parameter C FSTAB='filename' from the data. END IF KFS1(LFS)=I1 KFS2(LFS)=I2 KFS3(LFS)=I3 RETURN END C C======================================================================= C C C SUBROUTINE ONERAY(LTRACE,LRAYS,LEND,LURAYS,IREFL, * IPREDE,IEND1,IEND2,IEND3,IENDA,TEND,ERR,TERR,ISRC) LOGICAL LTRACE,LRAYS,LEND INTEGER LURAYS,IREFL,IPREDE,IEND1,IEND2,IEND3,IENDA,ISRC REAL TEND,ERR(*),TERR C C Subroutine called by ERRORS and TRACER, to perform backward ray C tracing. C C Input: C LTRACE,LRAYS,LEND,LURAYS,IREFL,IPREDE,IEND1,IEND2,IEND3,IENDA,TEND C IPREDE... Must be positive. C C Output: C TERR C C----------------------------------------------------------------------- C C Common blocks /GRID/ /NAMESR/ /POINTS/ /FILES/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C EXTERNAL INDX INTEGER INDX C INTEGER LU0,LU1,LU PARAMETER(LU0=10) PARAMETER(LU1=1) C INTEGER MTMP PARAMETER (MTMP=100) INTEGER KTMP(MTMP),NTMP,ITMP REAL TMPMIN(MTMP),TMPMAX(MTMP) C INTEGER IPREDI,IPNEW,NF,N123,L123,IPOS1,IPOS2,IPOS3,IPOS,IADR INTEGER IOLD1,IOLD2,IOLD3,ICBOLD,NFOLD,K REAL DPOS1,DPOS2,DPOS3,X1,X2,X3,T,TOLD,POLD,PAUX,ERRMIN,ERRMAX C C IPOS1,IPOS2,IPOS3,IPOS,IADR... Position and index (address) of the C nearest network node. C DPOS1,DPOS2,DPOS3... Offset from the nearest network node. C IPREDI..Predecessor of the current node. C IOLD1,IOLD2,IOLD3,ICBOLD,NFOLD,TOLD,POLD... Parameters of the last C network node along the ray. C X1,X2,X3,T,PAUX... Temporary storage locations. C C....................................................................... C IPREDI=IPREDE NTMP=0 IF(.NOT.LTRACE.OR.LEND) THEN IF(LICB) THEN ICBOLD=IRAM(ICB0+IENDA) ELSE ICBOLD=1 END IF IOLD1=IEND1 IOLD2=IEND2 IOLD3=IEND3 POLD=RAM(IP0+IENDA) C POLD is a constant approximation, not interpolation. TOLD=TEND NFOLD=IREFL ERRMIN=0. ERRMAX=0. END IF IF(IPREDI.GE.0) THEN DO 68 NF=IREFL,0,-1 IF(IPREDI.EQ.0) THEN C Endpoint of the ray at the reflector: C This is not possible at the receiver. C Thus, this is possible only if called by subroutine ERRORS. C Then IOLD1,IOLD2,IOLD3 are defined. IPREDI=1+IOLD1+(NL1*IOLD2+NL2*IOLD3) GO TO 66 ENDIF IF(NF.LT.IREFL) THEN LU=LU0+NF IF(L4.NE.0) THEN C Reading the index array: N123=N1*N2*N3 L123=L1*L2*L3 DO 11 K=1,N123 IRAM(IND0+K)=K 11 CONTINUE CALL RARRAI(LU1,FIND(NF),'FORMATTED', * .FALSE.,0,N123,IRAM(IND0+1)) DO 12 K=1,N123 IRAM(IND0+K)=(IRAM(IND0+K)-1)*L123+1 12 CONTINUE END IF END IF 50 CONTINUE IADR=INDX(IPREDI) IF(NF.LT.IREFL) THEN READ(LU,REC=IADR) T,IPNEW ELSE T=RAM(ITT0+IADR) IPNEW=IRAM(IPRED0+IADR) END IF IPOS1=IPREDI-1 IPOS2=IPOS1/NL1 IPOS3=IPOS2/NL2 IPOS1=IPOS1-IPOS2*NL1 IPOS2=IPOS2-IPOS3*NL2 IF(LTRACE.AND.LRAYS) THEN IF(LN1) THEN X1=X1MIN+(FLOAT(IPOS2)+0.5)*DSX1 ELSE X1=X1MIN+(FLOAT(IPOS1)+0.5)*DSX1 END IF X2=X2MIN+(FLOAT(IPOS2)+0.5)*DSX2 X3=X3MIN+(FLOAT(IPOS3)+0.5)*DSX3 WRITE(LURAYS,'(4F12.6,A)') X1,X2,X3,T,' /' END IF IF(.NOT.LTRACE.OR.LEND) THEN C Travel time error estimate: CALL SETERR(NF,IPREDI,IPOS1,IPOS2,IPOS3,IADR,T,NFOLD, * IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX) IF(.NOT.LTRACE) THEN IF(NF.EQ.IREFL) THEN IF(ERR(IADR).GE.0.) THEN ERRMIN=ERRMIN-ERR(IADR) ERRMAX=ERRMAX+ERR(IADR) GO TO 70 ELSE IF(NTMP.LT.MTMP) THEN NTMP=NTMP+1 KTMP(NTMP)=IADR TMPMIN(NTMP)=ERRMIN TMPMAX(NTMP)=ERRMAX END IF END IF END IF END IF END IF IF(IPNEW.EQ.0) THEN C Reflector: C The same point at the reflector is considered twice, C i.e. the ray segment of zero length is situated at the C reflector. GO TO 66 ELSE IF(IPNEW.LT.0) THEN C Source point: IPREDI=IPNEW GO TO 69 ELSE IPREDI=IPNEW ENDIF GO TO 50 66 CONTINUE C Node IPREDI is situated at a reflector. * AUX=POLD*AMAX1(ASX1,ASX2,ASX3)/2. * ERRMIN=ERRMIN-AUX * ERRMAX=ERRMIN+AUX 68 CONTINUE 69 CONTINUE END IF C Node IPREDI is a source point. ISRC=-IPREDI C C Source point: IF(.NOT.LTRACE.OR.LEND) THEN C Travel time error estimate: CALL SLOW(X1S(ISRC),X2S(ISRC),X3S(ISRC), * DPOS1,DPOS2,DPOS3,IPOS1,IPOS2,IPOS3,IPOS,IADR,PAUX) CALL SETERR(IREFL,IPOS,IPOS1,IPOS2,IPOS3,IADR,TTS(ISRC),NFOLD, * IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX) ERRMIN=ERRMIN-TTSERR(ISRC) ERRMAX=ERRMAX+TTSERR(ISRC) END IF C 70 CONTINUE C IF(.NOT.LTRACE) THEN C Rays are not traced, the end node is a gridpoint, C thus IENDA is defined. ERR(IENDA)=AMAX1(ERRMAX,-ERRMIN) DO 71 ITMP=1,NTMP ERR(KTMP(ITMP))= * AMAX1(ERRMAX-TMPMAX(ITMP),-ERRMIN+TMPMIN(ITMP)) 71 CONTINUE ELSE IF(LEND) THEN TERR=AMAX1(ERRMAX,-ERRMIN) END IF C RETURN END C C======================================================================= C C C SUBROUTINE SETERR(NF,IPOS,IPOS1,IPOS2,IPOS3,IADR,T,NFOLD, * IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX) INTEGER NF,IPOS,IPOS1,IPOS2,IPOS3,IADR INTEGER NFOLD,IOLD1,IOLD2,IOLD3,ICBOLD REAL T,POLD,TOLD,ERRMIN,ERRMAX C C Subroutine called by ONERAY, to estimate and accumulate the C arrival-time errors during backward ray tracing. C C Attention: C In this version, in a case of reflections, the slowness field is C correct only if it is the same for the direct wave and for all the C reflected waves. When the rays are traced, only time fields and C predecessors are alternated in the memory, unlike slowness fields. C C Input: C NF,IPOS,IPOS1,IPOS2,IPOS3,IADR,T, C NFOLD,IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD,ERRMIN,ERRMAX C C Output: C NFOLD,IOLD1,IOLD2,IOLD3,ICBOLD,POLD,TOLD... Values at the old node C are replaced by the values at the new node. C ERRMIN,ERRMAX... Input values increased by the increment of the C error bounds between the old node and the new node. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER ICBNEW,NFS1,IERR REAL UIJ(6),C1,C2MIN,C2MAX,ERR3,PMEAN REAL AUX1,AUX2,AUX3 C C....................................................................... C CALL OPTMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR,UIJ) IF(LNFS) THEN IF(NF.EQ.NREFL) THEN NFS1=IRAM(NFS0+IADR) ELSE CALL OPTNFS(IPOS,IPOS1,IPOS2,IPOS3,IADR,NFS1,C1,C2MIN,C2MAX) END IF ELSE NFS1=NFSMAX END IF C C Calculating the coefficient of the network-geometry error: IF(NL1.EQ.1) THEN C1=((ASX2**2+ASX3**2)/AMIN1(ASX2,ASX3)**2-1.)/8. ELSE IF(NL2.EQ.1) THEN C1=((ASX1**2+ASX3**2)/AMIN1(ASX1,ASX3)**2-1.)/8. ELSE IF(NL3.EQ.1) THEN C1=((ASX1**2+ASX2**2)/AMIN1(ASX1,ASX2)**2-1.)/8. ELSE C1=((ASX1**2+ASX2**2+ASX3**2)/AMIN1(ASX1,ASX2,ASX3)**2-1.)/8. END IF C C Calculating the relative heterogeneity error: AUX1=FLOAT(IPOS1-IOLD1) AUX2=FLOAT(IPOS2-IOLD2) AUX3=FLOAT(IPOS3-IOLD3) AUX1=AUX1*(UIJ(1)*AUX1+2.*(UIJ(2)*AUX2+UIJ(4)*AUX3)) * +AUX2*(UIJ(3)*AUX2+2.*UIJ(5)*AUX3)+AUX3*UIJ(6)*AUX3 AUX1=AUX1/(12.*RAM(IP0+IADR)) C C Geometry and heterogeneity increments of the error bounds: ERRMIN=ERRMIN+(TOLD-T)* AUX1 ERRMAX=ERRMAX+(TOLD-T)*(AUX1+C1/FLOAT(NFS1*NFS1+1)) C C Error due to structural interfaces: IF(LICB) THEN ICBNEW=IRAM(ICB0+IADR) ELSE ICBNEW=1 END IF IF(ICBNEW.NE.ICBOLD) THEN C For ABS(DSX1)=ABS(DSX2)=ABS(DSX3)=DSX: C D(ERRMIN)=-DSX* D(P)*SQRT(N)/2 C D(ERRMAX)=DSX*MIN(P, C (D(P)*SQRT(2)-SQRT(MIN(P)**2+((SQRT(N)-SQRT(2))*P)**2)) AUX1=AMAX1(ASX1,ASX2,ASX3) AUX3=AMIN1(ASX1,ASX2,ASX3) AUX2=ASX1+ASX2+ASX3-AUX1-AUX3 PMEAN=(RAM(IP0+IADR)+POLD)/2. ERRMIN=ERRMIN-SQRT(AUX1*AUX1+AUX2*AUX2+AUX3*AUX3) * *ABS(RAM(IP0+IADR)-POLD)/2. ERR3=AMIN1(RAM(IP0+IADR),POLD) C In the worst case, interface is perpendicular to the greatest C grid step AUX1. IF(AUX3.LE.0.) THEN C 2-D: C Now, ERR3 is the time derivative in the direction AUX2. ERR3=PMEAN*SQRT(AUX1*AUX1+AUX2*AUX2)-ERR3*AUX2 ELSE C 3-D: ERR3=ERR3*ERR3-(PMEAN*( SQRT(AUX1*AUX1+AUX2*AUX2+AUX3*AUX3) * -SQRT(AUX1*AUX1+AUX3*AUX3) )/AUX2)**2 IF(ERR3.GT.0.) THEN ERR3=SQRT(ERR3) ELSE ERR3=0. END IF C Now, ERR3 is the time derivative in the direction AUX3. ERR3=PMEAN*SQRT(AUX1*AUX1+AUX3*AUX3)-ERR3*AUX3 END IF ERR3=AMIN1(ERR3,PMEAN*AUX1) ERRMAX=ERRMAX+ERR3 END IF C C Error due to reflections: IF(NFOLD.NE.NF) THEN AUX1=(RAM(IP0+IADR)+POLD)*AMAX1(ASX1,ASX2,ASX3)/2. ERRMIN=ERRMIN-AUX1 ERRMAX=ERRMIN+AUX1 NFOLD=NF END IF C ICBOLD=ICBNEW IOLD1=IPOS1 IOLD2=IPOS2 IOLD3=IPOS3 POLD=RAM(IP0+IADR) TOLD=T C RETURN END C C======================================================================= C C C SUBROUTINE OPTNFS * (IPOS,IPOS1,IPOS2,IPOS3,IADR,NFSOPT,C1,C2MIN,C2MAX) INTEGER IPOS,IPOS1,IPOS2,IPOS3,IADR,NFSOPT REAL C1,C2MIN,C2MAX C C Subroutine called by SOURCE and SETERR, designed to estimate the C optimum size of a forward star at the given node. C C Input: C IPOS,IPOS1,IPOS2,IPOS3... Indices of the given gridpoint, C describing the position of the network node situated at C the gridpoint. C IADR... Index of the network node. C C Output: C NFSOPT..Optimum size of the forward star at the given node. C Default of NFSOPT=1 is set if the optimum size cannot be C determined. C C1,C2MIN,C2MAX... Coefficients of the travel time error estimate. C C----------------------------------------------------------------------- C C Common blocks /GRID/ /FS/ are required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C REAL GIANT PARAMETER (GIANT=1.E+10) C INTEGER IERR REAL C1SAVE,C2,UIJ(6),AUX0 SAVE C1SAVE DATA C1SAVE/0./ C C C1SAVE..Coefficient C1. C C2... Maximum absolute value of coefficients C2MIN, C2MAX. C UIJ... Symmetric matrix describing the level of heterogeneity. C AUX0... Dummy storage location. C C....................................................................... C C Calculating the coefficient of the network-geometry error: IF(C1SAVE.EQ.0.) THEN IF(NL1.EQ.1) THEN C1SAVE=((ASX2**2+ASX3**2)/AMIN1(ASX2,ASX3)**2-1.)/8. ELSE IF(NL2.EQ.1) THEN C1SAVE=((ASX1**2+ASX3**2)/AMIN1(ASX1,ASX3)**2-1.)/8. ELSE IF(NL3.EQ.1) THEN C1SAVE=((ASX1**2+ASX2**2)/AMIN1(ASX1,ASX2)**2-1.)/8. ELSE C1SAVE= * ((ASX1**2+ASX2**2+ASX3**2)/AMIN1(ASX1,ASX2,ASX3)**2-1.)/8. END IF END IF C C Coefficient of the network-geometry error (rel.error=C1/NFS**2): C1=C1SAVE C C Check for a free space: IF(RAM(IP0+IADR).GE.GIANT) THEN NFSOPT=0 C2MIN=0. C2MAX=0. RETURN END IF C CALL OPTMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR,UIJ) C IERR=0,1,...,16. Free space: IERR=16. IF(IERR.GE.16) THEN C Optimum size of a forward star cannot be determined: NFSOPT=1 RETURN END IF C C Coefficient of the heterogeneity error (rel.error=C2*NFS**2): CALL EIGEN(UIJ,AUX0,3,1) C2=12.*RAM(IP0+IADR) C2MAX=AMAX1(UIJ(1),0.)/C2 C2MIN=AMIN1(UIJ(6),0.)/C2 C2=AMAX1(C2MAX,-C2MIN-C2MAX) C C Maximum size of the forward star: IF(NFSMAX.GT.0) THEN NFSOPT=NFSMAX ELSE NFSOPT=999999 END IF C C Optimum size of the forward star: IF(C2.GT.0.) THEN NFSOPT=MAX0(1,MIN0(INT(SQRT(SQRT(C1/C2))+0.5),NFSOPT)) END IF C RETURN END C C======================================================================= C C C SUBROUTINE OPTMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR,UIJ) INTEGER IPOS,IPOS1,IPOS2,IPOS3,IADR,IERR REAL UIJ(6) C C Auxiliary subroutine to SETERR and OPTNFS, designed to calculate the C symmetric matrix describing the level of heterogeneity at the given C node. The matrix is composed of the first and second slowness C derivatives. C C Input: C IPOS,IPOS1,IPOS2,IPOS3... Indices of the given gridpoint, C describing the position of the network node situated at C the gridpoint. C IADR... Index of the network node. C C Output: C IERR... IERR=0 if the matrix is determined. C UIJ... Symmetric matrix describing the level of heterogeneity. C C----------------------------------------------------------------------- C C Common block /GRID/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C EXTERNAL INDX INTEGER INDX INTEGER ISHIFT(3),ICBOPT,I1,I2,I3,IP,IP1,IP2,IP3,IAUX INTEGER IP1OLD,IP2OLD,IP3OLD,IEROLD DATA ISHIFT/0,1,-1/ C C....................................................................... C ICBOPT=-999 CALL TRYMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,ICBOPT,IERR,UIJ) C ICBOPT unchanged: node not situated in the computational volume. C IERR.GT.0: matrix Uij has not been determined. C C Check whether matrix Uij is determined: IF(ICBOPT.NE.-999.AND.IERR.GT.0) THEN IEROLD=IERR IP1OLD=IPOS1 IP2OLD=IPOS2 IP3OLD=IPOS3 C C Looking for matrix Uij around: DO 13 I3=1,3 IP3=IPOS3+ISHIFT(I3) C Check for the model volume: IF(0.LE.IP3.AND.IP3.LT.NL3) THEN DO 12 I2=1,3 IP2=IPOS2+ISHIFT(I2) C Check for the model volume: IF(0.LE.IP2.AND.IP2.LT.NL2) THEN DO 11 I1=1,3 IF(I1.NE.1.OR.I2.NE.1.OR.I3.NE.1) THEN IP1=IPOS1+ISHIFT(I1) C Check for the model volume: IF(0.LE.IP1.AND.IP1.LT.NL1) THEN IP=1+IP1+NL1*(IP2+NL2*IP3) IAUX=INDX(IP) CALL TRYMAT(IP,IP1,IP2,IP3,IAUX,ICBOPT,IERR,UIJ) IF(IERR.EQ.0) THEN C Matrix Uij is determined RETURN ELSE IF(IERR.LT.IEROLD) THEN IEROLD=IERR IP1OLD=IP1 IP2OLD=IP2 IP3OLD=IP3 END IF END IF END IF 11 CONTINUE END IF 12 CONTINUE END IF 13 CONTINUE C IP=1+IP1OLD+NL1*(IP2OLD+NL2*IP3OLD) IAUX=INDX(IP) CALL TRYMAT(IP,IP1OLD,IP2OLD,IP3OLD,IAUX,ICBOPT,IERR,UIJ) END IF RETURN END C C======================================================================= C C C SUBROUTINE TRYMAT(IPOS,IPOS1,IPOS2,IPOS3,IADR,ICBOPT,IERR,UIJ) INTEGER IPOS,IPOS1,IPOS2,IPOS3,IADR,ICBOPT,IERR REAL UIJ(6) C C Auxiliary subroutine to OPTMAT, to try the calculation of the matrix C describing the level of heterogeneity at the given node. C The matrix is composed of the first and second slowness derivatives. C C Input: C IPOS,IPOS1,IPOS2,IPOS3... Indices of the given gridpoint, C describing the position of the network node situated at C the gridpoint. C IADR... Index of the network node. C C Output: C ICBOPT.. C IERR... IERR=0 if the matrix is determined. C UIJ... Symmetric matrix describing the level of heterogeneity. C C----------------------------------------------------------------------- C C Common block /GRID/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C REAL GIANT PARAMETER (GIANT=1.E+10) C EXTERNAL INDX INTEGER INDX INTEGER IA0,IA2,IERR12,IERR13,IERR23 REAL U1,U2,U3,U11,U12,U22,U13,U23,U33,AUX0,AUX1,AUX2,AUX3 C C IERR12,IERR13,IERR23... Switches for evaluation of mixed second C derivatives. C U1,U2,U3,U11,U12,U22,U13,U23,U33... Slowness derivatives. C AUX0,AUX1,AUX2,AUX3... Temporary storage locations. C C....................................................................... C C IERR=0,1,...,16. IERR=16 if slowness is not defined at the point. C Otherwise, slowness derivatives cannot be determined along IERR/4 C axes. In addition, MOD(IERR,4) mixed derivatives cannot be C determined. Initially: IERR=16 IERR12=2 IERR13=2 IERR23=2 C C 2-D model and default in 3-D: U1 =0. U2 =0. U3 =0. U11=0. U12=0. U22=0. U13=0. U23=0. U33=0. C C Check for the computational volume at the given node: IF(IADR.EQ.0) THEN RETURN END IF C C Determination of the geological block at the given node: IF(ICBOPT.EQ.-999) THEN IF(LICB) THEN ICBOPT=IRAM(ICB0+IADR) ELSE ICBOPT=1 END IF END IF C C Check for the geological block at the given node: IF(LICB) THEN IF(IRAM(ICB0+IADR).NE.ICBOPT) THEN RETURN END IF END IF C Check for a free space: AUX1=RAM(IP0+IADR) IF(AUX1.GE.GIANT) THEN RETURN END IF IERR=IERR-4 C C First and second partial slowness derivatives along the X1 axis: IF(NL1.GT.1) THEN C 3-D model: C Check for the model volume: IF(IPOS1.LT.1.OR.NL1-2.LT.IPOS1) THEN GO TO 20 END IF C Check for the computational volume: IA0=INDX(IPOS-1) IA2=INDX(IPOS+1) IF(IA0.EQ.0) THEN GO TO 20 END IF IF(IA2.EQ.0) THEN GO TO 20 END IF C Check for the geological block: IF(LICB) THEN IF(IRAM(ICB0+IA0).NE.ICBOPT) THEN GO TO 20 END IF IF(IRAM(ICB0+IA2).NE.ICBOPT) THEN GO TO 20 END IF END IF C Check for a free space: AUX0=RAM(IP0+IA0) AUX2=RAM(IP0+IA2) IF(AUX0.GE.GIANT) THEN GO TO 20 END IF IF(AUX2.GE.GIANT) THEN GO TO 20 END IF C Partial slowness derivatives scaled by the grid step: U1 =(AUX2-AUX0)/2. U11=AUX2-2.*AUX1+AUX0 IERR12=IERR12-1 IERR13=IERR13-1 END IF IERR=IERR-4 C C First and second partial slowness derivatives along the X2 axis: 20 CONTINUE IF(NL2.GT.1) THEN C 3-D model: C Check for the model volume: IF(IPOS2.LT.1.OR.NL2-2.LT.IPOS2) THEN GO TO 30 END IF C Check for the computational volume: IA0=INDX(IPOS-NL1) IA2=INDX(IPOS+NL1) IF(IA0.EQ.0) THEN GO TO 30 END IF IF(IA2.EQ.0) THEN GO TO 30 END IF C Check for the geological block: IF(LICB) THEN IF(IRAM(ICB0+IA0).NE.ICBOPT) THEN GO TO 30 END IF IF(IRAM(ICB0+IA2).NE.ICBOPT) THEN GO TO 30 END IF END IF C Check for a free space: AUX0=RAM(IP0+IA0) AUX2=RAM(IP0+IA2) IF(AUX0.GE.GIANT) THEN GO TO 30 END IF IF(AUX2.GE.GIANT) THEN GO TO 30 END IF C Partial slowness derivatives scaled by the grid step: U2 =(AUX2-AUX0)/2. U22=AUX2-2.*AUX1+AUX0 IERR12=IERR12-1 IERR23=IERR23-1 END IF IERR=IERR-4 C C First and second partial slowness derivatives along the X3 axis: 30 CONTINUE IF(NL3.GT.1) THEN C 3-D model: C Check for the model volume: IF(IPOS3.LT.1.OR.NL3-2.LT.IPOS3) THEN GO TO 40 END IF C Check for the computational volume: IA0=INDX(IPOS-NL1*NL2) IA2=INDX(IPOS+NL1*NL2) IF(IA0.EQ.0) THEN GO TO 40 END IF IF(IA2.EQ.0) THEN GO TO 40 END IF C Check for the geological block: IF(LICB) THEN IF(IRAM(ICB0+IA0).NE.ICBOPT) THEN GO TO 40 END IF IF(IRAM(ICB0+IA2).NE.ICBOPT) THEN GO TO 40 END IF END IF C Check for a free space: AUX0=RAM(IP0+IA0) AUX2=RAM(IP0+IA2) IF(AUX0.GE.GIANT) THEN GO TO 40 END IF IF(AUX2.GE.GIANT) THEN GO TO 40 END IF C Partial slowness derivatives scaled by the grid step: U3 =(AUX2-AUX0)/2. U33=AUX2-2.*AUX1+AUX0 IERR13=IERR13-1 IERR23=IERR23-1 END IF IERR=IERR-4 C C Mixed second partial slowness derivatives: 40 CONTINUE IF(NL1.GT.1.AND.NL2.GT.1.AND.IERR12.EQ.0) THEN CALL MIXDER(IPOS, 1,NL1 ,ICBOPT,U12,IERR) END IF IF(NL1.GT.1.AND.NL3.GT.1.AND.IERR13.EQ.0) THEN CALL MIXDER(IPOS, 1,NL1*NL2,ICBOPT,U13,IERR) END IF IF(NL2.GT.1.AND.NL3.GT.1.AND.IERR23.EQ.0) THEN CALL MIXDER(IPOS,NL1,NL1*NL2,ICBOPT,U23,IERR) END IF C C Symmetric matrix describing the level of heterogeneity: AUX0=0.5/RAM(IP0+IADR) AUX1=U1*AUX0 AUX2=U2*AUX0 AUX3=U3*AUX0 AUX0=AUX1*U1+AUX2*U2+AUX3*U3 UIJ(1)=U11-AUX1*U1+AUX0 UIJ(2)=U12-AUX1*U2 UIJ(3)=U22-AUX2*U2+AUX0 UIJ(4)=U13-AUX1*U3 UIJ(5)=U23-AUX2*U3 UIJ(6)=U33-AUX3*U3+AUX0 C RETURN END C C======================================================================= C C C SUBROUTINE MIXDER(IPOS0,IPOS1,IPOS2,ICBREF,DERMIX,IERR) INTEGER IPOS0,IPOS1,IPOS2,ICBREF,IERR REAL DERMIX C C Auxiliary subroutine to TRYMAT, to calculate mixed second partial C slowness derivatives. C C----------------------------------------------------------------------- C C Common block /GRID/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C REAL GIANT PARAMETER (GIANT=1.E+10) C EXTERNAL INDX INTEGER INDX INTEGER IA00,IA20,IA02,IA22 C C....................................................................... C C Addresses of the corner gridpoints: IA00=INDX(IPOS0-IPOS1-IPOS2) IA20=INDX(IPOS0+IPOS1-IPOS2) IA02=INDX(IPOS0-IPOS1+IPOS2) IA22=INDX(IPOS0+IPOS1+IPOS2) C C Check for the same geological block and for a free space: IF(IA00.GT.0) THEN IF(LICB) THEN IF(IRAM(ICB0+IA00).NE.ICBREF) THEN IA00=0 GO TO 11 END IF END IF IF(RAM(IP0+IA00).GE.GIANT) THEN IA00=0 END IF END IF 11 CONTINUE IF(IA20.GT.0) THEN IF(LICB) THEN IF(IRAM(ICB0+IA20).NE.ICBREF) THEN IA20=0 GO TO 12 END IF END IF IF(RAM(IP0+IA20).GE.GIANT) THEN IA20=0 END IF END IF 12 CONTINUE IF(IA02.GT.0) THEN IF(LICB) THEN IF(IRAM(ICB0+IA02).NE.ICBREF) THEN IA02=0 GO TO 13 END IF END IF IF(RAM(IP0+IA02).GE.GIANT) THEN IA02=0 END IF END IF 13 CONTINUE IF(IA22.GT.0) THEN IF(LICB) THEN IF(IRAM(ICB0+IA22).NE.ICBREF) THEN IA22=0 GO TO 14 END IF END IF IF(RAM(IP0+IA22).GE.GIANT) THEN IA22=0 END IF END IF 14 CONTINUE C C Calculating mixed second partial slowness derivatives: IF(IA22.GT.0) THEN IF(IA02.GT.0) THEN IF(IA20.GT.0) THEN IF(IA00.GT.0) THEN DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/4. ELSE IA00=INDX(IPOS0-IPOS1) IA20=INDX(IPOS0+IPOS1) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. END IF ELSE IF(IA00.GT.0) THEN IA20=INDX(IPOS0-IPOS2) IA22=INDX(IPOS0+IPOS2) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. ELSE IA00=INDX(IPOS0-IPOS1) IA20=INDX(IPOS0+IPOS1) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. END IF END IF ELSE IF(IA20.GT.0) THEN IF(IA00.GT.0) THEN IA02=INDX(IPOS0-IPOS1) IA22=INDX(IPOS0+IPOS1) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. ELSE IA00=INDX(IPOS0-IPOS2) IA02=INDX(IPOS0+IPOS2) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. END IF ELSE IF(IA00.GT.0) THEN IA02=INDX(IPOS0-IPOS1) IA20=INDX(IPOS0-IPOS2) IA22=INDX(IPOS0) DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22) ELSE IA00=INDX(IPOS0) IA20=INDX(IPOS0+IPOS1) IA02=INDX(IPOS0+IPOS2) DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22) END IF END IF END IF ELSE IF(IA02.GT.0) THEN IF(IA20.GT.0) THEN IF(IA00.GT.0) THEN IA20=INDX(IPOS0-IPOS2) IA22=INDX(IPOS0+IPOS2) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. ELSE IA00=INDX(IPOS0-IPOS1) IA20=INDX(IPOS0) IA22=INDX(IPOS0+IPOS2) DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22) END IF ELSE IF(IA00.GT.0) THEN IA20=INDX(IPOS0-IPOS2) IA22=INDX(IPOS0+IPOS2) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. ELSE IA00=INDX(IPOS0-IPOS1) IA20=INDX(IPOS0) IA22=INDX(IPOS0+IPOS2) DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22) END IF END IF ELSE IF(IA20.GT.0) THEN IF(IA00.GT.0) THEN IA02=INDX(IPOS0-IPOS1) IA22=INDX(IPOS0+IPOS1) DERMIX=(RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22))/2. ELSE IA00=INDX(IPOS0-IPOS2) IA02=INDX(IPOS0) IA22=INDX(IPOS0+IPOS1) DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22) END IF ELSE IF(IA00.GT.0) THEN IA02=INDX(IPOS0-IPOS1) IA20=INDX(IPOS0-IPOS2) IA22=INDX(IPOS0) DERMIX= RAM(IP0+IA00)-RAM(IP0+IA20) * -RAM(IP0+IA02)+RAM(IP0+IA22) ELSE IERR=IERR+1 END IF END IF END IF END IF C RETURN END C C======================================================================= C C C SUBROUTINE SLOW(X1,X2,X3,DPOS1,DPOS2,DPOS3, * IPOS1,IPOS2,IPOS3,IPOS,IADR,PS) REAL X1,X2,X3,DPOS1,DPOS2,DPOS3,PS INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR C C Subroutine called by SOURCE, RECVRS, TRACER, and ONERAY, to find C a close gridpoint within the computational volume, and to interpolate C the slowness within the same geological block. C If the distance of the given point from the model volume is greater C than a grid step, an error message is generated. C If the point is outside the computational volume or in a free space, C another error message is generated. C C Input: C X1,X2,X3... Coordinates of the point. C C Output: C DPOS1,DPOS2,DPOS3... Vectorial distance from the nearest C gridpoint situated in the same geological block. C IPOS1,IPOS2,IPOS3... Positional indices of the nearest gridpoint. C IPOS... Index of the nearest gridpoint in the index array. C IADR... Storage address of the nearest gridpoint if it is located C in the computational volume, otherwise zero. C PS... Interpolated slowness. C C----------------------------------------------------------------------- C C Common block /GRID/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C EXTERNAL INDX INTEGER INDX REAL GIANT PARAMETER (GIANT=1.E+10) C INTEGER I1,I2,I3,I1D,I2D,I3D,IAUX,ICBREF REAL DS,AUX C C I1,I2,I3... Loop variables - position of a gridpoint. C I1D,I2D,I3D... Loop increments. C IAUX... Index of the node at a gridpoint. C DS... Distance of the point from a gridpoint. C AUX... Temporary variable. C ICBREF... Index of the geological block. C C....................................................................... C C Determination of the nearest gridpoint: CALL POS(X1,X2,X3,DPOS1,DPOS2,DPOS3,IPOS1,IPOS2,IPOS3,IPOS,IADR) I1D=INT(SIGN(1.5,DSX1*DPOS1)) I2D=INT(SIGN(1.5,DSX2*DPOS2)) I3D=INT(SIGN(1.5,DSX3*DPOS3)) C C Slowness in the given point is determined by an interpolation C from neighbouring grid-point slownesses PS=0.0 DS=0.0 DO 23 I3=IPOS3,IPOS3+I3D,I3D DO 22 I2=IPOS2,IPOS2+I2D,I2D DO 21 I1=IPOS1,IPOS1+I1D,I1D IF(0.LE.I1.AND.I1.LE.NL1-1.AND. * 0.LE.I2.AND.I2.LE.NL2-1.AND. * 0.LE.I3.AND.I3.LE.NL3-1) THEN IAUX=1+I1+(I2+I3*NL2)*NL1 IAUX=INDX(IAUX) IF(IAUX.GT.0) THEN C C Test for nodes in a free space: IF(RAM(IP0+IAUX).GE.GIANT) THEN GO TO 20 ENDIF C C Check for the nearest gridpoint in the indexed volume IF(DS.EQ.0.) THEN C C Redefinition of the nearest gridpoint IF(IADR.NE.IAUX) THEN IF(LN1) THEN DPOS1=DPOS1+DSX1*FLOAT(IPOS2-I2) ELSE DPOS1=DPOS1+DSX1*FLOAT(IPOS1-I1) END IF DPOS2=DPOS2+DSX2*FLOAT(IPOS2-I2) DPOS3=DPOS3+DSX3*FLOAT(IPOS3-I3) IPOS1=I1 IPOS2=I2 IPOS3=I3 IPOS=1+I1+(I2+I3*NL2)*NL1 IADR=IAUX END IF C C Setting reference geological block IF(LICB) THEN ICBREF=IRAM(ICB0+IADR) END IF ELSE C C Check for the (geological) block: IF(LICB) THEN IF(IRAM(ICB0+IAUX).NE.ICBREF) THEN GO TO 20 ENDIF ENDIF ENDIF C C Interpolation: IF(LN1) THEN AUX=SQRT((FLOAT(I2-IPOS2)*DSX1-DPOS1)**2+ * (FLOAT(I2-IPOS2)*DSX2-DPOS2)**2+ * (FLOAT(I3-IPOS3)*DSX3-DPOS3)**2) ELSE AUX=SQRT((FLOAT(I1-IPOS1)*DSX1-DPOS1)**2+ * (FLOAT(I2-IPOS2)*DSX2-DPOS2)**2+ * (FLOAT(I3-IPOS3)*DSX3-DPOS3)**2) END IF IF(AUX.EQ.0.) THEN PS=RAM(IP0+IAUX) RETURN END IF PS=PS+RAM(IP0+IAUX)/AUX DS=DS+1.0/AUX C ENDIF ENDIF 20 CONTINUE 21 CONTINUE 22 CONTINUE 23 CONTINUE IF(DS.EQ.0.) THEN C NET-45 CALL ERROR('NET-45: Source or receiver point in a free space') C NET-45: Source or receiver point in a free space: C Source or receiver point is situated outside the indexed C computational volume or in a free space (zero velocity). ENDIF PS=PS/DS C RETURN END C C======================================================================= C C C SUBROUTINE POS(X1,X2,X3,DPOS1,DPOS2,DPOS3, * IPOS1,IPOS2,IPOS3,IPOS,IADR) REAL X1,X2,X3,DPOS1,DPOS2,DPOS3 INTEGER IPOS1,IPOS2,IPOS3,IPOS,IADR C C Auxiliary subroutine to SLOW, to find the nearest gridpoint. C This subroutine checks the position of the given point with respect to C the boundaries of the model volume. If the distance from the model C volume is greater than a grid step, the error message is generated. C C Input: C X1,X2,X3... Coordinates of the point. C C Output: C DPOS1,DPOS2,DPOS3... Vectorial distance from the nearest C gridpoint. C IPOS1,IPOS2,IPOS3... Positional indices of the nearest gridpoint. C IPOS... Index of the nearest gridpoint in the index array. C IADR... Storage address of the nearest gridpoint if it is located C in the computational volume, otherwise zero. C C----------------------------------------------------------------------- C C Common block /GRID/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C EXTERNAL INDX INTEGER INDX C C....................................................................... C IF(LN1) THEN IPOS1=0 IPOS2=INT( ((X1-X1MIN)*DSX1+(X2-X2MIN)*DSX2)/(DSX1**2+DSX2**2) ) IF(IPOS2.LT.0.OR.NL2.LT.IPOS2) THEN C NET-46 CALL ERROR * ('NET-46: Source or receiver point outside the model') C NET-46: Source or receiver point outside the model: C Source or receiver point is situated outside the model C volume. ELSE IF(IPOS2.GE.NL2) THEN IPOS2=NL2-1 END IF DPOS1=X1-X1MIN-(FLOAT(IPOS2)+0.5)*DSX1 DPOS2=X2-X2MIN-(FLOAT(IPOS2)+0.5)*DSX2 ELSE IF(NL1.EQ.1) THEN IPOS1=0 DPOS1=0. ELSE IPOS1=INT((X1-X1MIN)/DSX1) IF(IPOS1.LT.0.OR.NL1.LT.IPOS1) THEN C NET-47 CALL ERROR * ('NET-47: Source or receiver point outside the model') C NET-47: Source or receiver point outside the model: C Source or receiver point is situated outside the C model volume. ELSE IF(IPOS1.GE.NL1) THEN IPOS1=NL1-1 END IF DPOS1=X1-X1MIN-(FLOAT(IPOS1)+0.5)*DSX1 END IF IF(NL2.EQ.1) THEN IPOS2=0 DPOS2=0. ELSE IPOS2=INT((X2-X2MIN)/DSX2) IF(IPOS2.LT.0.OR.NL2.LT.IPOS2) THEN C NET-48 CALL ERROR * ('NET-48: Source or receiver point outside the model') C NET-48: Source or receiver point outside the model: C Source or receiver point is situated outside the C model volume. ELSE IF(IPOS2.GE.NL2) THEN IPOS2=NL2-1 END IF DPOS2=X2-X2MIN-(FLOAT(IPOS2)+0.5)*DSX2 END IF END IF IF(NL3.EQ.1) THEN IPOS3=0 DPOS3=0. ELSE IPOS3=INT((X3-X3MIN)/DSX3) IF(IPOS3.LT.0.OR.NL3.LT.IPOS3) THEN C NET-54 CALL ERROR * ('NET-54: Source or receiver point outside the model') C NET-54: Source or receiver point outside the model. ELSE IF(IPOS3.GE.NL3) THEN IPOS3=NL3-1 END IF DPOS3=X3-X3MIN-(FLOAT(IPOS3)+0.5)*DSX3 END IF C IPOS =1+IPOS1+(IPOS2+IPOS3*NL2)*NL1 IADR=INDX(IPOS) RETURN END C C======================================================================= C C C INTEGER FUNCTION INDX(IPOS) INTEGER IPOS C C Integer function, called by many subroutines, to return the index of C the node at the given gridpoint. The subroutine does not check the C validity of the input argument. C C----------------------------------------------------------------------- C C Common block /GRID/ is required here: INCLUDE 'net.inc' C net.inc C C----------------------------------------------------------------------- C INTEGER IPOS1,IPOS2,IPOS3,I1,I2,I3,J1,J2,J3,J C C....................................................................... C IF(L4.EQ.0) THEN INDX=IPOS ELSE IPOS1=IPOS-1 IPOS2=IPOS1/NL1 IPOS3=IPOS2/NL2 IPOS1=IPOS1-IPOS2*NL1 IPOS2=IPOS2-IPOS3*NL2 J1=IPOS1/L1 J2=IPOS2/L2 J3=IPOS3/L3 J=IRAM(IND0+1+J1+N1*(J2+N2*J3)) IF(J.LE.0) THEN INDX=0 ELSE I1=IPOS1-J1*L1 I2=IPOS2-J2*L2 I3=IPOS3-J3*L3 INDX=J+I1+L1*(I2+L2*I3) END IF END IF RETURN END C C======================================================================= C INCLUDE 'ttt.for' C ttt.for INCLUDE 'error.for' C error.for INCLUDE 'sep.for' c sep.for INCLUDE 'length.for' c length.for INCLUDE 'forms.for' c forms.for INCLUDE 'eigen.for' C eigen.for C C======================================================================= C