C Subroutine file 'writ.for' to create the output files of complete ray C tracing C C Date: 2001, January 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 (C.R.T.5.5.1). It is called with constant step STORE (see C input data in the file 'ray.for') of the independent C variable along the ray, and at the points of intersection C with interfaces either before and after the C transformation. C WRIT31 C WRIT32..Subroutine storing the computed quantities at specified C surfaces (C.R.T.5.5.2). It is called at the points of C intersection of the ray with the specified surfaces. C WRIT32 C WRIT33..Subroutine storing the computed quantities at the C endpoints of the individual elementary waves C (C.R.T.5.5.3). It is called at the endpoints of the C elements of rays, situated at structural interfaces. C WRIT33 C WRIT4...Subroutine storing the quantities at the initial point of C the ray (C.R.T.6.1). It is called after termination of C 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 WRITA...Subroutine writing the indices of rays forming the C homogeneous triangles in the ray-parameter domain. C WRITA 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 The data are read in by the list directed input (free format). In C the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). All C input variables are of the type CHARACTER. The data are read by C 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 significant. 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 Default: K2P=0. C (3) FN1A,FN1I,FN1T,/ C FN1A... String containing the first component of the name of the C output file with the computed quantities stored along the C rays (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 with the quantities at the initial points of C rays (C.R.T.6.1), related to the output file with the C computed quantities stored along the rays (C.R.T.5.5.1). C The other components of the filename are common with the C related file 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 specifying the homogeneous triangles C in the ray-parameter domain. The other components of the C filename are common with the related file containing the C computed quantities along rays. C Description of file with triangles C /... An obligatory slash. C Default: all filenames 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 with the computed quantities stored at the C specified surfaces (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 with the quantities at the initial points of C rays (C.R.T.6.1), related to the output file with the C computed quantities stored at the specified surfaces C (C.R.T.5.5.2). The other components of the filename are C common with the related file containing the computed C 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 with the computed quantities stored at the C I-th specified surface (C.R.T.5.5.2). Here NENDF is the C number of surfaces specified for storing the computed C quantities. C /... An obligatory slash. C Default: all filenames blank. C (5) Either (5.1) or (5.2): C (5.1) / C A slash is fully sufficient if MODCRT.LE.1, see input data C 'dcrt.dat' described in source code file 'ray.for'. 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 with the computed quantities stored at the C endpoints of the elements of rays (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 with the quantities at the initial points of C rays (C.R.T.6.1), related to the output file with the C computed quantities stored at the endpoints of the C elements of rays (C.R.T.5.5.3). The other components of C the filename are common with the related file containing C the computed 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 with the computed quantities stored at the C endpoints of the I-th element of rays (C.R.T.5.5.3). C These strings need not be specified for MODCRT.LE.1, see C input data 'dcrt.dat'. Otherwise, none of FN3B(i) C corresponding to the endpoint of a ray element at which C quantities are stored (see MODCRT in the input data C 'dcrt.dat') may equal spaces. C /... An obligatory slash. C Default: all filenames blank. C (6) For each elementary wave IWAVE the following data (5.1): C (6.1) FN1C,FN2C,FN3C,(FN2D(I),I=1,NENDF),/ C FN1C... String containing the second component of the name of the C output file with the computed quantities stored along the C rays (C.R.T.5.5.1). C If both the first and second component of the filename are C blank, the file is not created. C FN2C... String containing the third component of the name of the C output file with the computed quantities stored at the C specified surfaces (C.R.T.5.5.2). C If the second, third and fourth components of the filename C are all blank, the file is not created. C FN3C... String containing the third component of the name of the C output file with the computed quantities stored at the C endpoints of the elements of rays (C.R.T.5.5.3). C This string may not equal spaces for MODCRT.GE.2, see C input data 'dcrt.dat'. C If MODCRT.LE.1, see input data 'dcrt.dat', the file is not C created. C FN2D(I)... String containing the fourth component of the name of C the output file with the computed quantities stored at the C I-th specified surface (C.R.T.5.5.2). Here NENDF is the C number of surfaces specified for storing the computed C quantities. C /... An obligatory slash. C Default: all filenames blank. C Examples of data set WRIT: C writ.dat two-point rays stored. C writall.dat all rays stored. 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 consist 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'/ C 'CRT-R'.. Name of the unformatted CRT output file with the C quantities stored along rays (see C.R.T.5.5.1). C Description of file 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.R.T.5.5.2). C Description of file CRT-S C 'CRT-I'.. Name of the unformatted CRT output file with the C quantities at the initial points of rays (see C.R.T.6.1). C Description of file 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 CRT-T 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' C General defaults for the subsequent sets: C 'CRT-R'=' ', 'CRT-S'=' ', 'CRT-I'=' ', 'CRT-T'=' ' 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 rays C (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) IWAVE,IRAY,NY,ICB1,ISRF,X,(YL(I),I=1,6),(Y(I),I=1,NY) 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 (C.R.T.5.5.4). C Description of YL C Y... Array containing basic quantities computed along the ray 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 (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) 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.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 (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 (C.R.T.5.5.4). C Description of YL C Y... Array containing basic quantities computed along the ray. C (C.R.T.5.2.1). C Description of Y C YY... Array containing real auxiliary quantities computed along C the ray (C.R.T.5.2.2). C Description of YY C IY... Array containing integer auxiliary quantities computed C along the ray (C.R.T.5.2.2). 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 rays C (C.R.T.6.1): C Unformatted output. The quantities are stored at the initial C points of rays. C For whole rays (C.R.T.5.5.1), only initial points of two-point C rays are stored if only two-point rays are stored. C For intersections with surfaces (C.R.T.5.5.2) and endpoints of the C elements of rays (C.R.T.5.5.3) initial points of all rays C are stored. C For initial point - one record containing following 34 quantities: C (1) (-IWAVE),IRAY,ICB1I,IEND,ISHEET,IREC,(YLI(I),I=1,6),(YI(I),I=1,29) 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.R.T.6.1). C IEND... Reason of the termination of the computation of a ray (see C C.R.T.5.4). For a detailed description see subroutines C RAY2 (file 'ray.for') and INIT2 (file 'init.for'). C Description of IEND in RAY2 C Description of IEND in INIT2 C ISHEET..Ray-history index. The different ray histories are C consecutively indexed by positive integers 1,2,3,... C According to their appearance during ray tracing. C The ray histories are indexed independently within each C elementary wave. 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 RPAR4. C YLI... Array containing the values of the quantities YL(1)-YL(6), C see C.R.T.5.5.4, describing the local properties of the C model at the initial point of the ray. See the list of C YL(1) to YL(6) in the file 'ray.for'. 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.R.T.6.1. These quantities are C listed in the file 'init.for'. C Description of YI C C C Output files CRT-T specifying homogeneous triangles in the ray-parameter C domain: C Formatted output. C For each homogeneous triangle - one record containing 3 integers: C (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 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/ DATA JENDR/10,21,22,23,24,25,26,30,32,33,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: 2001, January 18 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=0 C C Reading the name of the input data 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=' ' READ(LUN,*) FN1A,FN1I,FN1T 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 CALL ERROR('553 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 ELSE C C Updating the source index IF(IWAVE.EQ.1) 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 CALL FNAME(2,NAME1,NAME2(1)) IF(NAME2(1).NE.' ') THEN 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 CALL FNAME(2,NAME1,NAME2(2)) IF(NAME2(2).NE.' ') THEN 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 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 NAME1(1)=FN2I IF(NAME1(2).EQ.' '.AND.NAME1(3).EQ.' '.AND.NAME1(4).EQ.' ') * THEN NAME1(1)=' ' END IF J=2*(NF1+I) 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 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 554 CALL ERROR('554 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 CALL FNAME(3,NAME1,NAME2(J)) IF(NAME2(J).NE.' ') THEN 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) CALL FNAME(3,NAME1,NAME2(J)) IF(NAME2(J).NE.' ') THEN 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 NAME2(J)=' ' IF(FN1T.NE.' ') THEN NAME1(1)=FN1T NAME1(2)=FN1C CALL FNAME(2,NAME1,NAME2(J)) C IF(NAME2(J).NE.' ') THEN 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 C 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 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 The type 1 to 7 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: 1996, June 15 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 Screen output CALL SCRO2(IRAY) RETURN END C C======================================================================= C C C SUBROUTINE WRIT31(YL,Y,YY,IY) REAL YL(6),Y(35),YY(5) INTEGER IY(12) C C This subroutine stores the computed quantities along a ray (see C C.R.T.5.5.1). It is called with constant step STORE of the C independent variable along the ray, and at the points of intersection C with interfaces 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: 2001, January 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C IF(NELEM(IY(2)).GT.KWRIT1) THEN IF(NAME2(1).NE.' ') THEN IF(K2P.EQ.0) THEN JPOINT(1)=JPOINT(1)+1 WRITE(LUW(1)) * ISRC,JWAVE,JRAY,IY(1),IY(5),IY(6),YY(1),YL,(Y(I),I=1,IY(1)) ELSE NTMP=NTMP+1 IF(NTMP.LE.MTMP) THEN KTMP(1,NTMP)=ISRC KTMP(2,NTMP)=JWAVE KTMP(3,NTMP)=JRAY KTMP(4,NTMP)=IY(1) KTMP(5,NTMP)=IY(5) KTMP(6,NTMP)=IY(6) TMP(1,NTMP)=YY(1) DO 11 I=1,6 TMP(1+I,NTMP)=YL(I) 11 CONTINUE DO 12 I=1,IY(1) TMP(7+I,NTMP)=Y(I) 12 CONTINUE ELSE C 552 CALL ERROR('552 in WRIT31: Too many points along a ray') C Dimension MTMP of array to temporarily store the traced C ray in the memory is smaller than the number of points C along the ray. END IF 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(12),NAMPL REAL YL(6),Y(27),YY(5),YC(NAMPL) C C This subroutine stores the computed quantities at specified surfaces C (see C.R.T.5.5.2). It is called at the points of intersection of the C ray with specified 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.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: 2001, January 18 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 WRITE(LUW(2*(NF1+ISTOR)-1)) * ISRC,JWAVE,JRAY,NY,IY(5),IY(6),YY(1),YL,Y,YC END IF END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE WRIT33(YL,Y,YY,IY) REAL YL(6),Y(35),YY(5) INTEGER IY(12) C C This subroutine stores the computed quantities at the endpoints of the C individual elementary waves (see C.R.T.5.5.3). It is called at the C endpoints of the elements of rays, situated 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: 1996, June 15 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage location: INTEGER 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 WRITE(LUW(2*(NF1+NF2+I)-1)) JRAY,IY,YL,(Y(I),I=1,IY(1)),YY 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(12),IEND,ISHEET,IREC REAL YL(6),Y(35),YY(5) C C This subroutine stores the quantities at the initial point of the ray C (SEE C.R.T.6.1). 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). For a detailed description see subroutines C RAY2 (file 'ray.for') and INIT2 (file 'init.for'). 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: 2001, January 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,I0,I1,I2,I3,I4,I5,I6,I7,I8,I9 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 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 at the initial point of the ray (C.R.T.6.1) I1=-JWAVE 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 WRITE(LUW(2*I)) ISRC,I1,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI END IF END IF 20 CONTINUE 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 DO 30 J=1,NTMP WRITE(LUW(1)) (KTMP(I,J),I=1,6),(TMP(I,J),I=1,7+KTMP(4,J)) 30 CONTINUE END IF END IF C C Writing to the output LOG file C Simplecticity of the propagator matrix (see 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 WRITA(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