C
C Subroutine file 'fdaux.for' to mediate input and output of finite C diference programs C C Version: 5.40 C Date: 1999, June 28 C Coded by Ludek Klimes C C....................................................................... C C This Fortran77 file consists of the following external procedures: C FDIN... Subroutine to read the effective elastic parameters for C 2-D finite difference programs. C FDIN C FDIN1...Subroutine to read the values of a single effective C elastic parameter for 2-D finite difference programs. C FDIN1 C FDOUT1..Subroutine to write the discrete values of a single C quantity for 2-D finite difference programs. C FDOUT1 C FDGRID..Subroutine to read grid dimensions for 2-D finite C differences. C FDGRID C FDREC...Subroutine to read receiver and source positions for 2-D C finite differences. C FDREC C C======================================================================= C C C SUBROUTINE FDIN(LU,NX,NZ,AHL,BHL,CHL,AVL,BVL,CVL, * AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC,WORK) C INTEGER LU,NX,NZ REAL AHL(NZ ,NX-1),BHL(NZ ,NX-1),CHL(NZ ,NX-1) REAL AVL(NZ-1,NX ),BVL(NZ-1,NX ),CVL(NZ-1,NX ) REAL AHS(NZ-1,NX-1),BHS(NZ-1,NX-1),CHS(NZ-1,NX-1) REAL AVS(NZ-1,NX-1),BVS(NZ-1,NX-1),CVS(NZ-1,NX-1) REAL DEN(NZ ,NX ),QFC(NZ ,NX ),WORK(NZ ,NX ) C C Subroutine FDIN reads the effective elastic parameters for 2-D finite C difference programs. C C Input: C LU... Logical unit number to be used to read all files. C NX... Number of grid points in the horizontal direction. C NZ... Number of grid points in the vertical direction. C AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC,WORK... C Arrays of the dimensions as declared above or greater. C WORK is a temporary working array, its content will be C destroyed. C C Output: C AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS,CVS,DEN,QFC... Values C read from formatted or unformatted files given by SEP C parameters AHL,BHL,CHL,AVL,BVL,CVL,AHS,BHS,CHS,AVS,BVS, C CVS,DEN,QFC, respectively. C WORK... Undefined. C C Date: 1998, November 11 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER NEWPAR C NEWPAR..NEWPAR=0: Average parameters along the legs are read. C NEWPAR=1: Average parameters over the grid facets are read C and the parameters corresponding to the legs are derived C from them. C C....................................................................... C CALL RSEP3I('NEWPAR',NEWPAR,0) IF(NEWPAR.EQ.0) THEN CALL FDIN1(LU,NX-1,NZ-1,'A13',AHS,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ-1,'B13',BHS,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ-1,'C13',CHS,WORK,0,CHS) CALL FDIN1(LU,NX-1,NZ-1,'A31',AVS,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ-1,'B31',BVS,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ-1,'C31',CVS,WORK,0,CHS) CALL FDIN1(LU,NX-1,NZ ,'A11',AHL,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ ,'B11',BHL,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ ,'C11',CHL,WORK,0,CHS) CALL FDIN1(LU,NX ,NZ-1,'A33',AVL,WORK,0,AVS) CALL FDIN1(LU,NX ,NZ-1,'B33',BVL,WORK,0,BVS) CALL FDIN1(LU,NX ,NZ-1,'C33',CVL,WORK,0,CVS) ELSE CALL FDIN1(LU,NX-1,NZ-1,'AA13',AHS,WORK,0,AHS) CALL FDIN1(LU,NX-1,NZ-1,'BB13',BHS,WORK,0,BHS) CALL FDIN1(LU,NX-1,NZ-1,'CC13',CHS,WORK,0,CHS) CALL FDIN1(LU,NX-1,NZ-1,' ',AVS,WORK,4,AHS) CALL FDIN1(LU,NX-1,NZ-1,' ',BVS,WORK,4,BHS) CALL FDIN1(LU,NX-1,NZ-1,' ',CVS,WORK,4,CHS) CALL FDIN1(LU,NX-1,NZ ,' ',AHL,WORK,2,AHS) CALL FDIN1(LU,NX-1,NZ ,' ',BHL,WORK,2,BHS) CALL FDIN1(LU,NX-1,NZ ,' ',CHL,WORK,2,CHS) CALL FDIN1(LU,NX ,NZ-1,' ',AVL,WORK,1,AVS) CALL FDIN1(LU,NX ,NZ-1,' ',BVL,WORK,1,BVS) CALL FDIN1(LU,NX ,NZ-1,' ',CVL,WORK,1,CVS) END IF CALL FDIN1(LU,NX ,NZ ,'DEN',DEN,WORK,0,DEN) CALL FDIN1(LU,NX ,NZ ,'QFC',QFC,WORK,0,QFC) RETURN END C C======================================================================= C C C SUBROUTINE FDIN1(LU,N1,N2,NAME,ARRAY,WORK,KDEF,DEF) C INTEGER LU,N1,N2,KDEF CHARACTER*(*) NAME REAL ARRAY(N2*N1),WORK(N1*N2),DEF(*) C C Subroutine FDIN1 reads the values of a single effective elastic C parameter for 2-D finite difference programs. C C Input: C LU... Logical unit number to be used for input. C N1... Number of values in the horizontal direction. C N2... Number of values in the vertical direction. C NAME... String containing the name of the parameter specifying the C name of the file containing the grid values to be read. C Except for its case, it should match the parameter name in C a previously read input SEP parameter file. If parameter C NAME is not defined in the input SEP parameter file, no C input is performed by this subroutine and the output C arrays are filled with zeros. C KDEF... Kind of default values: C KDEF=0: Zeros. Array DEF is not used. C KDEF=1: Array DEF should have dimensions DEF(N2,N1-1). C Then the default is ARRAY(I2,1)=DEF(I2,1), C ARRAY(I2,N1)=DEF(I2,N1-1), and for 1.LT.I1.LT.N1 C ARRAY(I2,I1)=(DEF(I2,I1-1)+DEF(I2,I1))/2. C KDEF=2: Array DEF should have dimensions DEF(N2-1,N1). C Then the default is ARRAY(1,I1)=0.5*DEF(1,I1), C ARRAY(N2,I1)=DEF(N2-1,I1), and for 1.LT.I2.LT.N2 C ARRAY(I2,I1)=(DEF(I2-1,I1)+DEF(I2,I1))/2. C KDEF=4: Array DEF should have dimensions DEF(N2,N1). C Then the default is ARRAY(I2,I1)=DEF(I2,I1). C DEF... Array to specify default values. C C Output: C ARRAY...Array of N2*N1 values. The values are read from the file C which name is given by parameter NAME in SEP input data C read previously to the invocation of this subroutine. C Since the inner loop in the input file is over the C horizontal grid lines (N1 values), whereas the inner loop C in 2-D FD programs is over the vertical grid lines C (N2 values), N1*N2 matrix WORK read from the file is C transposed to yield N2*N1 output array ARRAY. C WORK... Temporary storage array containing the grid values in the C same order as in the input file. C C Date: 1998, November 11 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*80 FILE CHARACTER*11 FORM INTEGER I1,I2,I,J C C....................................................................... C C Name and form of the input data file: IF(NAME.EQ.' ') THEN FILE=' ' ELSE CALL RSEP3T(NAME,FILE,' ') CALL RSEP3T('FORM',FORM,'FORMATTED') END IF C C Check whether the filename is specified in the input SEP data: IF(FILE.EQ.' ') THEN C C No input file specified: I=1 J=1 IF(KDEF.EQ.1) THEN C Averaging along the X1 axis: DO 12 I2=1,N2 ARRAY(I)=DEF(J) DO 11 I1=2,N1-1 I=I+N2 J=J+N2 ARRAY(I)=0.5*(DEF(J-N2)+DEF(J)) 11 CONTINUE I=I+N2 ARRAY(I)=DEF(J) I=I-N2*(N1-1)+1 J=J-N2*(N1-2)+1 12 CONTINUE ELSE IF(KDEF.EQ.2) THEN C Averaging along the X2 axis: DO 22 I1=1,N1 ARRAY(I)=0.5*DEF(J) DO 21 I2=2,N2-1 I=I+1 J=J+1 ARRAY(I)=0.5*(DEF(J-1)+DEF(J)) 21 CONTINUE I=I+1 ARRAY(I)=DEF(J) I=I+1 J=J+1 22 CONTINUE ELSE IF(KDEF.EQ.4) THEN C Copying: DO 40 I=1,N2*N1 ARRAY(I)=DEF(I) 40 CONTINUE ELSE C Zeros: DO 90 I=1,N1*N2 ARRAY(I)=0. 90 CONTINUE END IF C C Transposing the array: * CALL GMTRA(ARRAY,WORK,N2,N1) C ELSE C C Reading the array (undefined values are replaced by zeros): CALL RARRAY(LU,FILE,FORM,.TRUE.,0.,N1*N2,WORK) C C Transposing the array: CALL GMTRA(WORK,ARRAY,N1,N2) C END IF RETURN END C C======================================================================= C C C SUBROUTINE FDOUT1(LU,N1,N2,SIGNUM,NAME,ARRAY,WORK) C INTEGER LU,N1,N2 CHARACTER*(*) NAME REAL SIGNUM,ARRAY(N2*N1),WORK(N1*N2) C C Subroutine FDOUT1 to write the discrete values of a single quantity C for 2-D finite difference programs. C C Input: C LU... Logical unit number to be used for output. C N1... Number of values in the horizontal direction. C N2... Number of values in the vertical direction. C SIGNUM..If SIGNUM is negative, the array values will be written C with opposite signs. C NAME... String containing the name of the parameter specifying the C name of the output file to be written. C Except for its case, it should match the parameter name in C a previously read input SEP parameter file. If parameter C NAME is not defined in the input SEP parameter file, no C output is performed. C ARRAY...Array of N2*N1 values. The values will be written to the C file which name is given by parameter NAME in SEP input C data read previously to the invocation of this subroutine. C Since the inner loop in the input file is over the C horizontal grid lines (N1 values), whereas the inner loop C in 2-D FD programs is over the vertical grid lines C (N2 values), N2*N1 matrix ARRAY is transposed to yield C N1*N2 array WORK written to the output file. C WORK... Temporary working array, its content will be destroyed. C C Output: C WORK... Undefined. C C Date: 1999, April 6 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*80 FILE CHARACTER*11 FORM INTEGER I C C....................................................................... C C Name of the output data file: IF(NAME.EQ.' ') THEN C The output file is already open: FILE=' ' ELSE CALL RSEP3T(NAME,FILE,' ') C Check whether the filename is specified in the input SEP data: IF(FILE.EQ.' ') THEN RETURN END IF END IF C C Form of the output data file: CALL RSEP3T('FORM',FORM,'FORMATTED') C C Transposing the array: CALL GMTRA(ARRAY,WORK,N2,N1) C IF(SIGNUM.LT.0) THEN DO 10 I=1,N1*N2 WORK(I)=-WORK(I) 10 CONTINUE END IF C C Writing the array: CALL WARRAY(LU,FILE,FORM,.FALSE.,0.,.FALSE.,0.,N1*N2,WORK) C RETURN END C C======================================================================= C C C SUBROUTINE FDGRID(KMAX,LMAX,KLMAX,NX,NZ,DX) C INTEGER KMAX,LMAX,KLMAX,NX,NZ REAL DX C C Subroutine to read grid dimensions for 2-D finite differences. C C Input: C KMAX... Maximum number of vertical gridlines. C LMAX... Maximum number of horizontal gridlines. C KLMAX...Maximum number of gridpoints. C C Output: C NX... Number of vertical gridlines. C NZ... Number of horizontal gridlines. C DX... Absolute value of the grid interval. The absolute values C of the horizontal and vertical grid intervals must be C equal. C C Date: 1998, September 11 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER N1,N2,N3 REAL D1,D2,D3,DZ C C....................................................................... C C Recalling the data specifying grid dimensions C (arguments: Name of parameter in input data, Variable, Default): CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) IF (N1.LE.1) THEN NX=N2 NZ=N3 DX=ABS(D2) DZ=ABS(D3) ELSE IF (N2.LE.1) THEN NX=N1 NZ=N3 DX=ABS(D1) DZ=ABS(D3) ELSE IF (N3.LE.1) THEN NX=N1 NZ=N2 DX=ABS(D1) DZ=ABS(D2) ELSE C FDAUX-01 CALL ERROR('FDAUX-01: 3-D grid specified for 2-D FD') C Exactly one of input parameters N1, N2 ,N3 must equal 1 for 2-D C finite differences. END IF IF (NX.LE.1.OR.NZ.LE.1) THEN C FDAUX-02 CALL ERROR('FDAUX-02: 1-D grid specified for 2-D FD') C Exactly one of input parameters N1, N2 ,N3 must equal 1 for 2-D C finite differences. END IF IF (NX.GT.KMAX) THEN C FDAUX-03 CALL ERROR('FDAUX-03: Small dimension KMAX') C Number NX of vertical gridlines exceeds maximum number KMAX of C vertical gridlines. END IF IF (NZ.GT.LMAX) THEN C FDAUX-04 CALL ERROR('FDAUX-04: Small dimension LMAX') C Number NZ of horizontal gridlines exceeds maximum number LMAX of C horizontal gridlines. END IF IF (NX*NZ.GT.KLMAX) THEN C FDAUX-05 CALL ERROR('FDAUX-05: Small dimension KLMAX') C Number N1*N2*N3 of gridpoints exceeds maximum number KLMAX of C gridpoints. END IF IF (DX.NE.DZ) THEN C FDAUX-06 CALL ERROR('FDAUX-06: Different grid intervals') C Horizontal and vertical grid intervals must be equal for C 2-D finite-difference programs. END IF C RETURN END C C======================================================================= C C C SUBROUTINE FDREC(LU,MAXREC,NREC,REC,KRE,LRE,KSOUR,LSOUR) C INTEGER LU,MAXREC,NREC,KRE(MAXREC),LRE(MAXREC),KSOUR,LSOUR CHARACTER*(*) REC(MAXREC) C C Subroutine to read receiver and source positions for 2-D finite C differences. C C Input: C LU... Logical unit number to be used for input. C MAXREC..Maximum number of receivers. C C Output: C NREC... Number of receivers. C KREC,LREC... Horizontal and vertical indices of receiver C gridpoints. C KSOUR,LSOUR... Horizontal and vertical indices of the source C gridpoint. C C Date: 1998, October 7 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C REAL UNDEF PARAMETER (UNDEF=-999999.) C C Auxiliary storage locations: CHARACTER*80 FREC CHARACTER*6 TEXT INTEGER N1,N3 REAL D1,D3,O1,O3,X1,X2,X3 C C....................................................................... C CALL RSEP3I('N1' ,N1 ,1 ) CALL RSEP3I('N3' ,N3 ,1 ) CALL RSEP3R('D1' ,D1 ,1. ) CALL RSEP3R('D3' ,D3 ,1. ) CALL RSEP3R('O1' ,O1 ,0. ) CALL RSEP3R('O3' ,O3 ,0. ) C NREC=0 CALL RSEP3T('REC',FREC,' ') IF(FREC.NE.' ') THEN OPEN(LU,FILE=FREC,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) 11 CONTINUE X1=UNDEF X2=UNDEF X3=UNDEF READ(LU,*,END=12) TEXT,X1,X2,X3 IF(X1.EQ.UNDEF.AND.X2.EQ.UNDEF.AND.X3.EQ.UNDEF) THEN GO TO 12 END IF IF(NREC.GE.MAXREC) THEN C FDAUX-07 CALL ERROR('FDAUX-07: Small dimension MAXREC') C Dimension MRAM must be at least C 19*N1*N3+10*N1+12*N3+(2*NTFD+5)*MAX(1,NREC), C where N1, N3 and NTFD are parameters in the main input file C and NREC is the number of receivers specified in the file C which name is given by parameter REC. END IF NREC=NREC+1 REC(NREC)=TEXT IF(X1.EQ.UNDEF) X1=0. IF(X2.EQ.UNDEF) X2=0. IF(X3.EQ.UNDEF) X3=0. KRE(NREC)=NINT((X1-O1)/D1)+1 LRE(NREC)=NINT((X3-O3)/D3)+1 IF(KRE(NREC).LT.1.OR.KRE(NREC).GT.N1.OR. * LRE(NREC).LT.1.OR.LRE(NREC).GT.N3) THEN C FDAUX-08 CALL ERROR('FDAUX-08: Receiver point outside the FD grid') END IF GO TO 11 12 CONTINUE CLOSE(LU) END IF C CALL RSEP3T('SRC',FREC,' ') IF(FREC.NE.' ') THEN OPEN(LU,FILE=FREC,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) X1=0. X2=0. X3=0. READ(LU,*) TEXT,X1,X2,X3 CLOSE(LU) KSOUR=NINT((X1-O1)/D1)+1 LSOUR=NINT((X3-O3)/D3)+1 IF(KSOUR.LT.1.OR.KSOUR.GT.N1.OR. * LSOUR.LT.1.OR.LSOUR.GT.N3) THEN C FDAUX-09 CALL ERROR('FDAUX-09: Source point outside the FD grid') END IF END IF RETURN END C C======================================================================= C