C
C Program FDMOD to generate finite-difference effective parameters C (harmonic averages of elastic parameters) in 2-D models specified in C terms of homogeneous convex polygons C C Version: 5.30 C Date: 1999, April 2 C C Program written by Jiri Zahradnik C Department of Geophysics, Charles University Prague C V Holesovickach 2, 180 00 Praha 8, Czech Republic C Tel.: +420-2-21912546, Fax: +420-2-21912555 C E-mail: jz@karel.troja.mff.cuni.cz C Data input and output modified by Ludek Klimes C C....................................................................... C C Contents: C Description of the polygonal model C Description of data files C List of external procedures C C For the description of effective elastic parameters refer to file C modfd.for. C C....................................................................... C C Description of the polygonal model: C C The medium is composed of homogeneous polygonal blocks. The C blocks do not overlap each other and cover the computational region. C Convex blocks only are allowed (internal angles.le.180 degrees). C For this reason even some physically homogeneous parts C of the model must be subdivided in several (formal) blocks. C Numbering of the blocks is arbitrary. C Each block is characterized by the physical coordinates of its C nodes (in meters), S-wave velocity (m/s), P/S velocity ratio, C and density (kg/m**3). The nodes are numbered counterclockwise. C Also specified for each block segment is the serial number of the C neighbouring block. C See the accompanying example, or ask the author of the program C for more details. C C Description and remarks on C FD method C C....................................................................... C C Description of data files: C C Main input data read from the * input device: C The data are read in by the list directed input (free format), by C a single READ statement. The data consist of character strings C which have to be enclosed in apostrophes. Note that the C interactive * external unit may be piped or redirected to a file. C The following items are read: C (1) 'FDSEP' / C 'FDSEP'... String in apostrophes containing the name of the main C input data file with FD parameters in the SEP format. C The file specifies the filename of the polygonal model, C dimensions of the grid of N1*N2*N3 gridpoints, the form C and names of the files with effective elastic parameters, C etc. C /... It is recommended to terminate these data by a slash. C Default: 'FDSEP'='fd.h'. C C Data file 'FDSEP' has the form of the SEP (Stanford Exploration C Project) parameter file: C All the data are specified in the form of PARAMETER=VALUE, e.g. C N1=50, with PARAMETER directly preceding = without intervening C spaces and with VALUE directly following = without intervening C spaces. The PARAMETER=VALUE couple must be delimited by a space C or comma from both sides. C The PARAMETER string is not case-sensitive. C PARAMETER= followed by a space resets the default parameter value. C All other text in the input files is ignored. The file thus may C contain unused data or comments without leading comment character. C Everything between comment character # and the end of the C respective line is ignored, too. C The PARAMETER=VALUE couples may be specified in any order. C The last appearance takes precedence. C Data specifying the model by means of the MODEL package: C FDMOD='string'... Name of the input data file describing a 2-D C model composed of homogeneous convex polygons. C The data file contains the name of the model, the number C of blocks, and for each block: C number of edges of the block, shear-wave velocity C (e.g., in m/s) inside the block, P/S velocity ratio inside C the block, density (e.g., in kg/m**3) inside the block, C x and z coordinates (e.g., in m) of the initial point of C segment, sequential number of the neighbouring block C (outside put 0). C Parameter FDMOD is required by this program (but no other C one). C Example of data set C FDMOD C Default: FDMOD=' ' C This program understands the following parameters common to nearly all C programs dealing with SEP-like gridded data formats: C N1=positive integer... Number of gridpoints of the basic grid C along the X1 axis. Default: N1=1 C N3=positive integer... Number of gridpoints of the basic grid C along the X3 axis. Default: N3=1 C O1=real... X1 coordinate of the origin of the grid. Default: O1=0. C O3=real... X3 coordinate of the origin of the grid. C Note that FD program 'fd2d.for' assumes the free Earth C surface at vertical coordinate X3=O3. C Default: O3=0. C D1=real... Grid spacing along the X1 axis. Default: D1=1. C D3=real... Grid spacing along the X3 axis. Default: D3=1. C The grid intervals may also be negative. C And the following parameters specific to effective-parameter FD codes: C FORM='string'... Form of the output files: C FORM='formatted': Formatted ASCII output files, C FORM='unformatted': Unformatted output files. C A11='string', B11='string', C11='string'... Strings in apostrophes C containing the names of the output files with (N1-1)*N2*N3 C values of the elastic parameters averaged in the direction C the X1 axis, along the gridlines. C Defaults: A11=' ', B11=' ', C11=' '. C A33='string', B33='string', C33='string'... Strings in apostrophes C containing the names of the output files with N1*N2*(N3-1) C values of the elastic parameters averaged in the direction C the X3 axis, along the gridlines. C Defaults: A33=' ', B33=' ', C33=' '. C A13='string', B13='string', C13='string'... Strings in apostrophes C containing the names of the output files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X1 axis, along the gridlines shifted C half a grid interval in the direction of the X3 axis. C Defaults: A13=' ', B13=' ', C13=' '. C A31='string', B31='string', C31='string'... Strings in apostrophes C containing the names of the output files with C (N1-1)*N2*(N3-1) values of the elastic parameters averaged C in the direction the X3 axis, along the gridlines shifted C half a grid interval in the direction of the X1 axis. C Defaults: A31=' ', B31=' ', C31=' '. C DEN='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the C density. C Default: DEN=' '. C QFC='string'... String in apostrophes containing the name of the C output ASCII file with N1*N2*N3 grid values of the P-wave C quality factor. C Default: QFC=' ' C Example of data set FDSEP C C....................................................................... C C This file consists of the following external procedures: C FDMOD...Main program allocating the memory for FDPAR. C FDMOD C FDPAR...Subroutine to generate the parameters (originally the main C program). C FDPAR C XLEG... Subroutine for computation of effective (mean) squared C velocities for a grid leg. C XLEG C GRPT... Subroutine for location of grid point (XPT,ZPT) with C respect to block IBL. C GRPT C TWOPT...Subroutine for location of a grid leg with respect to C boundary of block IB. C TWOPT C C======================================================================= C C C PROGRAM FDMOD C C Main program allocating the memory. C C Subroutines and external functions required: EXTERNAL FDPAR,RSEP1,RSEP3I C FDPAR...This file. C RSEP1,RSEP3I... File 'sep.for'. C Note that the above subroutines reference many other external C procedures from various subroutine files. These indirectly C referenced procedures are not named here, but are listed in the C particular subroutine files. C C----------------------------------------------------------------------- C C Common block /RAMC/ to allocate large array RAM: INCLUDE 'ram.inc' C ram.inc C C----------------------------------------------------------------------- C C Input data filenames and logical unit numbers: INTEGER LU PARAMETER (LU=1) CHARACTER*80 FSEP C LU... Logical unit number to be used to read and write all data C files. C FSEP... The name of the input data file with FD parameters in SEP C format. C INTEGER N1,N3,N,NR,NS,NGRID,I1,I2,I3,I4 C C....................................................................... C C Main input data: C Default: FSEP='fd.h' C Reading main input data: WRITE(*,'(A)') ' FDMOD: Enter name of input SEP file [''fd.h'']: ' READ (*,*) FSEP WRITE(*,'(A)') '+FDMOD: Working... ' C C Reading all the data from input SEP parameter file: CALL RSEP1(LU1,FSEP) C C Data specifying numbers of gridpoints: CALL RSEP3I('N1' ,N1 ,1) CALL RSEP3I('N2' ,N3 ,1) C C Allocating the memory in array RAM: NGRID=15*N1*N3 IF(NGRID.GT.MRAM) THEN C FDMOD-01 CALL ERROR('FDMOD-01: Small dimension MRAM in ram.inc') C Dimension MRAM must be at least C 15*N1*N3 , C where N1 and N3 are parameters in the main input file. END IF N =N1*N3 C CALL FDPAR(LU,N, * RAM( 1),RAM( N + 1),RAM( 2*N + 1),RAM( 3*N + 1), * RAM( 4*N + 1),RAM( 5*N + 1),RAM( 6*N + 1),RAM( 7*N + 1), * RAM( 8*N + 1),RAM( 9*N + 1),RAM(10*N + 1),RAM(11*N + 1), * RAM(12*N + 1),RAM(13*N + 1),RAM(14*N + 1)) C WRITE(*,'(A)') '+FDMOD: Done. ' STOP END C C======================================================================= C C C SUBROUTINE FDPAR(LU,KLMAX, * AHS,BHS,CHS,AVS,BVS,CVS, * AHL,BHL,CHL,AVL,BVL,CVL,DEN,QFC,WORK) INTEGER LU,KLMAX REAL AHS(KLMAX),BHS(KLMAX),CHS(KLMAX) REAL AVS(KLMAX),BVS(KLMAX),CVS(KLMAX) REAL AHL(KLMAX),BHL(KLMAX),CHL(KLMAX) REAL AVL(KLMAX),BVL(KLMAX),CVL(KLMAX) REAL DEN(KLMAX),QFC(KLMAX),WORK(KLMAX) C C Subroutine (originally main program) to read the input data and C generate finite-difference effective parameters. C C Input: C LU... Logical unit number to be used to read and write data C files. C KLMAX...Array dimension. C AHL,BHL,CHL, C AVL,BVL,CVL, C AHS,BHS,CHS, C AVS,BVS,CVS, C DEN,QFC...arrays with effective parameters C WORK... temporary working array C C No output. C C Subroutines and external functions required: EXTERNAL XLEG,GRPT,TWOPT EXTERNAL RSEP3T,FDGRID,FDOUT1 C XLEG... This file. C GRPT... This file. C TWOPT.. This file. C RSEP3T..File 'sep.for'. C FDGRID..File 'fdaux.for'. C FDOUT1..File 'fdaux.for'. C Note that the above subroutines reference many other external C procedures from various subroutine files. These indirectly C referenced procedures are not named here, but are listed in the C particular subroutine files. C C C----------------------------------------------------------------------- C C Storage locations: C C Maximum array dimensions: PARAMETER (MAXBL=50,MAXEDG=10) C MAXBL.. Max number of polygonal blocks. C MAXEDG..Max number of edges for any block. C C Arrays: DIMENSION NOED(0:MAXBL),VS(0:MAXBL),PSRAT(0:MAXBL),HU(0:MAXBL) DIMENSION XED(MAXEDG,MAXBL),ZED(MAXEDG,MAXBL),NOVIC(MAXEDG,MAXBL) C COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C NOBL... Number of blocks. C NOED... Number of edges of the block. C VS... Shear-wave velocity (m/s) inside the block. C PSRAT.. P/S velocity ratio inside the block. C HU... Density (kg/m**3) inside the block. C XED,ZED. X and Z coordinates (m) of the initial point of segment. C NOVIC... Serial number of the neigh. block (outside put 0). C CHARACTER*80 FMOD CHARACTER *5 NAMDUM C FMOD... The name of the input data file with model parameters in C list directed (free) format. C NAMDUM... Auxiliary storage location. C C....................................................................... C C Data specifying grid dimensions CALL FDGRID(KLMAX,KLMAX,KLMAX,NX,NZ,DX) C C Reading data describing a model composed of homogeneous polygons: CALL RSEP3T('FDMOD',FMOD,' ') OPEN(LU,FILE=FMOD,STATUS='OLD') C READ(1,'(1X,A5)') NAMDUM C READ(1,'(25X,I5)') NSYM NSYM=0 C READ(1,'(25X,I5)') NOBL DO 150 I=1,NOBL READ(1,'(25X,I5,3F10.0)') NOED(I),VS(I),PSRAT(I),HU(I) DO 100 J=1,NOED(I) READ(1,'(25X,2F10.0,I5)') XED(J,I),ZED(J,I),NOVIC(J,I) 100 CONTINUE 150 CONTINUE C WRITE(*,'(A)') '+FDMOD: Computing effective parameters ' C C Auxiliary parameters. LACC=0 C DO 310 I=1,NOBL C Attention: Now VS no more contains shear velocity but rigidity C and PSRAT is squared Vp/Vs ratio. The corresponding C change was made in FLERIB,FBOB,FCORNB. C VS(I)=VS(I)*VS(I) VS(I)=VS(I)*VS(I)*HU(I) PSRAT(I)=PSRAT(I)*PSRAT(I) 310 CONTINUE C NOED(0)=0 VS(0)=0. PSRAT(0)=1. C DXA=DX C C Vertical and horizontal parameters of medium for INTERlegs C (i.e. for legs between grid lines) NPUP=(NX-1)*(NZ-1) DO 430 IPU=1,NPUP NVER=NZ-1 NP=IFIX((IPU-0.01)/NVER) XSOUR=DX/2.+FLOAT(NP)*DX NPOR=IPU-NP*NVER ZSOUR=FLOAT(NPOR-1)*DX CALL XLEG(0,XSOUR,ZSOUR,DXA,AL,BE,RH) AVS(IPU)=AL BVS(IPU)=BE CVS(IPU)=AL-2.*BE NHOR=NZ-1 NP=IFIX((IPU-0.01)/NHOR) XSOUR=FLOAT(NP)*DX NPOR=IPU-NP*NHOR ZSOUR=DX/2.+FLOAT(NPOR-1)*DX CALL XLEG(1,XSOUR,ZSOUR,DXA,AL,BE,RH) AHS(IPU)=AL BHS(IPU)=BE CHS(IPU)=AL-2.*BE 430 CONTINUE C C Vertical parameters of medium for LARGE squares (= grid legs) C (for ignoring detailed computations, choose LACC=0) DO 450 K=2,NX-1 XOLD=FLOAT(K-1)*DX LAU=(K-1)*NZ XNZ=FLOAT(NZ) DO 400 L=1,NZ-1 ICEN=L+LAU JCEN=ICEN-IFIX(FLOAT(ICEN)/XNZ) ZOLD=FLOAT(L-1)*DX IF(L.LT.LACC) THEN CALL XLEG(0,XOLD,ZOLD,DXA/2.,AL1,BE1,RH1) ZNEW=ZOLD+DX/2. CALL XLEG(0,XOLD,ZNEW,DXA/2.,AL2,BE2,RH2) AL=(AL1+AL2)/2. BE=(BE1+BE2)/2. ELSE CALL XLEG(0,XOLD,ZOLD,DXA,AL,BE,RH) ENDIF AVL(JCEN)=AL BVL(JCEN)=BE CVL(JCEN)=AL-2.*BE 400 CONTINUE 450 CONTINUE C C Horizontal parameters of medium for LARGE squares (=grid legs) C (incl. density) DO 550 K=1,NX-1 XOLD=FLOAT(K-1)*DX LAU=(K-1)*NZ DO 500 L=1,NZ-1 ICEN=L+LAU ZOLD=FLOAT(L-1)*DX IF(L.LE.LACC) THEN CALL XLEG(1,XOLD,ZOLD,DXA/2.,AL1,BE1,RH1) XNEW=XOLD+DX/2. CALL XLEG(1,XNEW,ZOLD,DXA/2.,AL2,BE2,RH2) AL=(AL1+AL2)/2. BE=(BE1+BE2)/2. ELSE CALL XLEG(1,XOLD,ZOLD,DXA,AL,BE,RH) ENDIF AHL(ICEN)=AL BHL(ICEN)=BE CHL(ICEN)=AL-2.*BE DEN(ICEN)=RH 500 CONTINUE 550 CONTINUE C C Vertical and horizontal parameters for legs in the strip between C the inner and outer nonreflecting boundaries (a continuation) C (incl. density for horizontal legs) DO 600 L=1,NZ-1 ICEN=L ICEN2=ICEN+NZ JCEN=L JCEN2=JCEN+NZ-1 AVL(JCEN)=AVL(JCEN2) BVL(JCEN)=BVL(JCEN2) CVL(JCEN)=CVL(JCEN2) 600 CONTINUE LAU=(NX-1)*NZ DO 610 L=1,NZ-1 ICEN=L+LAU JCEN=ICEN-IFIX(FLOAT(ICEN)/FLOAT(NZ)) JCEN2=JCEN-NZ+1 AVL(JCEN)=AVL(JCEN2) BVL(JCEN)=BVL(JCEN2) CVL(JCEN)=CVL(JCEN2) 610 CONTINUE DO 620 K=1,NX-1 ICEN=K*NZ ICEN2=ICEN-1 AHL(ICEN)=AHL(ICEN2) BHL(ICEN)=BHL(ICEN2) CHL(ICEN)=CHL(ICEN2) DEN(ICEN)=DEN(ICEN2) 620 CONTINUE C C Halving the density at the free surface IZ=1: DO 70 I=1,NX*NZ,NZ DEN(I)=0.5*DEN(I) 70 CONTINUE C C No attenuation is considered in this version DO 80 I=1,NX*NZ QFC(I)=0. 80 CONTINUE C C Writing the gridded effective elastic parameters: CALL FDOUT1(LU,NX-1,NZ ,1.,'A11',AHL,WORK) CALL FDOUT1(LU,NX-1,NZ ,1.,'B11',BHL,WORK) CALL FDOUT1(LU,NX-1,NZ ,1.,'C11',CHL,WORK) CALL FDOUT1(LU,NX ,NZ-1,1.,'A33',AVL,WORK) CALL FDOUT1(LU,NX ,NZ-1,1.,'B33',BVL,WORK) CALL FDOUT1(LU,NX ,NZ-1,1.,'C33',CVL,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'A13',AHS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'B13',BHS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'C13',CHS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'A31',AVS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'B31',BVS,WORK) CALL FDOUT1(LU,NX-1,NZ-1,1.,'C31',CVS,WORK) CALL FDOUT1(LU,NX ,NZ ,1.,'DEN',DEN,WORK) CALL FDOUT1(LU,NX ,NZ ,1.,'QFC',QFC,WORK) C RETURN END C C======================================================================= C C C SUBROUTINE XLEG(KEY,XPT,ZPT,DEL,ALF,BET,RHO) C C Purpose: C Effective (mean) squared velocities for a grid leg C C Input: C KEY=0.....vertical leg C =1.....horizontal leg C XPT,ZPT...x,z coordinates of the grid point C DEL.......length of the leg C C Output: C ALF,BET...P and S squared velocities for the leg C C DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50) DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50) COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C RHO=0. C C An artificial shortening of the leg to avoid numerical problems C !!!!!!Some models possible with 0.01 instead of this 0.03!!!!!! C DELSA=DEL-0.03*DEL cc DELSA=DEL cc DELSA=DEL-0.0001*DEL C C Location of the grid point with respect to a block C DO 100 J=1,NOBL DIST=0. CALL GRPT(XPT,ZPT,J,IPRES,IVIC,DIST,ISEG,IBLVC,IORT) DIST=DIST/DELSA IF(IPRES.EQ.0) GOTO 100 IBL=J GOTO 200 100 CONTINUE WRITE(*,'(A,2F10.3/A)') * '+FDMOD: Gridpoint not located successfully at x,z=',XPT,ZPT,' ' 200 CONTINUE IF(IPRES.EQ.2) GOTO 1000 C C The grid point inside block IBL C RHO=HU(IBL) CALL TWOPT(XPT,ZPT,IBL,DELSA,KEY,0,ISG,IVC,PART) IF(ISG.NE.0) GOTO 400 C C No intersection of the leg with any boundary segment C 300 BET=VS(IBL) ALF=BET*PSRAT(IBL) RETURN C C The leg is intersecting segment ISG C Locating the end point of the leg with respect to block IVC C 400 IF(KEY.EQ.0) THEN XEND=XPT ZEND=ZPT+DELSA ELSE XEND=XPT+DELSA ZEND=ZPT ENDIF DIS2=0. CALL GRPT(XEND,ZEND,IVC,IPR2,IVI2,DIS2,ISE2,IBL2,IOR2) DIS2=DIS2/DELSA IF(IPR2.EQ.0) GOTO 600 IF(IPR2.EQ.1) GOTO 450 IF(IPR2.EQ.2.AND.KEY.EQ.0.AND.IOR2.EQ.0) GOTO 500 IF(IPR2.EQ.2.AND.KEY.EQ.1.AND.IOR2.EQ.1) GOTO 500 C C The end point located inside block IVC C (or on its boundary segment not coinciding with the leg) C 450 VSIN=VS(IBL) VSOU=VS(IVC) VPIN=VSIN*PSRAT(IBL) VPOU=VSOU*PSRAT(IVC) c BET=PART*VSIN+(1.-PART)*VSOU c ALF=PART*VPIN+(1.-PART)*VPOU BET=1./(PART/VSIN+(1.-PART)/VSOU) ALF=1./(PART/VPIN+(1.-PART)/VPOU) RETURN C C The end point located on a vertical/horizontal boundary C between blocks IVC and IBL2 C 500 VSIN=VS(IBL) VSOU1=VS(IVC) VSOU2=VS(IBL2) VSOU=(VSOU1+VSOU2)/2. VPIN=VSIN*PSRAT(IBL) VPOU=(VSOU1*PSRAT(IVC)+VSOU2*PSRAT(IBL2))/2. c BET=PART*VSIN+(1.-PART)*VSOU c ALF=PART*VPIN+(1.-PART)*VPOU BET=1./(PART/VSIN+(1.-PART)/VSOU) ALF=1./(PART/VPIN+(1.-PART)/VPOU) RETURN C C The end point located outside block IVC,inside block IBLNEW C 600 DO 700 J=1,NOBL CALL GRPT(XEND,ZEND,J,IPR3,IVI3,DIS3,ISE3,IBL3,IOR3) IF(IPR3.EQ.0) GOTO 700 IBLNEW=J GOTO 800 700 CONTINUE WRITE(*,'(A,2F10.3/A)') *'+FDMOD: Gridpoint not located successfully at x,z=',XEND,ZEND,' ' C C The end point in block IBLNEW, the leg in blocks IBL, IVC, IBLNEW C 800 CALL TWOPT(XEND,ZEND,IBLNEW,DELSA,KEY,1,ISG2,IVC2, PART2) VSIN=VS(IBL) VSOU1=VS(IVC) VSOU2=VS(IBLNEW) VPIN=VSIN*PSRAT(IBL) VPOU1=VSOU1*PSRAT(IVC) VPOU2=VSOU2*PSRAT(IBLNEW) c BET=PART*VSIN+(1.-PART-PART2)*VSOU1+PART2*VSOU2 c ALF=PART*VPIN+(1.-PART-PART2)*VPOU1+PART2*VPOU2 BET=1./(PART/VSIN+(1.-PART-PART2)/VSOU1+PART2/VSOU2) ALF=1./(PART/VPIN+(1.-PART-PART2)/VPOU1+PART2/VPOU2) RETURN C C The grid point located on the boundary of block IBL C (on segment ISEG, neighbouring to block IBLVC) C Locating the end point with respect to block IBL C 1000 CONTINUE if(ibl.eq.0) rho=hu(iblvc) if(iblvc.eq.0) rho=hu(ibl) if(ibl.ne.0.and.iblvc.ne.0) RHO=(HU(IBL)+HU(IBLVC))/2. IF(KEY.EQ.0) THEN XEND=XPT ZEND=ZPT+DELSA ELSE XEND=XPT+DELSA ZEND=ZPT ENDIF DIS2=0. CALL GRPT(XEND,ZEND,IBL,IPR2,IVI2,DIS2,ISE2,IBL2,IOR2) DIS2=DIS2/DELSA IF(IPR2.EQ.0) GOTO 1800 IF(IPR2.EQ.1) GOTO 1600 IF(IPR2.EQ.2.AND.KEY.EQ.0.AND.IOR2.EQ.0) GOTO 1700 IF(IPR2.EQ.2.AND.KEY.EQ.1.AND.IOR2.EQ.1) GOTO 1700 C C The end point located inside block IBL C (or on its boundary segment not coinciding with the leg) C 1600 BET=VS(IBL) ALF=BET*PSRAT(IBL) RETURN C C The leg is coinciding with a segment C 1700 VSIN=VS(IBL) C VSOU=VS(IBLVC) VSOU=VS(IBL2) VPIN=VSIN*PSRAT(IBL) C VPOU=VSOU*PSRAT(IBLVC) VPOU=VSOU*PSRAT(IBL2) BET=(VSIN+VSOU)/2. ALF=(VPIN+VPOU)/2. RETURN C C The end point located outside block IBL C Locating the end point with respect to block IBLVC C 1800 DIS3=0. CALL GRPT(XEND,ZEND,IBLVC,IPR3,IVI3,DIS3,ISE3,IBL3,IOR3) DIS3=DIS3/DELSA IF(IPR3.EQ.0) GOTO 2100 IF(IPR3.EQ.1) GOTO 1850 IF(IPR3.EQ.2) GOTO 1900 C C The end point located inside block IBLVC C 1850 BET=VS(IBLVC) ALF=BET*PSRAT(IBLVC) RETURN C C The end point located on segment ISE3 of block IBLVC C 1900 IF(KEY.EQ.0.AND.IOR3.EQ.0) GOTO 2000 IF(KEY.EQ.1.AND.IOR3.EQ.1) GOTO 2000 C 1950 BET=VS(IBLVC) ALF=BET*PSRAT(IBLVC) RETURN C C The end point located on a vertical/horizontal boundary C between blocks IBLVC and IBL3 C 2000 VSIN=VS(IBL) VSOU=(VS(IBLVC)+VS(IBL3))/2. VPIN=VS(IBL)*PSRAT(IBL) VPOU=(VS(IBLVC)*PSRAT(IBLVC)+VS(IBL3)*PSRAT(IBL3))/2. c BET=DIST*VSIN+(1.-DIST)*VSOU c ALF=DIST*VPIN+(1.-DIST)*VPOU BET=1./(DIST/VSIN+(1.-DIST)/VSOU) ALF=1./(DIST/VPIN+(1.-DIST)/VPOU) RETURN C C The end point located outside block IBLVC, inside block IBLNEW C 2100 DO 2150 J=1,NOBL DIS5=0. CALL GRPT(XEND,ZEND,J,IPR5,IVI5,DIS5,ISE5,IBL5,IOR5) DIS5=DIS5/DELSA IF(IPR5.EQ.0) GOTO 2150 IBLNEW=J GOTO 2170 2150 CONTINUE 2170 CONTINUE IF(KEY.EQ.0.AND.IORT.EQ.0) GOTO 2200 IF(KEY.EQ.1.AND.IORT.EQ.1) GOTO 2200 GOTO 2400 2200 IF(IPR5.EQ.1) GOTO 2300 IF(IPR5.EQ.2) GOTO 2350 C C The grid point located on a vertical/horizontal boundary C between blocks IBL and IBLVC C The end point inside block IBLNEW C 2300 VSIN=(VS(IBL)+VS(IBLVC))/2. VSOU=VS(IBLNEW) VPIN=(VS(IBL)*PSRAT(IBL)+VS(IBLVC)*PSRAT(IBLVC))/2. VPOU=VSOU*PSRAT(IBLNEW) c BET=DIST*VSIN+(1.-DIST)*VSOU c ALF=DIST*VPIN+(1.-DIST)*VPOU BET=1./(DIST/VSIN+(1.-DIST)/VSOU) ALF=1./(DIST/VPIN+(1.-DIST)/VPOU) RETURN C C The grid point located on a vertical/horizontal boundary C between blocks IBL and IBLVC C The end point located on a vertical/horizontal boundary C between blocks IBLNEW and IBL5 C 2350 VSIN=(VS(IBL)+VS(IBLVC))/2. VSOU=(VS(IBLNEW)+VS(IBL5))/2. VPIN=(VS(IBL)*PSRAT(IBL)+VS(IBLVC)*PSRAT(IBLVC))/2. VPOU=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2. c BET=DIST*VSIN+(1.-DIST)*VSOU c ALF=DIST*VPIN+(1.-DIST)*VPOU BET=1./(DIST/VSIN+(1.-DIST)/VSOU) ALF=1./(DIST/VPIN+(1.-DIST)/VPOU) RETURN C C The end point located inside block IBLNEW, C the leg in blocks IBLNEW and IVC4 (where IVC4=IBL or IBLVC) C 2400 IF(IPR5.EQ.1) GOTO 2600 IF(IPR5.EQ.2.AND.KEY.EQ.0.AND.IOR5.EQ.0) GOTO 2800 IF(IPR5.EQ.2.AND.KEY.EQ.1.AND.IOR5.EQ.1) GOTO 2800 2600 CALL TWOPT(XEND,ZEND,IBLNEW,DELSA,KEY,1,ISG4,IVC4,PART4) VSIN=VS(IVC4) VSOU=VS(IBLNEW) VPIN=VSIN*PSRAT(IVC4) VPOU=VSOU*PSRAT(IBLNEW) c BET=(1.-PART4)*VSIN+PART4*VSOU c ALF=(1.-PART4)*VPIN+PART4*VPOU c some models with inclined boundaries ending at left or right sides c require the following modification: c if(vsin.eq.0.) then c BET=1./(PART4/VSOU) c ALF=1./(PART4/VPOU) c RETURN c endif BET=1./((1.-PART4)/VSIN+PART4/VSOU) ALF=1./((1.-PART4)/VPIN+PART4/VPOU) RETURN 2800 CALL TWOPT(XEND,ZEND,IBL,DELSA,KEY,1,ISG5,IVC5,PART5) IF(ISG5.NE.0) THEN IVXX=IBL PARTX=PART5 GOTO 2900 ELSE CALL TWOPT(XEND,ZEND,IBLVC,DELSA,KEY,1,ISG6,IVC6,PART6) IF(ISG6.NE.0) THEN IVXX=IBLVC PARTX=PART6 GOTO 2900 ELSE c 2850 BET=(VS(IBLNEW)+VS(IBL5))/2. ALF=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2. RETURN ENDIF ENDIF 2900 VSIN=VS(IVXX) VSOU=(VS(IBLNEW)+VS(IBL5))/2. VPIN=VSIN*PSRAT(IVXX) VPOU=(VS(IBLNEW)*PSRAT(IBLNEW)+VS(IBL5)*PSRAT(IBL5))/2. c BET=(1.-PARTX)*VSIN+PARTX*VSOU c ALF=(1.-PARTX)*VPIN+PARTX*VPOU BET=1./((1.-PARTX)/VSIN+PARTX/VSOU) ALF=1./((1.-PARTX)/VPIN+PARTX/VPOU) RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE GRPT(XPT,ZPT,IBL,IPRES,IVIC,DIST,ISEG,IBLVC,IORT) C C Purpose: C Location of grid point (XPT,ZPT) with respect to block IBL C C Input: C XPT,ZPT......x,z coordinates of the grid point C IBL..........the block to be tested C C Output: C IPRES=0......the grid point located outside the block C =1......the grid point located inside the block C =2......the grid point located on the block boundary C C IVIC=0.......the grid point not in vicinity (.LT.DX) of any edge C .ne.0.......the grid point in vicinity of edge IVIC C C DIST.........distance from the edge (for IVIC.ne.0) C C ISEG.........the grid point on segment ISEG (for IPRES=2) C C IBLVC........the block neighbouring to ISEG (for IPRES=2) C C IORT=0.......a vertical segment (for IPRES=2) C =1.......a horizontal segment (for IPRES=2) C =2.......an oblique segment (for IPRES=2) C C DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50) DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50) COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C I=IBL IF(I.EQ.0) THEN IPRES=0 RETURN ENDIF NEDG=NOED(I) DO 200 J=1,NEDG C C Unit position vector (PX,PZ) of grid point (XPT,ZPT) C with respect to the origin of J-th segment of I-th block C PX=XPT-XED(J,I) PZ=ZPT-ZED(J,I) SIZ=SQRT(PX*PX+PZ*PZ) IF(ABS(SIZ).LT.0.00001) GOTO 300 PX=PX/SIZ PZ=PZ/SIZ C C Unit vector (BX,BZ) of J-th segment of I-th block C IF(J.LT.NEDG) THEN BX=XED(J+1,I)-XED(J,I) BZ=ZED(J+1,I)-ZED(J,I) ELSE BX=XED(1,I)-XED(J,I) BZ=ZED(1,I)-ZED(J,I) ENDIF SIZ=SQRT(BX*BX+BZ*BZ) BX=BX/SIZ BZ=BZ/SIZ C C Y-component (=PROD) of the vector product (PX,PZ)x(BX,BZ) C PROD=PZ*BX-PX*BZ IF(PROD.LT.-0.00001) GOTO 200 IF(ABS(PROD).LT.0.00001) GOTO 300 C C The grid point located outside the block C 250 IPRES=0 GOTO 600 C C The grid point possibly located on J-th segment C 300 IPRES=2 ISEG=J IBLVC=NOVIC(J,I) C C End points of J-th segment, orientation of the segment C X1=XED(J,I) Z1=ZED(J,I) IF(J.LT.NEDG) THEN X2=XED(J+1,I) Z2=ZED(J+1,I) ELSE X2=XED(1,I) Z2=ZED(1,I) ENDIF SEGLEN=SQRT((X2-X1)**2+(Z2-Z1)**2) IORT=2 IF(ABS(X1-X2).LT.0.00001) IORT=0 IF(ABS(Z1-Z2).LT.0.00001) IORT=1 C C Distances of the grid point from the end points C DIST1=SQRT((XPT-X1)**2+(ZPT-Z1)**2) DIST2=SQRT((XPT-X2)**2+(ZPT-Z2)**2) ROZ1=DIST1-SEGLEN ROZ2=DIST2-SEGLEN C IF(DIST1.GT.SEGLEN.OR.DIST2.GT.SEGLEN) GOTO 250 IF(ROZ1.GT.0.0001.OR.ROZ2.GT.0.0001) GOTO 250 C C The grid point certainly located on J-th segment C IVIC=0 IF(DIST1.LT.DX) GOTO 400 IF(DIST2.LT.DX) GOTO 500 GOTO 600 400 IVIC=J DIST=DIST1 GOTO 600 500 IF(J.LT.NEDG) THEN IVIC=J+1 ELSE IVIC=1 ENDIF DIST=DIST2 GOTO 600 200 CONTINUE C C The grid point 'on the left' with respect to all segments, C i.e., the grid point located inside I-th block C IPRES=1 C C Distances from all edges of the block C (The program cannot detect proximity of two edges simultan.) C IVIC=0 DO 550 J=1,NEDG DISTX=SQRT((XPT-XED(J,I))**2+(ZPT-ZED(J,I))**2) IF(DISTX.GT.DX) GOTO 550 IVIC=J DIST=DISTX GOTO 600 550 CONTINUE 600 RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C C C SUBROUTINE TWOPT(X,Z,IB,DELSA,KY,KYOR,ISG,IVC,PART) C C Purpose: C Location of a grid leg with respect to boundary of block IB C C Input: C X,Z........x,z coordinates of the grid point C IB.........the block to which the grid point is belonging C DELSA........length of the leg C KY=0.......vertical leg C =1.......horizontal leg C KYOR=0.....leg oriented down, or right, of the grid point C =1.....up or left C C Output: C ISG=0......no intersection of the leg with the block boundary C .ne.0....intersection with ISG-th segment C IVC........the neighbouring block (for ISG.ne.0) C PART.......a relative distance of the grid point from C the intersection (for ISG.ne.0) C C DIMENSION NOED(0:50),VS(0:50),PSRAT(0:50),HU(0:50) DIMENSION XED(10,50),ZED(10,50),NOVIC(10,50) COMMON /BLO1/ NOBL,NOED,VS,PSRAT,DX,HU COMMON /BLO2/ XED,ZED,NOVIC C IF(IB.EQ.0) THEN ISG=0 RETURN ENDIF C C End point (XX,ZZ) of the leg (the leg from (X,Z) to (XX,ZZ)) C IF(KY.EQ.0) GOTO 100 IF(KY.EQ.1) GOTO 200 100 XX=X IF(KYOR.EQ.0) THEN ZZ=Z+DELSA ELSE ZZ=Z-DELSA ENDIF GOTO 300 200 ZZ=Z IF(KYOR.EQ.0) THEN XX=X+DELSA ELSE XX=X-DELSA ENDIF C Intersection of the leg with the block boundary C (legs coinciding with boundaries are not allowed) C 300 NEDG=NOED(IB) DO 500 J=1,NEDG C C The end points (X1,Z1) and (X2,Z2) of J-th segment C X1=XED(J,IB) Z1=ZED(J,IB) IF(J.LT.NEDG) THEN X2=XED(J+1,IB) Z2=ZED(J+1,IB) ELSE X2=XED(1,IB) Z2=ZED(1,IB) ENDIF IF(KY.EQ.0) GOTO 600 IF(KY.EQ.1) GOTO 700 600 IF(X1.LT.X2) THEN XMIN=X1 XMAX=X2 ELSE XMIN=X2 XMAX=X1 ENDIF IF(X.LT.XMIN.OR.X.GT.XMAX) GOTO 500 ZSEC=Z1+(X-X1)*(Z2-Z1)/(X2-X1) IF(KYOR.EQ.0) THEN ZBOT=Z ZTOP=ZZ ELSE ZBOT=ZZ ZTOP=Z ENDIF IF(ZSEC.LT.ZBOT.OR.ZSEC.GT.ZTOP) GOTO 500 ISG=J IVC=NOVIC(J,IB) PART=ABS(ZSEC-Z)/DELSA GOTO 1000 700 IF(Z1.LT.Z2) THEN ZMIN=Z1 ZMAX=Z2 ELSE ZMIN=Z2 ZMAX=Z1 ENDIF IF(Z.LT.ZMIN.OR.Z.GT.ZMAX) GOTO 500 XSEC=X1+(Z-Z1)*(X2-X1)/(Z2-Z1) IF(KYOR.EQ.0) THEN XBOT=X XTOP=XX ELSE XBOT=XX XTOP=X ENDIF IF(XSEC.LT.XBOT.OR.XSEC.GT.XTOP) GOTO 500 ISG=J IVC=NOVIC(J,IB) PART=ABS(XSEC-X)/DELSA GOTO 1000 500 CONTINUE ISG=0 1000 RETURN END C DEBUG SUBCHK C END DEBUG C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'fdaux.for' C fdaux.for INCLUDE 'gmtra.for' C gmtra.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= C