C
C Subroutine file 'writ.for' to create the output files of complete ray C tracing C C Version: 7.00 C Date: 2013, April 19 C Coded by Ludek Klimes C C....................................................................... C C This file consists of: C WRITB...Block data subroutine defining common blocks /WRITT/, C /WRITC/ and /WRITN/ to store the input data on the names C of output files. C WRITB C WRIT1...Subroutine called when starting the complete ray tracing C program, and when starting the computation of a new C elementary wave. C WRIT1 C WRIT2...Subroutine called when starting the complete tracing of a C new ray. C WRIT2 C WRIT31..Subroutine storing the computed quantities along a ray, C see C C.R.T.5.5.1. C It is called with constant step STORE (see input data C in the file 'ray.for') of the independent variable along C the ray, and at the points of intersection with interfaces C either before and after the transformation. C WRIT31 C WRIT32..Subroutine storing the computed quantities at specified C surfaces, see C C.R.T.5.5.2. C It is called at the points of intersection of the ray C with the specified surfaces. C WRIT32 C WRIT33..Subroutine storing the computed quantities at the C endpoints of the individual elementary waves, see C C.R.T.5.5.3. C It is called at the endpoints of the elements of rays, C situated at structural interfaces. C WRIT33 C WRIT4...Subroutine storing the quantities at the initial point of C the ray, see C C.R.T.6.1. C It is called after termination of tracing the ray. C WRIT4 C CHECK...Auxiliary subroutine to WRIT4. C CHECK C WRIT5...Subroutine called after termination of the computation of C an elementary wave, and when terminating the complete ray C tracing program. C WRIT5 C FNAME...Subroutine analyzing the specified strings and composing C the filename of them. Auxiliary subroutine to WRIT1. C FNAME C WRITTR..Subroutine writing the indices of rays forming the C homogeneous triangles in the ray-parameter domain. C WRITTR C WRITBR..Subroutine writing the indices of boundary rays. C WRITBR C WPTC2..Auxiliary subroutine to subroutine WRIT2. C Subroutine WPTC2 is called when starting a new ray. C It sets the number of points along a ray to IPT=0 C in common block /POINTC/. C WPTC2 C WPTC31..Auxiliary subroutine to subroutine WRIT31. C Subroutine WPTC31 stores the quantities at individual C points along the ray into common block /POINTC/. C The quantities are then written into the unformatted C output file by subroutine APW2, C or by subroutines APW3 and APW4. C WPTC31 C WPTC4..Auxiliary subroutine to subroutine WRIT4. C Subroutine WPTC4 stores the quantities at the initial C point of the ray into common block /POINTC/. C The quantities are then written into the unformatted C output file by subroutine APW1. C WPTC4 C C....................................................................... C C Description of data files: C C Input data WRIT to specify the names of the output files with the C computed quantities: C Data (1)-(5) common to all elementary waves should always be C located in the file specified by input SEP parameter C WRIT. C If optional input SEP parameter WRITEW=' ' (default), C data (6), repeated for all source points, should also be located C in file WRIT. C If optional input SEP parameter WRITEW is specified and is not C blank, data (6) for a single source should be located in file C WRITEW. C If there are several source points, data (6) are repeatedly C used for all source points. C The data are read in by the list directed input (free format). C The strings have thus to be enclosed in apostrophes. In the list C of input data below, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). All C input variables, except K2P, are of the type CHARACTER. C The data are read by the subroutine WRIT1. C The filenames are generated using several components (strings) in C the following way: Each component (string) is divided into words, C each word being followed by just one space. Here a word is a C substring containing no space, preceded and followed by spaces or C by the beginning or end of the string. An empty word lies between C two consecutive spaces. The filename is composed of C the 1-st word of the 1-st component, C the 1-st word of the 2-nd component, C ... C the 1-st word of the last component, C the 2-nd word of the 1-st component, C the 2-nd word of the 2-nd component, C ... C the last word of the last component. C Example: C 1-st component = 'd/ .out', C 2-nd component = 'a b', C 3-rd component = '1 2', C resulting filename = 'd/a1b2.out'. C Note that only first 16 characters of each filename component are C significant. C (1) 'TEXTW' C String describing the data. Only the first 80 characters of the C string are read. C Default: TEXTW=' '. C (2) K2P,/ C K2P... If non-zero, only two-point rays are written to the files C with quantities stored along rays. Storage of quantities C at specified surfaces is not affected by this parameter. C /... An obligatory slash. C Default: K2P=0. C (3) 'FN1A','FN1I','FN1T','FN1BR',/ C 'FN1A'... String containing the first component of the name of the C output file 'CRT-R' with the computed quantities stored C along the rays, see C C.R.T.5.5.1. C Description of output file C.R.T.5.5.1 C 'FN1I'... String containing the first component of the name of the C output file 'CRT-I' with the quantities at the initial C points of rays, see C C.R.T.6.1, C related to the output file 'CRT-R' with the computed C quantities stored along the rays. The other components C of the filename are common with the related file 'CRT-R' C containing the computed quantities. C Description of output file C.R.T.6.1 C 'FN1T'... String containing the first component of the name of the C formatted output file 'CRT-T' specifying the homogeneous C triangles in the ray-parameter domain. The other C components of the filename are common with the related C file 'CRT-R' with the computed quantities along rays. C Description of file with triangles C 'FN1BR'... String containing the first component of the name of C the formatted output file 'CRT-B' specifying the indices C of boundary rays. The other components of the filename C are common with the related file 'CRT-R' containing C the computed quantities along rays. C Description of file with boundary rays C /... An obligatory slash. C Default: all strings blank. C (4) 'FN2A','FN2I',('FN2B(I)',I=1,NENDF),/ C 'FN2A'... String containing the first component of the name of the C output file 'CRT-S' with the computed quantities stored at C the specified surfaces, see C C.R.T.5.5.2. C Description of output file C.R.T.5.5.2 C 'FN2I'... String containing the first component of the name of the C output file 'CRT-I' with the quantities at the initial C points of rays, see C C.R.T.6.1, C related to the output file 'CRT-S' with the computed C quantities stored at the specified surfaces. The other C components of the filename are common with the related C file 'CRT-S' containing the computed quantities. C Description of output file C.R.T.6.1 C 'FN2B(I)'... String containing the second component of the name of C the output file 'CRT-S' with the computed quantities C stored at the I-th specified surface. Here NENDF is C the number of surfaces specified for storing the computed C quantities. C /... An obligatory slash. C Default: all strings blank. C (5) Either (5.1) or (5.2): C (5.1) / C A slash is fully sufficient if MODCRT.LE.1, see C input data set DCRT. C If MODCRT.LE.1, the filenames (5.2) are ignored if specified. C If MODCRT.GE.2, the filenames (5.2) are obligatory. C (5.2) 'FN3A','FN3I',('FN3B(I)',I=1,NELEM),/ C 'FN3A'... String containing the first component of the name of the C output file 'CRT-E' with the computed quantities stored C at the endpoints of the elements of rays, see C C.R.T.5.5.3. C Description of output file C.R.T.5.5.3 C 'FN3I'... String containing the first component of the name of the C output file 'CRT-I' with the quantities at the initial C points of rays, see C C.R.T.6.1, C related to the output file 'CRT-E' with the computed C quantities stored at the endpoints of the elements C of rays. The other components of the filename are common C with the related file 'CRT-E' containing the computed C quantities. C Description of output file C.R.T.6.1 C 'FN3B(I)'... String containing the second component of the name of C the output file 'CRT-E' with the computed quantities C stored at the endpoints of the I-th element of rays. C None of the strings FN3B(I), corresponding to the endpoint C of a ray element at which quantities are stored, may equal C spaces for MODCRT.GE.2. C /... An obligatory slash. C Default: all strings blank. C (6) For each elementary wave IWAVE the following data (6.1): C (6.1) 'FN1C','FN2C','FN3C',('FN2D(I)',I=1,NENDF),/ C 'FN1C'... String containing the second component of the name of C the output file 'CRT-R' with the computed quantities C stored along the rays, see C C.R.T.5.5.1. C 'FN2C'... String containing the third component of the name of the C output file 'CRT-S' with the computed quantities stored at C the specified surfaces, see C C.R.T.5.5.2. C 'FN3C'... String containing the third component of the name of the C output file 'CRT-E' with the computed quantities stored at C the endpoints of the elements of rays, see C C.R.T.5.5.3. C This string may not equal spaces for MODCRT.GE.2. C If MODCRT.LE.1, the file is not created. C 'FN2D(I)'... String containing the fourth component of the name of C the output file 'CRT-S' with the computed quantities C stored at the I-th specified surface, see C C.R.T.5.5.2. C Here NENDF is the number of surfaces specified for storing C the computed quantities. C /... An obligatory slash. C Default: all strings blank. C Following table summarizes the components, of which the individual C filenames of the output files are composed: C Output file: Components: C 'CRT-R' FN1A + FN1C(IW) C corresponding 'CRT-I' FN1I + FN1C(IW) C corresponding 'CRT-T' FN1T + FN1C(IW) C corresponding 'CRT-B' FN1BR+ FN1C(IW) C 'CRT-S' FN2A + FN2B(IS) + FN2C(IW) + FN2D(IS,IW) C corresponding 'CRT-I' FN2I + FN2B(IS) + FN2C(IW) + FN2D(IS,IW) C 'CRT-E' FN3A + FN3B(IE) + FN3C(IW) C corresponding 'CRT-I' FN3I + FN3B(IE) + FN3C(IW) C IW is a total index of elementary wave, IW=NWAVE*(ISRC-1)+IWAVE, C where ISRC is index of the source, IWAVE is the index of C the elementary wave for source ISRC, and NWAVE is the number of C elementary waves to be calculated for each source. C IS is an index of a storing surface, IS=1,NENDF. C IE is an index of a storing element, IE=1,NELEM. C Notes: C If the first component of the name of any of the CRT output C files is blank, the file is not created. C File 'CRT-S' and corresponding file 'CRT-I' are not created if C FN2B(IS), FN2C(IW) and FN2D(IS,IW) are all blank. C The components FN1A, FN1I, FN1T and FN1BR must be mutually C different, if they are not blank. C The components FN2A and FN2I must differ, if they are not blank. C The components FN3A and FN3I must differ, if they are not blank. C C Examples of data set WRIT: C C writ.dat two-point rays stored. C C writall.dat all rays stored. C C writsrf.dat no rays stored. C In the above examples, the computed quantities at specified surfaces C are stored. C C The names of the output files of program 'crt.for' are specified here C in the input data WRIT. On the other hand, the programs processing C the results of complete ray tracing assume the names of the same C output files of program 'crt.for' to be specified in file CRTOUT C formatted as follows: C C General description of formatted file CRTOUT: C The file consists of character strings to be read by list directed C (free format) input. The strings have thus to be enclosed in C apostrophes. C (1) For each set of CRT output files data record (1.1): C (1.1) 'CRT-R','CRT-S','CRT-I','CRT-T','CRT-B',/ C 'CRT-R'.. Name of the unformatted CRT output file with the C quantities stored along rays, see C C.R.T.5.5.1. C Description of file C 'CRT-R'. C 'CRT-S'.. Name of the unformatted CRT output file with the C quantities stored at the points of intersections of rays C with the specified surfaces, see C C.R.T.5.5.2. C Description of file C 'CRT-S'. C 'CRT-I'.. Name of the unformatted CRT output file with the C quantities at the initial points of rays, see C C.R.T.6.1. C Description of file C 'CRT-I'. C 'CRT-T'.. Name of the formatted CRT output file with the C homogeneous triangles in the ray-parameter domain. C Description of file C 'CRT-T'. C 'CRT-B'.. Name of the formatted CRT output file with the C indices of boundary rays. C Description of file C 'CRT-B'. C /... An obligatory slash. C General defaults for the first set of the CRT output files: C The defaults may be blank or following: C 'CRT-R'='r01.out', 'CRT-S'='s01.out', 'CRT-I'='r01i.out' C or 'CRT-I'='s01i.out', 'CRT-T'='t01.out', 'CRT-B'=' ' C General defaults for the subsequent sets: C 'CRT-R'=' ', 'CRT-S'=' ', 'CRT-I'=' ', 'CRT-T'=' ', C 'CRT-B'=' ' C (2) / (a slash or the end of file). C The end of the data is also assumed if all filenames are blank. C C C Output files 'CRT-R' with the computed quantities stored along the C rays, see C.R.T.5.5.1. C Unformatted output. The computed quantities are stored with the C step STORE (see the input data in the subroutine file 'ray.for') C along all rays, and at all points of intersection of rays with C structural interfaces either before and after the transformation. C For each point - one record containing the following quantities: C (1) ISRC,IWAVE,IRAY,NY,ICB1,ISRF,X,(YL(I),I=1,6),(Y(I),I=1,NY) C ISRC .. Index of the source. C IWAVE...Index of the elementary wave (output of the subroutine C CODE1). Elementary waves are indexed by 1,2,3,... . C IRAY... Index of the ray (output of the subroutine RPAR2). Rays C within each elementary wave are indexed by 1,2,3,... . C Note that some of the indexed rays need not exist. C NY=IY(1)=27+NAMPL... Number of the basic quantities describing the C point of the ray, see the file 'ray.for'. C ICB1=IY(5)... Index of the complex block in which the stored point C of the ray is situated, supplemented by a sign '+' for P C wave and sign '-' for S wave. C ISRF=IY(6)... Index of the surface at which the endpoint of the C computed element of the ray is situated, supplemented by C a sign '+' or '-' for the endpoint situated at the C positive or negative side of the surface, respectively. C Zero inside the complex block. Note: the sign of this C quantity is the exception to the definition in the C original paper on C.R.T. C X... Independent variable along a ray, see YY(1) in C C.R.T.5.2.2. C YL... Array containing local quantities at the point of the ray, C see C C.R.T.5.5.4. C C Description of YL. C Y... Array containing basic quantities computed along the ray, C see C C.R.T.5.2.1. C Description of Y. C The detailed description of the quantities YL, YY and IY may be C found in the file 'ray.for'. C C C Output files 'CRT-S' with the computed quantities stored at the C specified surfaces, see C C.R.T.5.5.2. C Unformatted output. The computed quantities are stored at the C points of intersection of rays with the specified surfaces. C For each point - one record containing the following quantities: C (1) ISRC,IWAVE,IRAY,NY,ICB1,ISRF,X,(YL(I),I=1,6),(Y(I),I=1,NY) C The same quantities as in the data set described above, except C that the vectorial reduced amplitudes Y(28)-Y(IY(1)) may (but need C not) be replaced by the reduced amplitudes YC(1) to YC(NY-27) C involving appropriate conversion coefficients, see C C.R.T.5.5.4. C YC(1) to YC(NY-27) are described in the subroutine CONV (file C 'trans.for') and in the subroutine WRIT32 (this file). Note that C NY=33 for P wave and NY=39 for S wave at the initial point of the C ray. C C C Output files 'CRT-E' with the computed quantities stored at the C endpoints of the elements of rays, see C C.R.T.5.5.3. C Unformatted output. The files, generated only if KSTORE.GE.2 (see C input data in the subroutine file 'ray.for'), are auxiliary disk C storage locations to the program CRT and are not intended to be C the output of the complete ray tracing used by the user C application programs. The computed quantities are stored at the C endpoints of the elements of all rays before the transformation. C For each point - one record containing the following quantities: C (1) IRAY,(IY(I),I=1,12),(YL(I),I=1,6),(Y(I),I=1,IY(1)),(YY(I),I=1,5) C IRAY... Index of the ray (output of the subroutine RPAR2). Rays C within each elementary wave are indexed by 1,2,3,... . C Note that some of the indexed rays may not exist. C YL... Array containing local quantities at the point of the ray, C see C C.R.T.5.5.4. C C Description of YL. C Y... Array containing basic quantities computed along the ray, C see C C.R.T.5.2.1. C Description of Y. C YY... Array containing real auxiliary quantities computed along C the ray, see C C.R.T.5.2.2. C C Description of YY. C IY... Array containing integer auxiliary quantities computed C along the ray, see C C.R.T.5.2.2. C C Description of IY. C The detailed description of the quantities YL, Y, YY and IY may be C found in the file 'ray.for'. C C C Output files 'CRT-I' with the quantities at the initial points of C rays, see C.R.T.6.1. C Unformatted output. The quantities are stored at the initial C points of rays. C For whole rays ('CRT-R'), only initial points of two-point C rays are stored if only two-point rays are stored. C For intersections with surfaces ('CRT-S') and endpoints of the C elements of rays ('CRT-E') initial points of all rays C are stored. C For initial point - one record containing following 34 quantities: C (1) ISRC,-IWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,(YLI(I),I=1,6), C (YI(I),I=1,29) C ISRC .. Index of the source. C IWAVE...Index of the elementary wave (output of the subroutine C CODE1). Elementary waves are indexed by 1,2,3,... . C IRAY... Index of the ray (output of the subroutine RPAR2). Rays C within each elementary wave are indexed by 1,2,3,... . C Note that some of the indexed rays may not exist. C ICB1I...Index of the complex block in which the initial point of C the ray is situated, supplemented by a sign '+' for P wave C and sign '-' for S wave, see C C.R.T.6.1. C IEND... Reason of the termination of the computation of a ray, see C C.R.T.5.4. C See subroutine RAY in the subroutine file 'ray.for' for C detailed description of C IEND. C ISHEET..Ray-history index. The different ray histories are C consecutively indexed by positive integers 1,2,3,... C According to their appearance during ray tracing. C The ray histories are indexed independently within each C elementary wave. The indices are the output of subroutine C RPAR4 of the file 'rpar.for'. C The ray-history indices are complemented with sign: C Positive - successful ray (crossing reference surface), C negative - unsuccessful ray (terminating before crossing C reference surface). C IREC... Index of the receiver for a two-point ray, determined in C subroutine C RPAR4. C YLI... Array containing the values of the quantities YL(1)-YL(6), C see C C.R.T.5.5.4, C describing the local properties of the model at the C initial point of the ray. See the list of YL(1) to YL(6) C in the file 'ray.for'. C C Description of YL. C YI... Array containing the quantities describing the properties C of the rays and of the travel-time field at the initial C point of the ray, see C C.R.T.6.1. C These quantities are listed in the file 'init.for'. C C Description of YI. C C C Output files 'CRT-T' specifying homogeneous triangles in the C ray-parameter domain. Formatted output, each record containing C three integers: C (1) For each elementary wave IWAVE the following data (1.1) and (1.2): C (1.1) 0,ISRC,IWAVE C 0 ... Zero identifying new elementary wave. C ISRC .. Index of the source. C IWAVE . Index of the elementary wave. C (1.2) For each homogeneous triangle of the elementary wave C data (1.2.1): C (1.2.1) IRAY1,IRAY2,IRAY3 C IRAY1,IRAY2,IRAY3... Indices of rays forming the vertices of the C homogeneous triangle. C For one-parametric (2-D) systems of rays, triangles are C replaced with intervals and IRAY3=0. C C C Output file 'CRT-B' specifying indices of boundary rays. C Formatted output, each record containing four integers: C (1) For each elementary wave IWAVE the following data (1.1) and (1.2): C (1.1) 0,0,ISRC,IWAVE C 0 ... Zero identifying new elementary wave. C 0 ... Zero. C ISRC .. Index of the source. C IWAVE . Index of the elementary wave. C (1.2) For each boundary ray data (1.2.1): C (1.2.1) IB1,IB2,IB3,IB4 C IB1 ... Index of the boundary ray. C IB2 ... Index of the boundary ray coupled with IB1. C IB3,IB4 . Indices of two boundary rays equal in history with IB1, C coupled with rays equal in history with IB2, and connected C with boundary ray IB1 by sides of homogeneous triangles. C IB4=0 if there is only one ray of the described C properties. C IB3=0 and IB4=0 if there is no ray of the described C properties, and for one-parametric shooting in 2-D. C C....................................................................... C C Storage in the memory: C The input data WRIT (2) to (6) describing the names of the output C files with the computed quantities are stored in the common block C /WRITT/. Other important variables shared by the subroutines of C this file are stored in the common blocks /WRITC/ and /WRITN/. C The common blocks are defined in the following subroutine: C ------------------------------------------------------------------ C C BLOCK DATA WRITB INCLUDE 'writ.inc' C writ.inc DATA LUW/10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25, * 26,27,28,29,30,31,32,33,34,35,36,37,38,39/ DATA JENDR/10,21,22,23,24,25,26,30,32,33,37,39, * 40,50,60,61,71,72,73,74/ END C ------------------------------------------------------------------ C C======================================================================= C C C SUBROUTINE WRIT1(LUN,LULOG,IWAVE,IWAVE0,IKODE) INTEGER LUN,LULOG,IWAVE,IWAVE0,IKODE C C This subroutine is called when starting the complete ray tracing C program, and when starting the computation of a new elementary wave. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C LULOG...Logical unit number of the log output device. C IWAVE...Zero when starting the complete ray tracing program, C otherwise the index of the elementary wave which will be C computed (i.e. the output of the subroutine CODE1 from the C file 'code.for'). C IWAVE0..Index of the already computed elementary wave having the C most numerous common elements with the current elementary C wave. Need not be defined if IWAVE=0. C IKODE...The length of the common part of the codes of the IWAVE-th C and IWAVE0-th elementary waves. Need not be defined if C IWAVE=0. C C No output. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common blocks /WRITT/, /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C All the storage locations of the common blocks are defined in this C subroutine. C C Subroutines and external functions required: INTEGER LCODE,NELEM EXTERNAL WRITB,LCODE,NELEM,SCRO1,FNAME C WRITB.. Block data subroutine of this file. C LCODE,NELEM... File 'code.for'. C SCRO1...One of files 'scro*.for'. C FNAME... This file. C C Date: 2005, June 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER NKODE,IE1,IE2,I,J CHARACTER*16 NAME1(4) C IF(IWAVE.LE.0) THEN C C Initiating the source index ISRC=1 C C Reading the name of the input data TEXTW=' ' READ(LUN,*) TEXTW C K2P=0 READ(LUN,*) K2P C C Reading strings composing the names of output files (CRT.5.5.1) FN1A=' ' FN1I=' ' FN1T=' ' FN1BR=' ' READ(LUN,*) FN1A,FN1I,FN1T,FN1BR C Number of output files to be open IF(STORE.EQ.0.) THEN NF1=0 ELSE NF1=1 END IF C C Reading strings composing the names of output files (CRT.5.5.2) FN2A=' ' FN2I=' ' DO 21 I=1,MF FNB(I)=' ' 21 CONTINUE READ(LUN,*) FN2A,FN2I,(FNB(I),I=1,MF) C Number of output files to be open NF2=NSTOR IF(2*(NF1+NF2).GT.NUW) THEN C 553.1 CALL ERROR('553.1 in WRIT1: Few logical units available') C There are NUW logical units available, see the common C blocks /WRITC/ and /WRITN/ defined in the block data C WRITB. This number should be increased. END IF C C Reading strings composing the names of output files (CRT.5.5.3) FN3A=' ' FN3I=' ' DO 31 I=NF2+1,MF FNB(I)=' ' 31 CONTINUE READ(LUN,*) FN3A,FN3I,(FNB(I),I=NF2+1,MF) C Number of output files to be open IF(MODCRT.LE.1) THEN MF3=0 ELSE DO 32 I=NF2+1,MF IF(FNB(I).EQ.' ') THEN MF3=I-1-NF2 GO TO 33 END IF 32 CONTINUE MF3=MF-NF2 33 CONTINUE END IF C IF ( ((FN1A .NE.' ').AND. * ((FN1A .EQ.FN1I).OR.(FN1A .EQ.FN1T).OR.(FN1A .EQ.FN1BR))).OR. * ((FN1I .NE.' ').AND. * ((FN1I .EQ.FN1A).OR.(FN1I .EQ.FN1T).OR.(FN1I .EQ.FN1BR))).OR. * ((FN1T .NE.' ').AND. * ((FN1T .EQ.FN1A).OR.(FN1T .EQ.FN1I).OR.(FN1T .EQ.FN1BR))).OR. * ((FN1BR.NE.' ').AND. * ((FN1BR.EQ.FN1A).OR.(FN1BR.EQ.FN1I).OR.(FN1BR.EQ.FN1T ))).OR. * ((FN2A.NE.' ').AND.(FN2A.EQ.FN2I)).OR. * ((FN3A.NE.' ').AND.(FN3A.EQ.FN3I)) ) THEN C 554 CALL ERROR('554 in WRIT1: Equal first filename components.') C The components FN1A, FN1I, FN1T and FN1BR must be mutually C different, if they are not blank. Similarly FN2A and FN2I must C differ, and FN3A and FN3I must differ, if they are not blank. ENDIF C ELSE C C Updating the source index IF(IWAVE.LE.JWAVE) THEN ISRC=ISRC+1 END IF C C Numbers of initial ray elements along which nothing is stored IF(MODCRT.EQ.0) THEN KWRIT1=0 ELSE KWRIT1=NELEM(IKODE) END IF KWRIT2=KWRIT1 IF(MODCRT.LE.1) THEN KWRIT3=999999 ELSE KWRIT3=NELEM(IKODE) END IF C C Reading strings composing the names of output files DO 41 I=JWAVE+1,IWAVE-1 READ(LUN,*,END=98) FN1C,FN2C,FN3C,(FN2D(J),J=1,MF) 41 CONTINUE FN1C=' ' FN2C=' ' FN3C=' ' DO 42 I=1,MF FN2D(I)=' ' 42 CONTINUE READ(LUN,*,END=98) FN1C,FN2C,FN3C,(FN2D(I),I=1,MF) C C Opening the output files (CRT.5.5.1) IF(NF1.NE.0.) THEN NAME1(1)=FN1A NAME1(2)=FN1C NAME2(1)=' ' IF(FN1A.NE.' ') THEN CALL FNAME(2,NAME1,NAME2(1)) OPEN(LUW(1),FILE=NAME2(1),FORM='UNFORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C 551.1 C Open file error IE2=1 J=1 GO TO 91 END IF END IF NAME1(1)=FN1I NAME2(2)=' ' IF(FN1I.NE.' ') THEN CALL FNAME(2,NAME1,NAME2(2)) OPEN(LUW(2),FILE=NAME2(2),FORM='UNFORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C 551.4 C Open file error IE2=4 J=2 GO TO 91 END IF END IF END IF C C Opening the output files (CRT.5.5.2) DO 44 I=1,NF2 NAME1(1)=FN2A NAME1(2)=FNB(I) NAME1(3)=FN2C NAME1(4)=FN2D(I) IF(NAME1(2).EQ.' '.AND.NAME1(3).EQ.' '.AND.NAME1(4).EQ.' ') * THEN NAME1(1)=' ' END IF J=2*(NF1+I)-1 NAME2(J)=' ' IF(FN2A.NE.' ') THEN CALL FNAME(4,NAME1,NAME2(J)) IF(NAME2(J).NE.' ') THEN OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C 551.2 C Open file error IE2=2 GO TO 91 END IF END IF END IF NAME1(1)=FN2I IF(NAME1(2).EQ.' '.AND.NAME1(3).EQ.' '.AND.NAME1(4).EQ.' ') * THEN NAME1(1)=' ' END IF J=2*(NF1+I) NAME2(J)=' ' IF(FN2I.NE.' ') THEN CALL FNAME(4,NAME1,NAME2(J)) IF(NAME2(J).NE.' ') THEN OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C 551.5 C Open file error IE2=5 GO TO 91 END IF END IF END IF 44 CONTINUE C C Opening the output files (CRT.5.5.3) NF3=0 IF(MODCRT.GE.2) THEN IF(FN3C.EQ.' ') THEN C 557 CALL ERROR('557 in WRIT1: Few filenames specified') C There are more required elementary waves than the C corresponding specified filenames. ELSE NKODE=NELEM(LCODE()) C NKODE is the number of ray elements IF(MODCRT.GE.3) NKODE=MIN0(KWRIT3+1,NKODE) C NKODE is the index of the last ray element to be stored NF3=NKODE-KWRIT3 IF(NF3.GT.MF3) THEN C 556 CALL ERROR('556 in WRIT1: Few filenames specified') C There are more elements of the rays of an elementary wave C than the corresponding specified filenames. ELSE IF(2*(NF1+NF2+NF3).GT.NUW) THEN C 553.2 CALL ERROR('553.2 in WRIT1: Few logical units available') C There are NUW logical units available, see the common C blocks /WRITC/ and /WRITN/ defined in the block data C WRITB. This number should be increased. END IF DO 45 I=1,NF3 NAME1(1)=FN3A NAME1(2)=FNB(KWRIT3+I) NAME1(3)=FN3C J=2*(NF1+NF2+I)-1 NAME2(J)=' ' IF(FN3A.NE.' ') THEN CALL FNAME(3,NAME1,NAME2(J)) OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C C 551.3 C Open file error IE2=3 GO TO 91 END IF END IF NAME1(1)=FN3I J=2*(NF1+NF2+I) NAME2(J)=' ' IF(FN3I.NE.' ') THEN CALL FNAME(3,NAME1,NAME2(J)) OPEN(LUW(J),FILE=NAME2(J),FORM='UNFORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C C 551.6 C Open file error IE2=6 GO TO 91 END IF END IF 45 CONTINUE END IF END IF C C Opening the output file with homogeneous triangles: J=2*(NF1+NF2+NF3)+1 IF (J.GT.NUW) THEN C 553.3 CALL ERROR('553.3 in WRIT1: Few logical units available') C There are NUW logical units available, see the common C blocks /WRITC/ and /WRITN/ defined in the block data C WRITB. The number NUW should be increased. ENDIF NAME2(J)=' ' IF(FN1T.NE.' ') THEN NAME1(1)=FN1T NAME1(2)=FN1C CALL FNAME(2,NAME1,NAME2(J)) OPEN(LUW(J),FILE=NAME2(J),FORM='FORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C 551.7 C Open file error IE2=7 GO TO 91 END IF END IF C C Opening the output file with boundary rays: J=2*(NF1+NF2+NF3)+2 IF (J.GT.NUW) THEN C 553.4 CALL ERROR('553.4 in WRIT1: Few logical units available') C There are NUW logical units available, see the common C blocks /WRITC/ and /WRITN/ defined in the block data C WRITB. The number NUW should be increased. ENDIF NAME2(J)=' ' IF(FN1BR.NE.' ') THEN NAME1(1)=FN1BR NAME1(2)=FN1C CALL FNAME(2,NAME1,NAME2(J)) OPEN(LUW(J),FILE=NAME2(J),FORM='FORMATTED',IOSTAT=IE1) IF(IE1.NE.0) THEN C 551.8 C Open file error IE2=8 GO TO 91 END IF END IF C C Initial values for the elementary wave JRAY=0 JUEB=0 JFCT=0 JOUTP=0 JTRANS=0 DO 51 I=1,MENDR+2 NENDR(I)=0 51 CONTINUE DO 52 I=1,2*(NF1+NF2+NF3) JPOINT(I)=0 52 CONTINUE C C Writing to the output LOG file WRITE(LULOG,81) IWAVE 81 FORMAT(' WAVE',I5,':') C C Writing to the file 'CRT-T' with homogeneous triangles: CALL WRITTR(0,ISRC,IWAVE) C C Writing to the file 'CRT-B' with boundary rays: CALL WRITBR(0,0,ISRC,IWAVE) C END IF C C Index of the last wave: JWAVE=IWAVE C C Screen output CALL SCRO1(ISRC,IWAVE) RETURN C 91 CONTINUE C 551 WRITE(*,'('' STATUS'',I6,''.'',I1,'': '',A)') IE1,IE2,NAME2(J) CALL ERROR('551 in WRIT1: Open file error') C Open file error when opening the file to store: C 1 computed quantities along the rays. C Location of error 551.1 in the code C 2 computed quantities at the specified surfaces. C Location of error 551.2 in the code C 3 computed quantities at the endpoints of the elements of C the rays of the elementary wave. C Location of error 551.3 in the code C 4 quantities at the initial points of rays, corresponding C to the above file 1. C Location of error 551.4 in the code C 5 quantities at the initial points of rays, corresponding C to the above file 2. C Location of error 551.5 in the code C 6 quantities at the initial points of rays, corresponding C to the above file 3. C Location of error 551.6 in the code C 7 homogeneous triangles in the ray-parameter domain. C Location of error 551.7 in the code C 8 boundary rays. C Location of error 551.8 in the code C The type 1 to 8 of the file is given by the first decimal C place of the status. 98 CONTINUE C 558 CALL ERROR('558 in WRIT1: End of input file') C End of input data file WRIT specifying the names of the C output files encountered before all required elementary C waves are computed. The number of elementary waves C exceeds the number of specified output filenames. END C C======================================================================= C C C SUBROUTINE WRIT2(LULOG,IRAY) INTEGER LULOG,IRAY C C This subroutine is called when starting the complete tracing of a new C ray. C C Input: C LULOG...Logical unit number of the log output device. C IRAY... The index of the ray which will be computed (i.e. the C output of the subroutine RPAR2 from the file 'rpar.for'). C C No output. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C JRAY is redefined in this subroutine. No other of the storage C locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL SCRO2 C SCRO2...One of files 'scro*.for'. C C Date: 2013, February 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C NTMP=0 JRAY=IRAY DO 1 I=1,NF2 JSTOR(I)=0 1 CONTINUE C C Setting the number of points along a ray in common block /POINTC/ C to IPT=0: CALL WPTC2() C C Screen output CALL SCRO2(IRAY) RETURN END C C======================================================================= C C C SUBROUTINE WRIT31(YL,Y,YY,IY) REAL YL(6),Y(44),YY(5) INTEGER IY(13) C C This subroutine stores the computed quantities along a ray, see C C.R.T.5.5.1. C It is called with constant step STORE of the independent variable C along the ray, and at the points of intersection with interfaces C either before and after the transformation. C C Input: C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. C None of the input parameters are altered. C C No output. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C Array JPOINT is modified in this subroutine. C C Subroutines and external functions required: EXTERNAL NELEM,SCRO3 INTEGER NELEM C NELEM.. File 'code.for'. C SCRO3...One of files 'scro*.for'. C C Date: 2013, February 21 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C IF(NELEM(IY(2)).GT.KWRIT1) THEN IF(NAME2(1).NE.' ') THEN CALL WPTC31(ISRC,JWAVE,JRAY,IY(1),IY(5),IY(6),IY(12), * YY(1),YY(2),YL,Y,IY(13),Y(36)) IF(K2P.EQ.0) THEN JPOINT(1)=JPOINT(1)+1 CALL APW2(LUW(1)) ELSE NTMP=NTMP+1 CALL APW3(KTMP,MKTMP,TMP,MTMP) END IF END IF END IF C C Screen output CALL SCRO3(YL,Y,YY,IY) RETURN END C C======================================================================= C C C SUBROUTINE WRIT32(ISTOR,YL,Y,YY,IY,NAMPL,YC) INTEGER ISTOR,IY(13),NAMPL REAL YL(6),Y(44),YY(5),YC(NAMPL) C C This subroutine stores the computed quantities at specified surfaces, C see C.R.T.5.5.2. C It is called at the points of intersection of the ray with specified C surfaces. C C Input: C ISTOR...The sequential number in the input data of the specified C surface. C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C Quantities Y(28) to Y(IY(1)) representing reduced C amplitudes are ignored. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. C NAMPL...Number of real quantities representing complex-valued C vectorial reduced amplitudes. If no conversion C coefficients are applied NAMPL=IY(1)-27, otherwise C NAMPL=6 or 12, see C C.R.T.5.5.4. C YC... Array containing real quantities representing complex- C -valued vectorial reduced amplitudes. If no conversion C coefficients are applied, YC is a copy of Y(28) to C Y(IY(1)). Otherwise, YC represents the vectorial reduced C amplitudes involving appropriate conversion coefficients, C expressed in ray-centred coordinate system, see C C.R.T.5.5.4. C P wave at the initial point of the ray: C NAMPL=6, C YC(1)=REAL(A13), YC(2)=AIMAG(A13), C YC(3)=REAL(A23), YC(4)=AIMAG(A23), C YC(5)=REAL(A33), YC(6)=AIMAG(A33). C S wave at the initial point of the ray: C NAMPL=12, C YC(1)=REAL(A11), YC(2)=AIMAG(A11), C YC(3)=REAL(A21), YC(4)=AIMAG(A21), C YC(5)=REAL(A31), YC(6)=AIMAG(A31), C YC(7)=REAL(A12), YC(8)=AIMAG(A12), C YC(9)=REAL(A22), YC(10)=AIMAG(A22), C YC(11)=REAL(A32), YC(12)=AIMAG(A32). C None of the input parameters are altered. C C No output. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C Array JPOINT is modified in this subroutine. C C Subroutines and external functions required: EXTERNAL NELEM INTEGER NELEM C NELEM...File 'code.for'. C C Date: 2013, February 21 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER NY C JSTOR(ISTOR)=JSTOR(ISTOR)+1 IF(KSTOR(NSTOR+ISTOR).LE.JSTOR(ISTOR) * .AND.JSTOR(ISTOR).LE.KSTOR(2*NSTOR+ISTOR)) THEN IF(NELEM(IY(2)).GT.KWRIT2) THEN IF(NAME2(2*(NF1+ISTOR)-1).NE.' ') THEN JPOINT(2*(NF1+ISTOR)-1)=JPOINT(2*(NF1+ISTOR)-1)+1 NY=27+NAMPL CALL WPTC32(ISRC,JWAVE,JRAY,NY,IY(5),IY(6),IY(12), * YY(1),YY(2),YL,Y,YC,IY(13),Y(36)) CALL APW2(LUW(2*(NF1+ISTOR)-1)) END IF END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE WRIT33(YL,Y,YY,IY) REAL YL(6),Y(44),YY(5) INTEGER IY(13) C C This subroutine stores the computed quantities at the endpoints of the C individual elementary waves, see C C.R.T.5.5.3. C It is called at the endpoints of the elements of rays, situated C at structural interfaces. C C Input: C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. C None of the input parameters are altered. C C No output. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL NELEM INTEGER NELEM C NELEM...File 'code.for'. C C Date: 2013, February 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage location: INTEGER I1,I2,I C I=NELEM(IY(2))-KWRIT3 IF(1.LE.I.AND.I.LE.NF3) THEN IF(NAME2(2*(NF1+NF2+I)-1).NE.' ') THEN JPOINT(2*(NF1+NF2+I)-1)=JPOINT(2*(NF1+NF2+I)-1)+1 IF(IY(13).EQ.3) THEN I1=7 I2=9 ELSE I1=1 I2=IY(13) END IF WRITE(LUW(2*(NF1+NF2+I)-1)) ISRC,JWAVE,JRAY,IY, * YL,(Y(I),I=1,IY(1)), * (Y(I),I=35+I1,35+I2),YY C These data are not read in the current versin of package CRT. END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE WRIT4(LULOG,IRAY,YL,Y,YY,IY,IEND,ISHEET,IREC) C INTEGER LULOG,IRAY,IY(13),IEND,ISHEET,IREC REAL YL(6),Y(44),YY(5) C C This subroutine stores the quantities at the initial point of the ray, C see C.R.T.6.1. C It is called after termination of tracing the ray. C C Input: C LULOG...Logical unit number of the log output device. C IRAY... The index of the ray which has been computed (i.e. the C output of the subroutine RPAR2 from the file 'rpar.for'). C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. C IEND... Reason of the termination of the computation of a ray, see C C.R.T.5.4. C See subroutine RAY in the subroutine file 'ray.for' for C detailed description of C IEND. C ISHEET..Ray-history index. The different ray histories are C consecutively indexed by positive integers 1,2,3,... C According to their appearance during ray tracing. C The ray histories are indexed independently within each C elementary wave. C The ray-history indices are complemented with sign: C positive - successful ray (crossing reference surface), C negative - unsuccessful ray (terminating before crossing C reference surface). C IREC... Index of the receiver for a two-point ray, C otherwise zero. C C No output. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /INITC/ (see subroutine file 'init.for'): INCLUDE 'initc.inc' C initc.inc C None of the storage locations are altered. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C JUEB, JFCT, JOUTP, JTRANS, arrays JPOINT and NENDR are modified in C this subroutine. C C Subroutines and external functions required: EXTERNAL SCRO4,CHECK C SCRO4...One of files 'scro*.for'. C CHECK...This file. C C Date: 2013, April 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,I0,I1,I2,I3,I4,I5,I6,I7,I8,I9,J REAL S12,S13,S23,S14,S24,S34 C C....................................................................... C C Screen output CALL SCRO4(IRAY,YL,Y,YY,IY,IEND,ISHEET) C C Statistics DO 11 I=1,MENDR IF(IEND.EQ.JENDR(I)) THEN NENDR(I)=NENDR(I)+1 GO TO 12 END IF 11 CONTINUE NENDR(MENDR+1)=NENDR(MENDR+1)+1 12 CONTINUE NENDR(MENDR+2)=NENDR(MENDR+2)+1 C C Writing the quantities at the initial point of the ray, see C C.R.T.6.1. DO 20 I=1,NF1+NF2+NF3 IF(NAME2(2*I).NE.' ') THEN C For whole rays, only initial points of two-point rays are C stored if only two-point rays are stored. C For intersections with surfaces and endpoints of elementary C waves, initial points of all rays are stored. IF(K2P.EQ.0.OR.IREC.GT.0.OR.I.GT.NF1) THEN JPOINT(2*I)=JPOINT(2*I)+1 CALL WPTC4(ISRC,JWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI, * NPVI,YI(30)) CALL APW1(LUW(2*I)) END IF END IF 20 CONTINUE C C Check for non-existing rays: IF(IEND.GE.70) THEN RETURN END IF C C Statistics JFCT=JFCT+IY(9) JOUTP=JOUTP+IY(10) JTRANS=JTRANS+IY(11) C C Writing the quantities along the two-point ray: IF(NAME2(1).NE.' ') THEN IF(K2P.NE.0.AND.IREC.GT.0) THEN JPOINT(1)=JPOINT(1)+NTMP CALL APW4(LUW(1),KTMP,TMP) END IF KTMP(1)=0 END IF C C Writing to the output LOG file C Simplecticity of the propagator matrix, see C C.R.T.5.6l. IF(UEBDRT.GT.0.) THEN S12=Y(12)*Y(18)+Y(13)*Y(19)-Y(14)*Y(16)-Y(15)*Y(17) S13=Y(12)*Y(22)+Y(13)*Y(23)-Y(14)*Y(20)-Y(15)*Y(21)-1. S23=Y(16)*Y(22)+Y(17)*Y(23)-Y(18)*Y(20)-Y(19)*Y(21) S14=Y(12)*Y(26)+Y(13)*Y(27)-Y(14)*Y(24)-Y(15)*Y(25) S24=Y(16)*Y(26)+Y(17)*Y(27)-Y(18)*Y(24)-Y(19)*Y(25)-1. S34=Y(20)*Y(26)+Y(21)*Y(27)-Y(22)*Y(24)-Y(23)*Y(25) END IF C Checking exceeding the specified limits I=0 CALL CHECK(YY(2),UEB,I0,I) CALL CHECK(YY(3),UEBPP,I1,I) CALL CHECK(YY(4),UEBPH,I2,I) CALL CHECK(YY(5),UEBHH,I3,I) CALL CHECK(S12,UEBDRT,I4,I) CALL CHECK(S13,UEBDRT,I5,I) CALL CHECK(S23,UEBDRT,I6,I) CALL CHECK(S14,UEBDRT,I7,I) CALL CHECK(S24,UEBDRT,I8,I) CALL CHECK(S34,UEBDRT,I9,I) C I is the number of checked quantities exceeding their specified C Upper limits IF(I.GT.0) THEN C Writing report on this ray IF(JUEB.EQ.0) THEN WRITE(LULOG,'(2A)') ' RAY: ', * 'CHECKED QUANTITIES IN PERCENTS OF THEIR SPECIFIED LIMITS' END IF JUEB=JUEB+1 WRITE(LULOG,'(I10,A,10I6)') * IRAY,': ',I0,I1,I2,I3,I4,I5,I6,I7,I8,I9 END IF C RETURN END C C----------------------------------------------------------------------- C C C SUBROUTINE CHECK(Q,QUEB,IRATE,I) REAL Q,QUEB INTEGER IRATE,I C C Auxiliary subroutine to WRIT4 C IF(QUEB.LE.0.) THEN IRATE=0 ELSE IF(ABS(Q).GT.QUEB) THEN I=I+1 END IF IF(ABS(Q).GT.10000.*QUEB) THEN IRATE=1000000. ELSE IRATE=INT(100.*Q/QUEB+0.5) END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE WRIT5(LULOG,IWAVE) INTEGER LULOG,IWAVE C C This subroutine is called after termination of the computation of an C elementary wave, and when terminating the complete ray tracing C program. C C Input: C LULOG...Logical unit number of the log output device. C IWAVE...Zero when terminating the complete ray tracing program, C otherwise the index of the elementary wave which has been C computed (i.e. the output of the subroutine CODE1 from the C file 'code.for'). C C No output. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL SCRO5 C SCRO5...One of files 'scro*.for'. C C Date: 1992, December 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,J,N,K(4) C C Writing to the output LOG file: IF(IWAVE.NE.0) THEN IF(JUEB.GT.0) THEN WRITE(LULOG,10) 10 FORMAT(11X,'ABOVE RAYS ARE COMPUTED WITH DECREASED ACCURACY') END IF C Reasons of termination of tracing rays WRITE(LULOG,21) NENDR(MENDR+2),JFCT,JOUTP,JTRANS N=0 DO 12 J=1,MENDR IF(NENDR(J).GT.0) THEN N=N+1 K(N)=J END IF IF(N.EQ.4.OR.(J.EQ.MENDR.AND.N.GT.0)) THEN WRITE(LULOG,22) (NENDR(K(I)),JENDR(K(I)),I=1,N) N=0 END IF 12 CONTINUE IF(NENDR(MENDR+1).NE.0) THEN WRITE(LULOG,23) NENDR(MENDR+1) END IF C List of output files DO 14 I=1,2*(NF1+NF2+NF3) IF(NAME2(I).NE.' ') THEN WRITE(LULOG,24) JPOINT(I),NAME2(I) END IF 14 CONTINUE C formats 21 FORMAT( I10,' RAYS ',I12,' FCT ',I12,' STEPS',I12,' TRANS') 22 FORMAT(4(I10,' IEND=',I2)) 23 FORMAT( I10,' IEND=OTHER') 24 FORMAT( I10,' POINTS IN FILE: ',A) END IF C C Screen output: CALL SCRO5(IWAVE) RETURN END C C======================================================================= C C C SUBROUTINE FNAME(NUM1,NAME1,NAME2) INTEGER NUM1 CHARACTER*(*) NAME1(NUM1) CHARACTER*(*) NAME2 C C This subroutine is designed to analyze the specified strings and to C compose the filename of them. C C Input: C NAME1...Character array of character strings to be analyzed. C NUM1... Number of filenames to be analyzed. C C Output: C NAME2...Filename composed of the analyzed input components C (strings). C C No subroutines and external functions required. C C Date: 1992, December 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER M1 PARAMETER (M1=4) INTEGER J1(M1),L1(M1),I1,J2,L2,I2,N,I C C M1... Maximum number of analyzed strings. C J1(N)...Number of the currently analyzed characters of NAME1(N), C i.e. the position of the space after the last analyzed C word of NAME1(N). A word is a substring between two C consecutive spaces. C L1(N)...Length of NAME1(N). C I1... Beginning of the last analyzed word of NAME1(N), i.e. the C previous value of J1(N) increased by 1. C J2... Number of the currently defined characters of NAME2. C L2... Length of NAME2. C I2... Auxiliary variable - previous value of J2 increased by 1. C N... Loop variable - index of input filename. C I... Loop variable - position in a string. C C....................................................................... C IF(M1.LT.NUM1) THEN C 559 CALL ERROR('559 in FNAME: Too many strings') C The number NUM1 of input strings exceeds the dimension M1 C of arrays J1 and L1. Parameter M1 must be increased. END IF DO 11 N=1,NUM1 J1(N)=0 L1(N)=LEN(NAME1(N)) 11 CONTINUE J2=0 L2=LEN(NAME2) C 20 CONTINUE DO 29 N=1,NUM1 C (a) Analyzing NAME1(N): I1=J1(N)+1 DO 21 I=I1,L1(N) IF(NAME1(N)(I:I).EQ.' ') THEN J1(N)=I GO TO 22 END IF 21 CONTINUE J1(N)=L1(N)+1 22 CONTINUE C J1(N)-th character of NAME1(N) is found to be blank. C (b) Copying the word from NAME1(N) to NAME2: IF(J1(N).GT.I1) THEN I2=J2+1 J2=MIN0(J2+J1(N)-I1,L2) NAME2(I2:J2)=NAME1(N)(I1:J1(N)-1) END IF 29 CONTINUE C DO 31 N=1,NUM1 IF(J1(N).LT.L1(N)) GO TO 20 31 CONTINUE DO 32 I=J2+1,L2 NAME2(I:I)=' ' 32 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE WRITTR(IRAY1,IRAY2,IRAY3) INTEGER IRAY1,IRAY2,IRAY3 C C This subroutine is designed to write the indices of rays forming C the homogeneous triangles in the ray-parameter domain. C C Input: C IRAY1,IRAY2,IRAY3... Indices of rays forming the homogeneous C triangle. IRAY3=0 in 2-D. C C No output. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 1997, October 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER IFORM1,IFORM3 CHARACTER*10 FORMAT SAVE IFORM1,IFORM3,FORMAT DATA IFORM1/99999/,IFORM3/0/,FORMAT/'(I6,I6,I2)'/ C C Writing to the output file with homogeneous triangles: IF(NAME2(2*(NF1+NF2+NF3)+1).NE.' ') THEN IF(IRAY1.GT.IFORM1.OR.IRAY2.GT.IFORM1.OR.IRAY3.GT.IFORM1) THEN IFORM1=IFORM1*10+9 FORMAT(3:3)=CHAR(ICHAR(FORMAT(3:3))+1) FORMAT(6:6)=FORMAT(3:3) IFORM3=0 END IF IF(IRAY3.GT.IFORM3) THEN IFORM3=IFORM1 FORMAT(9:9)=FORMAT(3:3) END IF WRITE(LUW(2*(NF1+NF2+NF3)+1),FORMAT) IRAY1,IRAY2,IRAY3 END IF C RETURN END C C======================================================================= C C C SUBROUTINE WRITBR(IB1,IB2,IB3,IB4) INTEGER IB1,IB2,IB3,IB4 C C This subroutine is designed to write the indices of boundary rays. C C Input: C IB1 ... Index of the boundary ray. C IB2 ... Index of the boundary ray coupled with IB1. C IB3,IB4 . Indices of two boundary rays equal in history with IB1, C coupled with rays equal in history with IB2, and connected C with boundary ray IB1 by sides of homogeneous triangles. C IB4=0 if there is only one ray of the described C properties. C IB3=0 and IB4=0 if there is no ray of the described C properties and for one-parametric shooting in 2-D. C C No output. C C Common blocks /WRITC/ and /WRITN/: INCLUDE 'writ.inc' C writ.inc C None of the storage locations of the common block are altered. C C No subroutines and external functions required. C C Date: 2001, October 29 C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER IFORM1,IFORM3 CHARACTER*13 FORMAT SAVE IFORM1,IFORM3,FORMAT DATA IFORM1/99999/,IFORM3/0/,FORMAT/'(I6,I6,I6,I2)'/ C C Writing to the output file with boundary rays: IF(NAME2(2*(NF1+NF2+NF3)+2).NE.' ') THEN IF(IB1.GT.IFORM1.OR.IB2.GT.IFORM1.OR.IB3.GT.IFORM1.OR. * IB4.GT.IFORM1) THEN IFORM1=IFORM1*10+9 FORMAT(3:3)=CHAR(ICHAR(FORMAT(3:3))+1) FORMAT(6:6)=FORMAT(3:3) FORMAT(9:9)=FORMAT(3:3) IFORM3=0 END IF IF(IB4.GT.IFORM3) THEN IFORM3=IFORM1 FORMAT(12:12)=FORMAT(3:3) END IF WRITE(LUW(2*(NF1+NF2+NF3)+2),FORMAT) IB1,IB2,IB3,IB4 END IF C RETURN END C C======================================================================= C C C SUBROUTINE WPTC2() C C This auxiliary subroutine to subroutine WRIT2 is called when starting C a new ray. It sets the number of points along a ray to IPT=0 C in common block /POINTC/. C C No input. C C No output. C C Common block /POINTC/: INCLUDE 'pointc.inc' C pointc.inc C C Date: 2013, February 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C IPT=0 RETURN END C C======================================================================= C C C SUBROUTINE WPTC31(JSRC,JWAVE,JRAY,NZ,JCB1,JSRF,KMAHJ,XJ,UEBJ,ZL,Z, * NQV,QV) INTEGER JSRC,JWAVE,JRAY,NZ,JCB1,JSRF,KMAHJ,NQV REAL XJ,UEBJ,ZL(6),Z(NZ),QV(9) C C This auxiliary subroutine to subroutine WRIT31 stores the quantities C at individual points along the ray into common block /POINTC/. C The quantities are then written into the unformatted output file by C subroutine APW2, or by subroutines APW3 and APW4. C C Input: C JSRC,JWAVE,JRAY,NZ,JCB1,JSRF,KMAHJ,XJ,UEBJ,ZL,Z,NQV,QV... C Equivalent to C ISRC,IWAVE,IRAY,NY,ICB1,ISRF,KMAH,X,UEBRAY,YL,Y,NPV,PV C in common block /POINTC/. C C No output. C C Common block /POINTC/: INCLUDE 'pointc.inc' C pointc.inc C C Date: 2013, February 21 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C C....................................................................... C C Sequential index of the point on a ray: IPT=IPT+1 C C Settings overlapping with WPTC4: ISRC =JSRC IWAVE =JWAVE IRAY =JRAY C C Settings unique for WPTC31 and WPTC32: NY =NZ ICB1 =JCB1 ISRF =JSRF KMAH =KMAHJ X =XJ UEBRAY=UEBJ DO 11 I=1,6 YL(I)=ZL(I) 11 CONTINUE DO 12 I=1,NZ Y(I) =Z(I) 12 CONTINUE NPV=NQV IF(NPV.GE.3) THEN DO 13 I=1,9 PV(I)=QV(I) 13 CONTINUE END IF C RETURN END C C======================================================================= C C C SUBROUTINE WPTC32(JSRC,JWAVE,JRAY,NZ,JCB1,JSRF,KMAHJ,XJ,UEBJ,ZL,Z, * ZC,NQV,QV) INTEGER JSRC,JWAVE,JRAY,NZ,JCB1,JSRF,KMAHJ,NQV REAL XJ,UEBJ,ZL(6),Z(27),ZC(NZ-27),QV(9) C C This auxiliary subroutine to subroutine WRIT32 stores the computed C quantities at specified surfaces into common block /POINTC/. C It is called by subroutine WRIT32 at the points of intersection of the C ray with specified surfaces. C The quantities are then written into the unformatted output file by C subroutine APW2, or by subroutines APW3 and APW4. C C Input: C JSRC,JWAVE,JRAY,NZ,JCB1,JSRF,KMAHJ,XJ,UEBJ,ZL,Z,NQV,QV... C Equivalent to C ISRC,IWAVE,IRAY,NY,ICB1,ISRF,KMAH,X,UEBRAY,YL,Y,NPV,PV C ZC(1:NZ-27)... Array equivalent to quantities Y(28:NZ) C in common block /POINTC/. C C No output. C C Common block /POINTC/: INCLUDE 'pointc.inc' C pointc.inc C C Date: 2013, February 21 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C C....................................................................... C C Settings overlapping with WPTC4: ISRC =JSRC IWAVE =JWAVE IRAY =JRAY C C Settings unique for WPTC31 and WPTC32: NY =NZ ICB1 =JCB1 ISRF =JSRF KMAH =KMAHJ X =XJ UEBRAY=UEBJ DO 11 I=1,6 YL(I)=ZL(I) 11 CONTINUE DO 12 I=1,27 Y(I) =Z(I) 12 CONTINUE DO 13 I=28,NZ Y(I) =ZC(I-27) 13 CONTINUE NPV=NQV IF(NPV.GE.3) THEN DO 14 I=1,9 PV(I)=QV(I) 14 CONTINUE END IF RETURN END C C======================================================================= C C C SUBROUTINE WPTC4(JSRC,JWAVE,JRAY,JCB1I,JEND,JSHEET,JREC,ZLI,ZI, * NQVI,QVI) INTEGER JSRC,JWAVE,JRAY,JCB1I,JEND,JSHEET,JREC,NQVI REAL ZLI(6),ZI(29),QVI(9) C C This auxiliary subroutine to subroutine WRIT4 stores the quantities at C the initial point of the ray into common block /POINTC/. C The quantities are then written into the unformatted output file by C subroutine APW1. C C Input: C JSRC,JWAVE,JRAY,JCB1I,JEND,JSHEET,JREC,ZLI,ZI,NQVI,QVI... C Equivalent to C ISRC,IWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI,NPVI,PVI C in common block /POINTC/. C C No output. C C Common block /POINTC/: INCLUDE 'pointc.inc' C pointc.inc C C Date: 2013, February 21 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C C....................................................................... C C Settings overlapping with WPTC31 and WPTC32: ISRC =JSRC IWAVE =JWAVE IRAY =JRAY C C Settings unique for WPTC4: ICB1I =JCB1I IEND =JEND ISHEET=JSHEET IREC =JREC DO 11 I=1,6 YLI(I)=ZLI(I) 11 CONTINUE DO 12 I=1,29 YI(I) =ZI(I) 12 CONTINUE NPVI=NQVI IF(NPVI.GE.3) THEN DO 13 I=1,9 PVI(I)=QVI(I) 13 CONTINUE END IF C RETURN END C C======================================================================= C