C
C Program MODSRF for coverage of structural interfaces by polygons C and for computation of 2-D slices through a model. C C Version: 7.61 C Date: 2020, March 21 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C....................................................................... C C Structural interfaces: C Structural interfaces are in MODEL package defined by implicit C functions. This is very useful for computation purposes, but not very C useful for visualization of models. For visualization purposes the C interfaces must be expressed in explicit form, e.g. as sets of C polygons composed of points with known coordinates. C Program MODSRF covers structural interfaces by polygons composed C of points with known cartesian coordinates. C Method: The model is decomposed into a set of cubes given by input C parameters Oi, Ni, Di. Points of intersections of cube edges with C structural interfaces are computed together with points of C intersection of structural edges with faces of the cubes. Points C which belong to the same interface are then connected into C polygons for each cube. C This mode is started if all N1, N2, N3 are greater than 1. C C 2-D sections: C If one of the values N1, N2, N3 equal 1, the 2-D sections are C generated. Above mentioned points of intersection are computed C along rectangles given by Oi, Ni, Di. Points which belong to the C same complex block are then connected into polygons for each C rectangle. C C Problems occur, if the model is not consistent, i.e. if it contains C overlapping simple blocks. It is thus recommended to perform the model C consistency check by means of program C MODCHK C before running MODSRF. C C Problems may also occur, when some gridpoints of the basic grid are C located exactly at the structural interfaces or even at the structural C edges, or even when cube edges coincide with structural interfaces or C even with the structural edges. In such situations the program may C collapse or its results may contain huge numerical errors. C To solve such problems, it is reccomended to slightly shift the basic C grid (if possible), or to decrease input parameter ERRSRF. C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'... String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Data specifying the model by means of the MODEL package: C MODEL='string'... Name of the input file with the data specifying C the model. Cartesian coordinates must be used C in file MODEL. C Description of MODEL. C Example of file MODEL. C No default, 'MODEL' must be specified and cannot be blank. C Parameters defining the basic regular rectangular grid: C N1=positive integer... Number of gridpoints of the basic grid C along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints of the basic grid C along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints of the basic grid C along the X3 axis. C Default: N3=1 C The values of N1, N2, N3 also decide about the mode in C which MODSRF runs: C all of N1,N2,N3 greater 1: Computation of polygons along C structural interfaces. C one of N1,N2,N3 equal 1: Computation of polygons along C 2-D section through the model. C O1=real... X1 coordinate of the origin of the grid. C Default: O1=0. C O2=real... X2 coordinate of the origin of the grid. C Default: O2=0. C O3=real... X3 coordinate of the origin of the grid. C Default: O3=0. C D1=real... Grid spacing along the X1 axis. C Default: D1=1. C D2=real... Grid spacing along the X2 axis. C Default: D2=1. C D3=real... Grid spacing along the X3 axis. C Default: D3=1. C The grid intervals may also be negative. C Names of the output files: C VRTX='string'... Name of the output file with vertices of the C polygons. Description of file VRTX. C Default: VRTX='vrtx.out' C PLGN='string'... Name of the output file describing the polygons. C If blank, the file is not generated. C Description of file PLGN. C Default: PLGN='plgn.out' C PLGNS='string'... Name of the file describing the polygons in C terms of the names of the vertices. C Description of file PLGNS. C Default: PLGNS=' ' (the file is not generated) C Parameters specifying the quantities to be written into the file VRTX. C TEXTP='string'... First part of names of vertices. The second C part of the name of a vertex is formed by number giving C its position in the file VRTX. C Default: TEXTP='V' C COLUMN01 to COLUMN69, POWER01 to POWER69, IVALUE01 to IVALUE69: C IVALUEii=integer... An integer value required for some special C values of COLUMNii. See, e.g., description of COLUMNii C in case that COLUMNii='SRF'. C POWERii=real... Power of the quantity to be written in column ii. C COLUMNii='string'... String which specifies the quantity to be C written to the column ii of the file VRTX. First six C columns usually contain coordinates of the vertices and C the normals. Column zero is reserved for names of the C vertices. Following strings are allowed: C ' ' (a space)... Nothing is to be written to the column C ii and to all the following columns. C 'X1'... First coordinate of the vertex. C 'X2'... Second coordinate of the vertex. C 'X3'... Third coordinate of the vertex. C 'NORM1'... First covariant component of the normal to the C interface (or to the 2-D slice) at the vertex. C 'NORM2'... Second covariant component of the normal. C 'NORM3'... Third covariant component of the normal. C 'ISRF'... Index of the interface. C '+ISB'... Index of simple block at positive side of the C interface. C '-ISB'... Index of simple block at negative side of the C interface. C '+ICB','-ICB'... Index of complex block. C '+VP','-VP'... P-wave velocity. C '+VS','-VS'... S-wave velocity. C '+DEN','-DEN'... Density. C '+QP','-QP'... P wave loss factor. C '+QS','-QS'... S wave loss factor. C '+A11','-A11','+A12','-A12','+A22','-A22','+A13','-A13' C '+A23','-A23','+A33','-A33','+A14','-A14','+A24','-A24' C '+A34','-A34','+A44','-A44','+A15','-A15','+A25','-A25' C '+A35','-A35','+A45','-A45','+A55','-A55','+A16','-A16' C '+A26','-A26','+A36','-A36','+A46','-A46','+A56','-A56' C '+A66','-A66'... Reduced (i.e. divided by the density) C anisotropic elastic parameters C (components of the real part of the symmetric 6*6 C stiffness matrix divided by the density). C '+Q11','-Q11','+Q12','-Q12','+Q22','-Q22','+Q13','-Q13' C '+Q23','-Q23','+Q33','-Q33','+Q14','-Q14','+Q24','-Q24' C '+Q34','-Q34','+Q44','-Q44','+Q15','-Q15','+Q25','-Q25' C '+Q35','-Q35','+Q45','-Q45','+Q55','-Q55','+Q16','-Q16' C '+Q26','-Q26','+Q36','-Q36','+Q46','-Q46','+Q56','-Q56' C '+Q66','-Q66'... Reduced (i.e. divided by the density) C imaginary anisotropic elastic parameters C (components of the imaginary part of the C symmetric 6*6 stiffness matrix divided by the density). C 'SRF'... Value of the function describing a surface C with index IVALUEii (IVALUEii must be specified in C addition to the COLUMNii='SRF'). C All strings may be entered either in uppercase or in C lowercase letters. C Defaults: COLUMN01='X1', COLUMN02='X2', COLUMN03='X3', C COLUMN04='NORM1', COLUMN05='NORM2', COLUMN06='NORM3', C COLUMN07='ISRF', COLUMN08 to COLUMN69=' ', C POWER01 to POWER69=1, IVALUE01 to IVALUE69=1. C Upper error bound in the position of points at interfaces: C ERRSRF=real... The upper error bound. C Default: ERRSRF equals the maximum of (D1+D2+D3)/3000. C and MAX(ABS(COOR))/1000000. where COOR are possible C values of coordinates of gridpoints. C The default value is mostly sufficient. C Parameter to decide, whether the polygons composed of points C situated in free space are to be written to the files PLGN and/or C PLGNS. Important mainly for 2-D slices. C FREESRF=real... the polygons in free space are written to the C output files only if FREESRF=1. C Default: FREESRF=0. (polygons in free space not written) C Optional parameters specifying the form of the real quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers... See the description in file C forms.for. C C C Output file VRTX with the vertices: C (1) / (a slash) C (2) For each vertex data (2.1): C (2.1) 'NAME',R1,R2,... / C 'NAME'... Name of the vertex. See parameter TEXTP above. C R1,R2,... /... None to several values terminated by a slash. See C parameters COLUMN01 to COLUMN69 above. C (3) /... A slash at the end of file. C C C Output file PLGN with the polygons specified in terms of C indices of vertices: C (1) For each polygon data (1.1): C (1.1) I1,I2,...,IN,/ C I1,I2,...,IN... Indices of N vertices of the polygon. C The vertices in file VRTX are indexed by positive integers C according to their order. C /... List of vertices is terminated by a slash. C (2) /... A slash at the end of file. C C C Output file PLGNS with the polygons specified in terms of C names of vertices. This enables vertices from several files C to be combined in further application (the names of the vertices must C differ, see the parameter TEXTP above). C (1) For each polygon data (1.1): C (1.1) 'VRTX1','VRTX2',...,'VRTXN',/ C 'VRTX1','VRTX2',...,'VRTXN'... Strings containing the names of N C vertices of the polygon. The names correspond to the C names in file VRTX. C /... List of vertices is terminated by a slash. C (2) /... A slash at the end of file. C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL MSLPIF,MSLPIL,FCTMS,OUTMS,MSGLEG,MSGFAC,MSGCUB,MSGCGP, *MSXFAC,MSCP,MSEROR, *ERROR,WARN,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM1,LOWER,LENGTH, *MODEL1,BLOCK,BLOCKS,SEPAR,CDE,SRFC2,PARM2,PARM3,RKGS LOGICAL MSLPIF,MSLPIL INTEGER LENGTH C MSLPIF,MSLPIL,FCTMS,OUTMS,MSGLEG,MSGFAC,MSGCUB,MSGCGP, C MSXFAC,MSCP,MSEROR... This file. C ERROR,WARN... File error.for. C RSEP1,RSEP3T,RSEP3R,RSEP3I... File sep.for. C FORM1... File forms.for. C LENGTH,LOWER... File length.for. C MODEL1,BLOCK,BLOCKS,SEPAR... File model.for. C CDE... File means.for. C SRFC2... File srfc.for. C PARM2,PARM3... File parm.for. C RKGS... File rkgs.for. C C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Input and output data files: CHARACTER*80 FSEP,FMOD INTEGER LU,LU1,LU2 PARAMETER (LU=1,LU1=1,LU2=2) C Auxiliary storage locations: INTEGER I1,I2,I3,I4,I5,I6,I7,J1,J2,J3,ITER,NPOINT,LEN1,LEN2,LENG INTEGER IL1,IL2,IL3,IL4,IF1,IF2,IF3,IF4,IF5,IF6 INTEGER IX,IP1,IP2,ISHIFT INTEGER ISBOLD,ICBOLD,ISBNEW,ICBNEW,ISRF,ICBPOS,ICBNEG REAL X1,X2,X3,XX(3),Y1,Y2,Y3,YY(3) EQUIVALENCE (X1,XX(1)),(X2,XX(2)),(X3,XX(3)) EQUIVALENCE (Y1,YY(1)),(Y2,YY(2)),(Y3,YY(3)) REAL F(10),FF(10) INTEGER IDUMMY(1),IY(8),IB REAL ERRSRF,DUMMY(1),XOLD1,XOLD(3),DOLD(3),XNEW1,XNEW(3),DNEW(3), * XTMP1,XTMP(3),DTMP(3),XINT1,XINT(3),DINT(3),FRESRF,RAUX,RAUX1 LOGICAL LFREE REAL PRMT(5),AUX(8,3) INTEGER MJPT,NJPT PARAMETER (MJPT=8) INTEGER JPT(6,MJPT) C JPT(1 to 6,i)... Points along the gridface: address of the point, C ICB,ISB,ISRF,ISB,ICB. JPT(1,i).LT.0 means that the point is C stored in reverse order than it is in IRAM. INTEGER MJCN,NJCN,NNJCN(0:7),MJPOL,NJPOL,MLJPOL,NLJPOL,NIND PARAMETER (MJCN=50,MJPOL=10,MLJPOL=10) INTEGER JCN(4,MJCN),JPOL(0:MLJPOL,MJPOL),IND(MJCN) C JCN(1 to 4,i)... Connection i of the face: address of the first C point, address of the second point, index of the interface, C status of the connection. C NNJCN(1:7)... Number of connections for first gridface, second C gridface, ... , sixth gridface, between edges. C JPOL(0,i)... Number of points in polygon i. C JPOL(1 to JPOL(0,i),i)... Points in the polygon. INTEGER MIPTE PARAMETER (MIPTE=5) INTEGER IPTE(MIPTE),NIPTE C IPTE(i)... address of the point to be used as starting point C when searching for the edge. IPTE(i).LT.0 means that the point C is to be used in reverse order than it is stored in array IRAM. INTEGER IVALUE(69) REAL Z(69),POWER(69),VP(10),VS(10),RHO,QP,QS,A(10,21),Q(10,21), * OUTMIN,OUTMAX LOGICAL LCN CHARACTER*20 FORMAT,FORMA1,FORMA2,FORMA3,VRTX,PLGN,PLGNS,TEXTP, * TEXTC(69)*5,TEXTS(MLJPOL),TEXT C C....................................................................... C C Reading name of SEP file with input data: WRITE(*,'(A)') '+MODSRF: Enter input filename: ' FSEP=' ' READ(*,*) FSEP C C Reading all data from the SEP file into the memory: IF (FSEP.EQ.' ') THEN C MODSRF-01 CALL ERROR('MODSRF-01: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C C Reading all the data from file FSEP to the memory C (SEP parameter file form): CALL RSEP1(LU,FSEP) C C Recalling the data specifying grid dimensions C (arguments: Name of value in input data, Variable, Default): CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) IF (D1.LT.0.) THEN O1=O1+(N1-1)*D1 D1=-D1 ENDIF IF (D2.LT.0.) THEN O2=O2+(N2-1)*D2 D2=-D2 ENDIF IF (D3.LT.0.) THEN O3=O3+(N3-1)*D3 D3=-D3 ENDIF IF ((N1.LE.0).OR.(N2.LE.0).OR.(N3.LE.0).OR. * ((D1.EQ.0.).AND.(N1.NE.1)).OR. * ((D2.EQ.0.).AND.(N2.NE.1)).OR. * ((D3.EQ.0.).AND.(N3.NE.1))) THEN C MODSRF-02 CALL ERROR('MODSRF-02: Wrong specification of the grid.') C This specification of the grid may cause problems. Please, C specify D1,D2,D3 nonzero and N1,N2,N3 greater than 0. Di may C equal zero in case that corresponding Ni equals 1. ENDIF C C Reading the data for the model: CALL RSEP3T('MODEL',FMOD,' ') IF (FMOD.EQ.' ') THEN C MODSRF-03 CALL ERROR('MODSRF-03: MODEL not given') C Input file MODEL with the model must be specified. C There is no default filename. ENDIF OPEN(LU,FILE=FMOD,STATUS='OLD') CALL MODEL1(LU) CLOSE(LU) C C Recalling the output filenames: CALL RSEP3T('VRTX',VRTX,'vrtx.out') CALL RSEP3T('PLGN',PLGN,'plgn.out') CALL RSEP3T('PLGNS',PLGNS,' ') C C Recalling the first part of names of points in output file VRTX: CALL RSEP3T('TEXTP',TEXTP,'V') C C Recalling the data specifying the quantities to be written into C the output file with points at structural interfaces: FORMA1='COLUMN00' FORMA2='POWER00' FORMA3='IVALUE00' I1=1 5 CONTINUE FORMA1(8:8)=CHAR(ICHAR('0')+MOD(I1,10)) FORMA2(7:7)=FORMA1(8:8) FORMA3(8:8)=FORMA1(8:8) FORMA1(7:7)=CHAR(ICHAR('0')+I1/10) FORMA2(6:6)=FORMA1(7:7) FORMA3(7:7)=FORMA1(7:7) IF (I1.EQ.1) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'X1') ELSEIF (I1.EQ.2) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'X2') ELSEIF (I1.EQ.3) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'X3') ELSEIF (I1.EQ.4) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'NORM1') ELSEIF (I1.EQ.5) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'NORM2') ELSEIF (I1.EQ.6) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'NORM3') ELSEIF (I1.EQ.7) THEN CALL RSEP3T(FORMA1,TEXTC(I1),'ISRF') ELSE CALL RSEP3T(FORMA1,TEXTC(I1),' ') ENDIF CALL RSEP3R(FORMA2,POWER(I1),1.) CALL RSEP3I(FORMA3,IVALUE(I1),1) IF (TEXTC(I1).NE.' ') THEN CALL LOWER(TEXTC(I1)) IF ((TEXTC(I1).NE.'x1').AND.(TEXTC(I1).NE.'x2').AND. * (TEXTC(I1).NE.'x3').AND.(TEXTC(I1).NE.'norm1').AND. * (TEXTC(I1).NE.'norm2').AND.(TEXTC(I1).NE.'norm3').AND. * (TEXTC(I1).NE.'isrf').AND. * (TEXTC(I1).NE.'+isb').AND.(TEXTC(I1).NE.'-isb').AND. * (TEXTC(I1).NE.'+icb').AND.(TEXTC(I1).NE.'-icb').AND. * (TEXTC(I1).NE.'+vp') .AND.(TEXTC(I1).NE.'-vp') .AND. * (TEXTC(I1).NE.'+vs') .AND.(TEXTC(I1).NE.'-vs') .AND. * (TEXTC(I1).NE.'+den').AND.(TEXTC(I1).NE.'-den').AND. * (TEXTC(I1).NE.'+qp') .AND.(TEXTC(I1).NE.'-qp') .AND. * (TEXTC(I1).NE.'+qs') .AND.(TEXTC(I1).NE.'-qs') .AND. * (TEXTC(I1).NE.'+a11').AND.(TEXTC(I1).NE.'-a11').AND. * (TEXTC(I1).NE.'+a12').AND.(TEXTC(I1).NE.'-a12').AND. * (TEXTC(I1).NE.'+a22').AND.(TEXTC(I1).NE.'-a22').AND. * (TEXTC(I1).NE.'+a13').AND.(TEXTC(I1).NE.'-a13').AND. * (TEXTC(I1).NE.'+a23').AND.(TEXTC(I1).NE.'-a23').AND. * (TEXTC(I1).NE.'+a33').AND.(TEXTC(I1).NE.'-a33').AND. * (TEXTC(I1).NE.'+a14').AND.(TEXTC(I1).NE.'-a14').AND. * (TEXTC(I1).NE.'+a24').AND.(TEXTC(I1).NE.'-a24').AND. * (TEXTC(I1).NE.'+a34').AND.(TEXTC(I1).NE.'-a34').AND. * (TEXTC(I1).NE.'+a44').AND.(TEXTC(I1).NE.'-a44').AND. * (TEXTC(I1).NE.'+a15').AND.(TEXTC(I1).NE.'-a15').AND. * (TEXTC(I1).NE.'+a25').AND.(TEXTC(I1).NE.'-a25').AND. * (TEXTC(I1).NE.'+a35').AND.(TEXTC(I1).NE.'-a35').AND. * (TEXTC(I1).NE.'+a45').AND.(TEXTC(I1).NE.'-a45').AND. * (TEXTC(I1).NE.'+a55').AND.(TEXTC(I1).NE.'-a55').AND. * (TEXTC(I1).NE.'+a16').AND.(TEXTC(I1).NE.'-a16').AND. * (TEXTC(I1).NE.'+a26').AND.(TEXTC(I1).NE.'-a26').AND. * (TEXTC(I1).NE.'+a36').AND.(TEXTC(I1).NE.'-a36').AND. * (TEXTC(I1).NE.'+a46').AND.(TEXTC(I1).NE.'-a46').AND. * (TEXTC(I1).NE.'+a56').AND.(TEXTC(I1).NE.'-a56').AND. * (TEXTC(I1).NE.'+a66').AND.(TEXTC(I1).NE.'-a66').AND. * (TEXTC(I1).NE.'+q11').AND.(TEXTC(I1).NE.'-q11').AND. * (TEXTC(I1).NE.'+q12').AND.(TEXTC(I1).NE.'-q12').AND. * (TEXTC(I1).NE.'+q22').AND.(TEXTC(I1).NE.'-q22').AND. * (TEXTC(I1).NE.'+q13').AND.(TEXTC(I1).NE.'-q13').AND. * (TEXTC(I1).NE.'+q23').AND.(TEXTC(I1).NE.'-q23').AND. * (TEXTC(I1).NE.'+q33').AND.(TEXTC(I1).NE.'-q33').AND. * (TEXTC(I1).NE.'+q14').AND.(TEXTC(I1).NE.'-q14').AND. * (TEXTC(I1).NE.'+q24').AND.(TEXTC(I1).NE.'-q24').AND. * (TEXTC(I1).NE.'+q34').AND.(TEXTC(I1).NE.'-q34').AND. * (TEXTC(I1).NE.'+q44').AND.(TEXTC(I1).NE.'-q44').AND. * (TEXTC(I1).NE.'+q15').AND.(TEXTC(I1).NE.'-q15').AND. * (TEXTC(I1).NE.'+q25').AND.(TEXTC(I1).NE.'-q25').AND. * (TEXTC(I1).NE.'+q35').AND.(TEXTC(I1).NE.'-q35').AND. * (TEXTC(I1).NE.'+q45').AND.(TEXTC(I1).NE.'-q45').AND. * (TEXTC(I1).NE.'+q55').AND.(TEXTC(I1).NE.'-q55').AND. * (TEXTC(I1).NE.'+q16').AND.(TEXTC(I1).NE.'-q16').AND. * (TEXTC(I1).NE.'+q26').AND.(TEXTC(I1).NE.'-q26').AND. * (TEXTC(I1).NE.'+q36').AND.(TEXTC(I1).NE.'-q36').AND. * (TEXTC(I1).NE.'+q46').AND.(TEXTC(I1).NE.'-q46').AND. * (TEXTC(I1).NE.'+q56').AND.(TEXTC(I1).NE.'-q56').AND. * (TEXTC(I1).NE.'+q66').AND.(TEXTC(I1).NE.'-q66').AND. * (TEXTC(I1).NE.'srf')) THEN C MODSRF-04 CALL ERROR('MODSRF-04: Wrong value of COLUMN.') C See allowed values of COLUMNii in the C description of file SEP. ENDIF I1=I1+1 IF (I1.GT.69) THEN CALL RSEP3T('COLUMN70',TEXT,' ') IF (TEXT.NE.' ') THEN C MODSRF-05 CALL ERROR('MODSRF-05: More than 69 COLUMNs.') C Currently up to 69 values of COLUMNii may be specified. C See allowed values of COLUMNii in the C description of file SEP. ENDIF ELSE GOTO 5 ENDIF ENDIF C End of the loop over future columns of the output file. C C Maximum allowed error in the position of interfaces: RAUX=AMAX1(ABS(O1),ABS(O2),ABS(O3), * ABS(O1+FLOAT(N1-1)*D1), * ABS(O2+FLOAT(N2-1)*D2), * ABS(O3+FLOAT(N3-1)*D3)) RAUX=RAUX/1000000. RAUX1=AMAX1(RAUX,(D1+D2+D3)/3000.) CALL RSEP3R('ERRSRF',ERRSRF,RAUX1) IF (ERRSRF.LT.RAUX) THEN C MODSRF-06 CALL WARN('MODSRF-06: Small ERRSRF.') C Too small value of ERRSRF may cause problems with calculation C of the intersections of gridlegs with interfaces. C If the program works, small ERRSRF does not mind. ENDIF C C Parameter to decide about writing of the polygons situated C in free space: CALL RSEP3R('FREESRF',FRESRF,0.) IF (FRESRF.EQ.1.) THEN LFREE=.TRUE. ELSE LFREE=.FALSE. ENDIF C C Preparing numbers of gridpoints, gridlegs, gridfaces C and gridcubes: N11=N1-1 N21=N2-1 N1N2= N1 * N2 N11N2= (N1-1)* N2 N1N21= N1 *(N2-1) N11N21=(N1-1)*(N2-1) NGPS=N1*N2*N3 OLEG=NGPS*2+1 NLEG1=(N1-1)*N2*N3 NLEG2=N1*(N2-1)*N3 NLEG3=N1*N2*(N3-1) NLEG12=NLEG1+NLEG2 NLEG =NLEG12+NLEG3 OFAC=OLEG+NLEG+1 NFAC1=N1*(N2-1)*(N3-1) NFAC2=(N1-1)*N2*(N3-1) NFAC3=(N1-1)*(N2-1)*N3 NFAC12=NFAC1+NFAC2 NFAC=NFAC12+NFAC3 NCUB=(N1-1)*(N2-1)*(N3-1) OPOI=OFAC+NFAC NPOI=OPOI IF (NPOI.GT.MRAM) CALL MSEROR(1) C DO 10, I1=1,OPOI IRAM(I1)=0 10 CONTINUE C C C Loop along all gridpoints, recording indices of blocks: WRITE(*,'(A)') *'+MODSRF: Computing indices of blocks.' DO 13, I3=1,N3 DO 12, I2=1,N2 DO 11, I1=1,N1 IX=(I3-1)*N1N2+(I2-1)*N1+I1 IX=2*(IX-1) X1=O1+FLOAT(I1-1)*D1 X2=O2+FLOAT(I2-1)*D2 X3=O3+FLOAT(I3-1)*D3 CALL BLOCK(XX,0,0,ISRF,ISBNEW,ICBNEW) IRAM(IX+1)=ISBNEW IRAM(IX+2)=ICBNEW 11 CONTINUE 12 CONTINUE 13 CONTINUE C C C Loop along all gridlegs, C searching for intersections with structural interfaces: WRITE(*,'(A)') *'+MODSRF: Computing intersections along gridlegs.' C Address of intersection points on gridleg 0: IRAM(OLEG)=NPOI C Index of coordinate in the direction of gridlegs: I4=1 C Auxiliary quantities: DOLD(1)=1. DOLD(2)=0. DOLD(3)=0. DNEW(1)=1. DNEW(2)=0. DNEW(3)=0. DTMP(1)=1. DTMP(2)=0. DTMP(3)=0. DO 29, I1=1,NLEG IF (I1.EQ.NLEG1+1) THEN I4=2 DOLD(1)=0. DOLD(2)=1. DNEW(1)=0. DNEW(2)=1. DTMP(1)=0. DTMP(2)=1. ENDIF IF (I1.EQ.NLEG12+1) THEN I4=3 DOLD(2)=0. DOLD(3)=1. DNEW(2)=0. DNEW(3)=1. DTMP(2)=0. DTMP(3)=1. ENDIF CALL MSGLEG(I1,IP1,IP2) ISBOLD=IRAM(2*(IP1-1)+1) ISBNEW=IRAM(2*(IP2-1)+1) IF (ISBOLD.NE.ISBNEW) THEN C The points of the gridleg are in different simple blocks: ICBOLD=IRAM(2*(IP1-1)+2) ICBNEW=IRAM(2*(IP2-1)+2) CALL MSCP(IP1,XOLD(1),XOLD(2),XOLD(3)) CALL MSCP(IP2,XNEW(1),XNEW(2),XNEW(3)) XTMP(1)=XNEW(1) XTMP(2)=XNEW(2) XTMP(3)=XNEW(3) C Loop for intersection points along the gridleg: ITER=0 21 CONTINUE ITER=ITER+1 IF (ITER.GT.100) THEN C MODSRF-07 CALL ERROR ('MODSRF-07: More than 100 intersections.') C More than 100 points of intersection of the C gridline element with the structural interfaces. C Check the input data and then contact the author. ENDIF C C Determining the interface between points XOLD, XNEW: IY(4)=ISBOLD IY(5)=ICBOLD IY(6)=0 XTMP(I4)=XNEW(I4) XOLD1=XOLD(I4) XNEW1=XNEW(I4) XTMP1=XNEW1 CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERRSRF, * XOLD1,XOLD1,XOLD,DOLD,XNEW1,XNEW,DNEW, * XTMP1,XTMP,DTMP,XINT1,XINT,DINT) DO 23, I5=1,3 IF (I5.NE.I4) THEN XTMP(I5)=XOLD(I5) XINT(I5)=XOLD(I5) ENDIF 23 CONTINUE IF ((ICBOLD.NE.ICBNEW).AND.(IY(6).EQ.0)) THEN C The interface was not found between the points with C different ICB. This may happen when the NEW point is C situated directly at the interface. C Repeating the search in reverse direction: DOLD(I4)=-DOLD(I4) DNEW(I4)=-DNEW(I4) DTMP(I4)=-DTMP(I4) XTMP1=XOLD1 XTMP(1)=XOLD(1) XTMP(2)=XOLD(2) XTMP(3)=XOLD(3) IY(4)=ISBNEW IY(5)=ICBNEW CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERRSRF, * (XTMP1-ERRSRF),XNEW1,XNEW,DNEW,XOLD1,XOLD,DOLD, * XTMP1,XTMP,DTMP,XINT1,XINT,DINT) DO 25, I5=1,3 IF (I5.NE.I4) THEN XTMP(I5)=XOLD(I5) XINT(I5)=XOLD(I5) ENDIF 25 CONTINUE IF (IY(6).EQ.0) THEN C The interface was not found even in reverse direction. C This may happen when the whole gridleg coincides with an C interface: CALL BLOCK(XOLD,0,ISBOLD,IDUMMY,IB,IDUMMY) IF (IB.NE.ISBOLD) CALL MSEROR(8) CALL BLOCK(XOLD,0,ISBNEW,IDUMMY,IB,IDUMMY) IF (IB.NE.ISBNEW) CALL MSEROR(8) CALL BLOCK(XNEW,0,ISBNEW,IDUMMY,IB,IDUMMY) IF (IB.NE.ISBNEW) CALL MSEROR(8) CALL BLOCK(XNEW,0,ISBOLD,IDUMMY,IB,IDUMMY) IF (IB.NE.ISBOLD) CALL MSEROR(8) C Whole gridleg coincides with an interface: CALL SEPAR(ISBNEW,ISBOLD,IDUMMY,IY(6)) IF (IDUMMY(1).LT.1) CALL MSEROR(8) IY(7)=ISBOLD IY(8)=ICBOLD ENDIF IF (IY(8).NE.ICBOLD) THEN C There should be only one interface between the points. C C MODSRF-08 CALL ERROR('MODSRF-08: More than one interface.') C This error should not appear. Contact the author. ENDIF C Reverting the points to the original order: DOLD(I4)=-DOLD(I4) DNEW(I4)=-DNEW(I4) DTMP(I4)=-DTMP(I4) XTMP1=XNEW1 XTMP(1)=XNEW(1) XTMP(2)=XNEW(2) XTMP(3)=XNEW(3) ISBOLD=IY(7) ICBOLD=IY(8) ISBNEW=IY(4) ICBNEW=IY(5) IY(4)=ISBOLD IY(5)=ICBOLD IY(6)=-IY(6) IY(7)=ISBNEW IY(8)=ICBNEW ENDIF IF (IY(6).NE.0) THEN C Structural interface, recording the point of intersection: C XTMP and XINT are the points of intersection with the C interface, IY(4) and IY(5) are the indices of the simple C block and the complex block before the interface, C IY(6) is the index of the interface, IY(7) and IY(8) are C the indices of the simple block and the complex block C behind the interface: IF (.NOT.MSLPIL(XTMP,I1)) THEN C MODSRF-09 CALL ERROR('MODSRF-09: Point not on the gridleg.') C This error should not appear. Contact the author. ENDIF IF (NPOI+NQPOI.GT.MRAM) CALL MSEROR(1) RAM(NPOI+1)=XTMP(1) RAM(NPOI+2)=XTMP(2) RAM(NPOI+3)=XTMP(3) IRAM(NPOI+4)=IY(5) IRAM(NPOI+5)=IY(4) IRAM(NPOI+6)=IY(6) IRAM(NPOI+7)=IY(7) IRAM(NPOI+8)=IY(8) NPOI=NPOI+NQPOI IF (IY(7).NE.ISBNEW) THEN C IY(7) is the index of simple block behind the interface, C the search for intersection points must continue. C Shifting point XOLD to the point of intersection: XOLD(I4)=XTMP(I4) XOLD1=XOLD(I4) ISBOLD=IY(7) ICBOLD=IY(8) GOTO 21 ENDIF ENDIF C End of the loop for intersection points along the gridleg. ENDIF IRAM(OLEG+I1)=NPOI 29 CONTINUE NPOIL=NPOI C C C Loop along all gridfaces, C connecting the points of intersection of structural interfaces C with gridlegs, if necessary computing points of intersection C of structural edges with gridfaces. WRITE(*,'(A)') *'+MODSRF: Connecting points along gridfaces. ' IF (NCUB.NE.0) THEN OCON=NPOI+(NPOI-OPOI)/2 ELSE OCON=NPOI + (NPOI-OPOI)*3 + NGPS*NQPOI ENDIF NCON=OCON IF (NCON.GT.MRAM) CALL MSEROR(1) C Initialization for RKGS: PRMT(1)=0. PRMT(2)=(D1+D2+D3) PRMT(3)=ERRSRF PRMT(4)=ERRSRF DTMP(1)=1./3. DTMP(2)=DTMP(1) DTMP(3)=DTMP(1) C Address of connections on gridface 0: IRAM(OFAC)=NCON DO 49, I1=1,NFAC C Indices of gridlegs of the gridface: CALL MSGFAC(I1,IL1,IL2,IL3,IL4) C C Forming array with intersection points for the gridface: NJPT=0 DO 31, I2=IRAM(OLEG+IL1-1),IRAM(OLEG+IL1)-NQPOI,NQPOI C I2 points to the start of records for intersection point. NJPT=NJPT+1 IF (NJPT.GT.MJPT) CALL MSEROR(2) JPT(1,NJPT)=I2+NQPOI JPT(2,NJPT)=IRAM(I2+4) JPT(3,NJPT)=IRAM(I2+5) JPT(4,NJPT)=IRAM(I2+6) JPT(5,NJPT)=IRAM(I2+7) JPT(6,NJPT)=IRAM(I2+8) 31 CONTINUE DO 32, I2=IRAM(OLEG+IL2-1),IRAM(OLEG+IL2)-NQPOI,NQPOI C I2 points to the start of records for intersection point. NJPT=NJPT+1 IF (NJPT.GT.MJPT) CALL MSEROR(2) JPT(1,NJPT)=I2+NQPOI JPT(2,NJPT)=IRAM(I2+4) JPT(3,NJPT)=IRAM(I2+5) JPT(4,NJPT)=IRAM(I2+6) JPT(5,NJPT)=IRAM(I2+7) JPT(6,NJPT)=IRAM(I2+8) 32 CONTINUE DO 33, I2=IRAM(OLEG+IL3)-NQPOI,IRAM(OLEG+IL3-1),-NQPOI C I2 points to the start of records for intersection point. NJPT=NJPT+1 IF (NJPT.GT.MJPT) CALL MSEROR(2) JPT(1,NJPT)=-(I2+NQPOI) JPT(2,NJPT)=IRAM(I2+8) JPT(3,NJPT)=IRAM(I2+7) JPT(4,NJPT)=-IRAM(I2+6) JPT(5,NJPT)=IRAM(I2+5) JPT(6,NJPT)=IRAM(I2+4) 33 CONTINUE DO 34, I2=IRAM(OLEG+IL4)-NQPOI,IRAM(OLEG+IL4-1),-NQPOI C I2 points to the start of records for intersection point. NJPT=NJPT+1 IF (NJPT.GT.MJPT) CALL MSEROR(2) JPT(1,NJPT)=-(I2+NQPOI) JPT(2,NJPT)=IRAM(I2+8) JPT(3,NJPT)=IRAM(I2+7) JPT(4,NJPT)=-IRAM(I2+6) JPT(5,NJPT)=IRAM(I2+5) JPT(6,NJPT)=IRAM(I2+4) 34 CONTINUE C C Making connections between points of the gridface: DO 46, I2=1,NJPT DO 36, I3=I2+1,NJPT IF ((JPT(2,I2).EQ.JPT(6,I3)).AND. * (JPT(6,I2).EQ.JPT(2,I3))) THEN C Same indices of complex blocks: IF ((JPT(3,I2).EQ.JPT(5,I3)).AND. * (JPT(5,I2).EQ.JPT(3,I3))) THEN C Same indices of simple blocks: IF (IABS(JPT(4,I2)).EQ.IABS(JPT(4,I3))) THEN C Same interface - connection found: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (JPT(4,I2).LT.0) THEN IRAM(NCON+1)=IABS(JPT(1,I2)) IRAM(NCON+2)=IABS(JPT(1,I3)) IRAM(NCON+3)=JPT(4,I2) NCON=NCON+3 ELSE IRAM(NCON+1)=IABS(JPT(1,I3)) IRAM(NCON+2)=IABS(JPT(1,I2)) IRAM(NCON+3)=JPT(4,I3) NCON=NCON+3 ENDIF GOTO 45 ENDIF ENDIF ENDIF 36 CONTINUE DO 37, I3=IRAM(OFAC+I1-1),NCON-3,3 IF ((IRAM(I3+1).EQ.IABS(JPT(1,I2))).OR. * (IRAM(I3+2).EQ.IABS(JPT(1,I2)))) C The point I2 is already connected: * GOTO 45 37 CONTINUE C Point I2 is not connected, looking for point of intersection C of structural edges with gridface. (Following the intersection C of the interface with the gridface to the edge.) IFACE=I1 IPTE(1)=JPT(1,I2) NIPTE=1 C Loop for identification of all edges connected with point I2: 39 CONTINUE C Initiating the search for the edge: I3=IABS(IPTE(NIPTE))-NQPOI X1= RAM(I3+1) X2= RAM(I3+2) X3= RAM(I3+3) IF (IPTE(NIPTE).LT.0) THEN ICBO1=IRAM(I3+8) ISBO1=IRAM(I3+7) ISRFO=-IRAM(I3+6) ISBO2=IRAM(I3+5) ICBO2=IRAM(I3+4) ELSE ICBO1=IRAM(I3+4) ISBO1=IRAM(I3+5) ISRFO=IRAM(I3+6) ISBO2=IRAM(I3+7) ICBO2=IRAM(I3+8) ENDIF NIPTE=NIPTE-1 C Computing the point of intersection of the edge C with the gridface: CALL RKGS(PRMT,XX,DTMP,3,I7,FCTMS,OUTMS,AUX) IF (PRMT(5).EQ.2.) THEN C An intersection point was not found within the gridface. C This may happen e.g. when the interface is crossed by C the interface, which separates two (four) simple blocks C of the same complex block(s). I3=I3+NQPOI C Looking for the connection along the given interface: DO 394, I4=1,NJPT IF ((ICBO2.EQ.JPT(2,I4)).AND. * (ICBO1.EQ.JPT(6,I4)).AND. * (IABS(ISRFO).EQ.IABS(JPT(4,I4)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (JPT(4,I4).LT.0) THEN IRAM(NCON+1)=IABS(JPT(1,I4)) IRAM(NCON+2)=I3 IRAM(NCON+3)=JPT(4,I4) ELSE IRAM(NCON+1)=I3 IRAM(NCON+2)=IABS(JPT(1,I4)) IRAM(NCON+3)=-JPT(4,I4) ENDIF NCON=NCON+3 GOTO 398 ENDIF 394 CONTINUE C MODSRF-10 CALL ERROR('MODSRF-10: Point cannot be connected.') C An intersection point was not found within the gridface, C and the point cannot be connected with the points on C the gridlegs. C This error should not appear. Contact the author. 398 CONTINUE ELSEIF (PRMT(5).EQ.1.) THEN C The intersection point was found by RKGS. C Storing the first point on the edge: IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1) RAM(NPOI+1)=X1 RAM(NPOI+2)=X2 RAM(NPOI+3)=X3 IRAM(NPOI+4)=ICBO1 IRAM(NPOI+5)=ISBO1 IRAM(NPOI+6)=ISRFO IRAM(NPOI+7)=ISBO2 IRAM(NPOI+8)=ICBO2 NPOI=NPOI+NQPOI J1=NPOI C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (ISRFO.LT.0) THEN IRAM(NCON+1)=IABS(IPTE(NIPTE+1)) IRAM(NCON+2)=NPOI IRAM(NCON+3)=ISRFO ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=IABS(IPTE(NIPTE+1)) IRAM(NCON+3)=-ISRFO ENDIF NCON=NCON+3 C IF (ICBO1.NE.ICBN1) THEN C Storing the second point on the edge: IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1) RAM(NPOI+1)=X1 RAM(NPOI+2)=X2 RAM(NPOI+3)=X3 IRAM(NPOI+4)=ICBN1 IRAM(NPOI+5)=ISBN1 IRAM(NPOI+6)=-ISRFN1 IRAM(NPOI+7)=ISBO1 IRAM(NPOI+8)=ICBO1 NPOI=NPOI+NQPOI IF (NCUB.EQ.0) THEN C 'Connection' between first and second point on C the edge: IF (NCON+3.GT.MRAM) CALL MSEROR(1) J2=NPOI-NQPOI IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2 J3=NPOI IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3 IRAM(NCON+1)=J2 IRAM(NCON+2)=J3 IRAM(NCON+3)=0 NCON=NCON+3 ENDIF DO 40, I3=1,NJPT IF ((ICBN1.EQ.JPT(2,I3)).AND. * (ISBN1.EQ.JPT(3,I3)).AND. * (ISBO1.EQ.JPT(5,I3)).AND. * (ICBO1.EQ.JPT(6,I3)).AND. * (IABS(ISRFN1).EQ.IABS(JPT(4,I3)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (JPT(4,I3).LT.0) THEN IRAM(NCON+1)=IABS(JPT(1,I3)) IRAM(NCON+2)=NPOI IRAM(NCON+3)=JPT(4,I3) ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=IABS(JPT(1,I3)) IRAM(NCON+3)=-JPT(4,I3) ENDIF NCON=NCON+3 GOTO 405 ENDIF 40 CONTINUE DO 401, I3=1,NIPTE I4=IABS(IPTE(I3))-NQPOI IF ((ICBN1.EQ.IRAM(I4+8)).AND. * (ISBN1.EQ.IRAM(I4+7)).AND. * (ISBO1.EQ.IRAM(I4+5)).AND. * (ICBO1.EQ.IRAM(I4+4)).AND. * (IABS(ISRFN1).EQ.IABS(IRAM(I4+6)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (IRAM(I4+6).LT.0) THEN IRAM(NCON+1)=I4+NQPOI IRAM(NCON+2)=NPOI IRAM(NCON+3)=IRAM(I4+6) ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=I4+NQPOI IRAM(NCON+3)=-IRAM(I4+6) ENDIF NCON=NCON+3 GOTO 405 ENDIF 401 CONTINUE C More than one edge, current point will become a starting C point of the search for new edge: NIPTE=NIPTE+1 IF (NIPTE.GT.MIPTE) CALL MSEROR(3) IPTE(NIPTE)=-NPOI 405 CONTINUE ENDIF C IF (ICBN2.NE.ICBN1) THEN C Storing the third point on the edge: IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1) RAM(NPOI+1)=X1 RAM(NPOI+2)=X2 RAM(NPOI+3)=X3 IRAM(NPOI+4)=ICBN2 IRAM(NPOI+5)=ISBN2 IRAM(NPOI+6)=-ISRFO IRAM(NPOI+7)=ISBN1 IRAM(NPOI+8)=ICBN1 NPOI=NPOI+NQPOI IF (NCUB.EQ.0) THEN C 'Connection' between second and third point on C the edge: IF (NCON+3.GT.MRAM) CALL MSEROR(1) J2=NPOI-NQPOI IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2 J3=NPOI IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3 IRAM(NCON+1)=J2 IRAM(NCON+2)=J3 IRAM(NCON+3)=0 NCON=NCON+3 ENDIF DO 42, I3=1,NJPT IF ((ICBN2.EQ.JPT(2,I3)).AND. * (ISBN2.EQ.JPT(3,I3)).AND. * (ISBN1.EQ.JPT(5,I3)).AND. * (ICBN1.EQ.JPT(6,I3)).AND. * (IABS(ISRFO).EQ.IABS(JPT(4,I3)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (JPT(4,I3).LT.0) THEN IRAM(NCON+1)=IABS(JPT(1,I3)) IRAM(NCON+2)=NPOI IRAM(NCON+3)=JPT(4,I3) ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=IABS(JPT(1,I3)) IRAM(NCON+3)=-JPT(4,I3) ENDIF NCON=NCON+3 GOTO 425 ENDIF 42 CONTINUE DO 421, I3=1,NIPTE I4=IABS(IPTE(I3))-NQPOI IF ((ICBN2.EQ.IRAM(I4+8)).AND. * (ISBN2.EQ.IRAM(I4+7)).AND. * (ISBN1.EQ.IRAM(I4+5)).AND. * (ICBN1.EQ.IRAM(I4+4)).AND. * (IABS(ISRFO).EQ.IABS(IRAM(I4+6)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (IRAM(I4+6).LT.0) THEN IRAM(NCON+1)=I4+NQPOI IRAM(NCON+2)=NPOI IRAM(NCON+3)=IRAM(I4+6) ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=I4+NQPOI IRAM(NCON+3)=-IRAM(I4+6) ENDIF NCON=NCON+3 GOTO 425 ENDIF 421 CONTINUE C More than one edge, current point will become a starting C point of the search for new edge: NIPTE=NIPTE+1 IF (NIPTE.GT.MIPTE) CALL MSEROR(3) IPTE(NIPTE)=-NPOI 425 CONTINUE ENDIF C IF (ICBO2.NE.ICBN2) THEN C Storing the fourth point on the edge: IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1) RAM(NPOI+1)=X1 RAM(NPOI+2)=X2 RAM(NPOI+3)=X3 IRAM(NPOI+4)=ICBO2 IRAM(NPOI+5)=ISBO2 IRAM(NPOI+6)=ISRFN2 IRAM(NPOI+7)=ISBN2 IRAM(NPOI+8)=ICBN2 NPOI=NPOI+NQPOI IF (NCUB.EQ.0) THEN C 'Connection' between third and fourth point on C the edge: IF (NCON+3.GT.MRAM) CALL MSEROR(1) J2=NPOI-NQPOI IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2 J3=NPOI IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3 IRAM(NCON+1)=J2 IRAM(NCON+2)=J3 IRAM(NCON+3)=0 NCON=NCON+3 ENDIF DO 41, I3=1,NJPT IF ((ICBO2.EQ.JPT(2,I3)).AND. * (ISBO2.EQ.JPT(3,I3)).AND. * (ISBN2.EQ.JPT(5,I3)).AND. * (ICBN2.EQ.JPT(6,I3)).AND. * (IABS(ISRFN2).EQ.IABS(JPT(4,I3)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (JPT(4,I3).LT.0) THEN IRAM(NCON+1)=IABS(JPT(1,I3)) IRAM(NCON+2)=NPOI IRAM(NCON+3)=JPT(4,I3) ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=IABS(JPT(1,I3)) IRAM(NCON+3)=-JPT(4,I3) ENDIF NCON=NCON+3 GOTO 415 ENDIF 41 CONTINUE DO 411, I3=1,NIPTE I4=IABS(IPTE(I3))-NQPOI IF ((ICBO2.EQ.IRAM(I4+8)).AND. * (ISBO2.EQ.IRAM(I4+7)).AND. * (ISBN2.EQ.IRAM(I4+5)).AND. * (ICBN2.EQ.IRAM(I4+4)).AND. * (IABS(ISRFN2).EQ.IABS(IRAM(I4+6)))) THEN C Recording the connection: IF (NCON+3.GT.MRAM) CALL MSEROR(1) IF (IRAM(I4+6).LT.0) THEN IRAM(NCON+1)=I4+NQPOI IRAM(NCON+2)=NPOI IRAM(NCON+3)=IRAM(I4+6) ELSE IRAM(NCON+1)=NPOI IRAM(NCON+2)=I4+NQPOI IRAM(NCON+3)=-IRAM(I4+6) ENDIF NCON=NCON+3 GOTO 415 ENDIF 411 CONTINUE C More than one edge, current point will become a starting C point of the search for new edge: NIPTE=NIPTE+1 IF (NIPTE.GT.MIPTE) CALL MSEROR(3) IPTE(NIPTE)=-NPOI 415 CONTINUE ENDIF C IF (NCUB.EQ.0) THEN C 'Connection' between fourth and first point on the edge: IF (NCON+3.GT.MRAM) CALL MSEROR(1) J2=NPOI IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2 J3=J1 IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3 IRAM(NCON+1)=J2 IRAM(NCON+2)=J3 IRAM(NCON+3)=0 NCON=NCON+3 ENDIF ELSE C MODSRF-11 CALL ERROR('MODSRF-11: Wrong value of PRMT(5).') C PRMT(5) should equal either 1. or 2. after RKGS is called. C This error should not appear. Contact the author. ENDIF IF (NIPTE.GE.1) C Continuing with the search for next edge: * GOTO 39 C End of the loop for identification of all edges connected with C point I2. C C No more edges, the point is connected. 45 CONTINUE C Continuing with next point of the gridface: 46 CONTINUE IRAM(OFAC+I1)=NCON 49 CONTINUE NPOIE=NPOI C C IF (NCUB.NE.0) THEN C Loop along all gridcubes, C merging the connections between the points on structural C interfaces into polygons approximating the interfaces. WRITE(*,'(A)') *'+MODSRF: Making polygons approximating the interfaces.' NPOL=MRAM+1 DO 100, I1=1,NCUB C Indices of gridfaces of the gridcube: CALL MSGCUB(I1,IF1,IF2,IF3,IF4,IF5,IF6) C C Forming array with connections for the gridcube: NJCN=0 NNJCN(0)=NJCN DO 51, I2=IRAM(OFAC+IF1-1),IRAM(OFAC+IF1)-3,3 C Connections of the gridface IF1: NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+1) JCN(2,NJCN)=IRAM(I2+2) JCN(3,NJCN)=IRAM(I2+3) JCN(4,NJCN)=0 51 CONTINUE NNJCN(1)=NJCN DO 52, I2=IRAM(OFAC+IF2-1),IRAM(OFAC+IF2)-3,3 C Connections of the gridface IF2: NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+1) JCN(2,NJCN)=IRAM(I2+2) JCN(3,NJCN)=IRAM(I2+3) JCN(4,NJCN)=0 52 CONTINUE NNJCN(2)=NJCN DO 53, I2=IRAM(OFAC+IF3-1),IRAM(OFAC+IF3)-3,3 C Connections of the gridface IF3: NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+1) JCN(2,NJCN)=IRAM(I2+2) JCN(3,NJCN)=IRAM(I2+3) JCN(4,NJCN)=0 53 CONTINUE NNJCN(3)=NJCN DO 54, I2=IRAM(OFAC+IF4-1),IRAM(OFAC+IF4)-3,3 C Connections of the gridface IF4: NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+2) JCN(2,NJCN)=IRAM(I2+1) JCN(3,NJCN)=IRAM(I2+3) JCN(4,NJCN)=0 54 CONTINUE NNJCN(4)=NJCN DO 55, I2=IRAM(OFAC+IF5-1),IRAM(OFAC+IF5)-3,3 C Connections of the gridface IF5: NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+2) JCN(2,NJCN)=IRAM(I2+1) JCN(3,NJCN)=IRAM(I2+3) JCN(4,NJCN)=0 55 CONTINUE NNJCN(5)=NJCN DO 56, I2=IRAM(OFAC+IF6-1),IRAM(OFAC+IF6)-3,3 C Connections of the gridface IF6: NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+2) JCN(2,NJCN)=IRAM(I2+1) JCN(3,NJCN)=IRAM(I2+3) JCN(4,NJCN)=0 56 CONTINUE NNJCN(6)=NJCN DO 58, I2=1,NNJCN(6) C Connections between points on edges: IF (JCN(2,I2).GT.NPOIL) THEN ICBO1=IRAM(JCN(2,I2)-NQPOI+4) ICBO2=IRAM(JCN(2,I2)-NQPOI+8) LCN=.FALSE. DO 57, I3=1,NNJCN(6) IF (I3.NE.I2) THEN IF ((JCN(1,I3).GT.NPOIL).AND. * (IABS(JCN(3,I2)).EQ.IABS(JCN(3,I3)))) THEN ICBN1=IRAM(JCN(1,I3)-NQPOI+4) ICBN2=IRAM(JCN(1,I3)-NQPOI+8) IF (((ICBO1.EQ.ICBN1).AND.(ICBO2.EQ.ICBN2)).OR. * ((ICBO2.EQ.ICBN1).AND.(ICBO1.EQ.ICBN2))) THEN NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=JCN(2,I2) JCN(2,NJCN)=JCN(1,I3) JCN(3,NJCN)=JCN(3,I3) JCN(4,NJCN)=0 LCN=.TRUE. ENDIF ENDIF ENDIF 57 CONTINUE IF (.NOT.LCN) THEN C Connection for edge not found, this connection is to be C canceled by means of creation of linear 'triangle': NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=JCN(2,I2) JCN(2,NJCN)=JCN(1,I2) JCN(3,NJCN)=JCN(3,I2) JCN(4,NJCN)=0 ENDIF ENDIF 58 CONTINUE NNJCN(7)=NJCN C C Looking for polygons on the faces of the cube. NJPOL=0 C Loop for faces of the cube: DO 70, I2=1,6 C Initiating indices of connections: DO 61, I3=1,MJCN IND(I3)=0 61 CONTINUE NIND=NNJCN(I2)-NNJCN(I2-1) C Loop along connections: I3=1 62 CONTINUE IF (NJPOL+1.GT.MJPOL) CALL MSEROR(6) C Search for first unused connection: IF (I3.LE.NIND-2) THEN IF (IND(I3).GE.0) THEN C Unused connection, recording first two points C of the polygon: NLJPOL=2 JPOL(1,NJPOL+1)=JCN(1,NNJCN(I2-1)+I3) JPOL(2,NJPOL+1)=JCN(2,NNJCN(I2-1)+I3) IND(I3)=-1 ELSE I3=I3+1 GOTO 62 ENDIF C Search for next point of polygon: 63 CONTINUE DO 69, I4=I3+1,NIND IF ( JCN(1,NNJCN(I2-1)+I4) .EQ. * JPOL(NLJPOL,NJPOL+1) ) THEN C Recording next point of polygon: NLJPOL=NLJPOL+1 IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7) JPOL(NLJPOL,NJPOL+1)=JCN(2,NNJCN(I2-1)+I4) IND(I4)=-1 DO 68, I5=NLJPOL-3,1,-1 C Examining whether the polygon is closed: IF (JPOL(I5,NJPOL+1).EQ.JPOL(NLJPOL,NJPOL+1)) THEN C Polygon is closed. C Recording the polygon: NJPOL=NJPOL+1 NLJPOL=NLJPOL-I5 JPOL(0,NJPOL)=NLJPOL DO 65, I6=1,NLJPOL JPOL(I6,NJPOL)=JPOL(I5+I6,NJPOL) 65 CONTINUE C Removing used starting points: DO 67, I6=1,NLJPOL IF (I6.EQ.1) THEN J1=JPOL(NLJPOL,NJPOL) ELSE J1=JPOL(I6-1,NJPOL) ENDIF J2=JPOL(I6,NJPOL) IF (I6.EQ.NJPOL) THEN J3=JPOL(1,NJPOL) ELSE J3=JPOL(I6+1,NJPOL) ENDIF 67 CONTINUE C Continuing with next polygon: GOTO 62 ENDIF 68 CONTINUE GOTO 63 ENDIF 69 CONTINUE C Polygon is opened, searching for other polygon: GOTO 62 ENDIF C No more unused connections, C no more polygons on this face of the cube. C End of the loop along connections. 70 CONTINUE C C Two or more polygons on the faces of the cube: IF (NJPOL.GE.2) THEN C MODSRF-12 CALL ERROR('MODSRF-12: Two or more polygons on gridface.') C This situation cannot be handled by current version of C MODSRF. Try to reduce the gridstep (D1,D2,D3). ENDIF C C If there is one polygon on the faces of the cube, C the polygon is to be recorded as a polygon approximating C a structural interface: IF (NJPOL.EQ.1) THEN C Searching, whether the polygon is to be recorded: C Loop for already recorded polygons: I2=MRAM+1 73 CONTINUE IF (I2.GT.NPOL) THEN J1=IRAM(I2-1) I2=I2-J1-1 IF (J1.EQ.JPOL(0,1)) THEN DO 74, I3=1,J1 IF (IRAM(I2+I3).NE.JPOL(I3,1)) GOTO 73 74 CONTINUE C The polygon is already recorded: GOTO 77 ELSE GOTO 73 ENDIF ENDIF C End of the loop for already recorded polygons. C Recording the polygon: IF (NPOL-JPOL(0,1)-1.LE.NCON) CALL MSEROR(1) NPOL=NPOL-1 IRAM(NPOL)=JPOL(0,1) DO 75, I2=1,JPOL(0,1) NPOL=NPOL-1 IRAM(NPOL)=JPOL(I2,1) 75 CONTINUE 77 CONTINUE ENDIF C C Recording polygons approximating structural interfaces: C Marking wrong and sure connections: CALL MSJCN(JCN,NJCN,MJCN) C C Recording the polygons: NJPOL=0 80 CONTINUE C Loop for adding new polygons C Initializing a polygon: DO 81, I2=1,NJCN IF (JCN(4,I2).EQ.1) THEN NJPOL=NJPOL+1 IF (NJPOL.GT.MJPOL) CALL MSEROR(6) NLJPOL=2 JPOL(0,NJPOL)=0 JPOL(1,NJPOL)=JCN(1,I2) JPOL(2,NJPOL)=JCN(2,I2) JCN(4,I2)=-2 GOTO 83 ENDIF 81 CONTINUE C No other point to initialize a polygon: GOTO 97 C 83 CONTINUE C Loop for adding points to the polygon: C C Removing connections, which would cause overlooping C of the polygon in the current position: DO 87, I2=1,NJCN IF (JCN(4,I2).EQ.0) THEN IF (JCN(1,I2).EQ.JPOL(NLJPOL,NJPOL)) THEN DO 84, I3=2,NLJPOL IF (JCN(2,I2).EQ.JPOL(I3,NJPOL)) THEN C This connection would cause overlooping C of the end of the polygon: JCN(4,I2)=-1 CALL MSJCN(JCN,NJCN,MJCN) GOTO 86 ENDIF 84 CONTINUE ELSEIF (JCN(2,I2).EQ.JPOL(1,NJPOL)) THEN DO 85, I3=2,NLJPOL IF (JCN(1,I2).EQ.JPOL(I3,NJPOL)) THEN C This connection would cause overlooping C of the beginning of the polygon: JCN(4,I2)=-1 CALL MSJCN(JCN,NJCN,MJCN) GOTO 86 ENDIF 85 CONTINUE ENDIF ENDIF 86 CONTINUE 87 CONTINUE C C Looking for the point, which might be added C to the polygon in the current position: DO 92, I2=1,NJCN C Trying to add a point to the end of the polygon: IF ((JCN(1,I2).EQ.JPOL(NLJPOL,NJPOL)).AND. * (JCN(4,I2).EQ.1)) THEN C Point I2 may be added to the polygon. IF (JCN(2,I2).EQ.JPOL(1,NJPOL)) THEN C This point closes the polygon: JPOL(0,NJPOL)=NLJPOL JCN(4,I2)=-2 GOTO 80 ELSE C Adding a point to the polygon: NLJPOL=NLJPOL+1 IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7) JPOL(NLJPOL,NJPOL)=JCN(2,I2) JCN(4,I2)=-2 GOTO 83 ENDIF ENDIF 92 CONTINUE DO 94, I2=1,NJCN C Trying to add a point to the beginning of the polygon: IF ((JCN(2,I2).EQ.JPOL(1,NJPOL)).AND. * (JCN(4,I2).EQ.1)) THEN C Point I2 may be added to the polygon. IF (JCN(1,I2).EQ.JPOL(NLJPOL,NJPOL)) THEN C This point closes the polygon: JPOL(0,NJPOL)=NLJPOL JCN(4,I2)=-2 GOTO 80 ELSE C Adding a point to the polygon: NLJPOL=NLJPOL+1 IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7) DO 93, I3=NLJPOL,2,-1 JPOL(I3,NJPOL)=JPOL(I3-1,NJPOL) 93 CONTINUE JPOL(1,NJPOL)=JCN(1,I2) JCN(4,I2)=-2 GOTO 83 ENDIF ENDIF 94 CONTINUE DO 95, I2=1,NJCN C Trying to find and cancel the connection between C the last and the first point of the polygon: IF ((JCN(1,I2).EQ.JPOL(NLJPOL,NJPOL)).AND. * (JCN(2,I2).EQ.JPOL(1,NJPOL)).AND. * (JCN(4,I2).EQ.0)) THEN JCN(4,I2)=-1 CALL MSJCN(JCN,NJCN,MJCN) GOTO 83 ENDIF 95 CONTINUE C MODSRF-13 CALL ERROR ('MODSRF-13: Polygon opened.') C This error should not appear. Contact the author. C End of the loop for adding new points to the polygon. C End of the loop for adding new polygons. C 97 CONTINUE C Writing polygons to IRAM: DO 99, I2=1,NJPOL IF (JPOL(0,I2).GT.2) THEN IF (NPOL-JPOL(0,I2)-1.LE.NCON) CALL MSEROR(1) NPOL=NPOL-1 IRAM(NPOL)=JPOL(0,I2) DO 98, I3=1,JPOL(0,I2) NPOL=NPOL-1 IRAM(NPOL)=JPOL(I3,I2) 98 CONTINUE ENDIF 99 CONTINUE NJPOL=0 C C End of the loop over gridcubes. 100 CONTINUE C C ELSE C Loop over rectangles of the 2-D grid, making polygons along the C 2-D slice through the model. WRITE(*,'(A)') *'+MODSRF: Making polygons along the 2-D slice. ' NPOL=MRAM+1 C C Rewriting points on interfaces as points on + side of interface: ISHIFT=NPOIE-OPOI IF (NPOI+ISHIFT.GT.OCON) CALL MSEROR(1) DO 102, I1=OPOI,NPOIE-NQPOI,NQPOI RAM(NPOI+1)=RAM(I1+1) RAM(NPOI+2)=RAM(I1+2) RAM(NPOI+3)=RAM(I1+3) IF (IRAM(I1+6).LT.0) THEN IRAM(NPOI+4)=IRAM(I1+8) IRAM(NPOI+5)=IRAM(I1+7) IRAM(NPOI+6)=0 IRAM(NPOI+7)=IRAM(I1+7) IRAM(NPOI+8)=IRAM(I1+8) ELSE IRAM(NPOI+4)=IRAM(I1+4) IRAM(NPOI+5)=IRAM(I1+5) IRAM(NPOI+6)=0 IRAM(NPOI+7)=IRAM(I1+5) IRAM(NPOI+8)=IRAM(I1+4) ENDIF NPOI=NPOI+NQPOI 102 CONTINUE NPOILP=NPOIL+ISHIFT NPOIEP=NPOIE+ISHIFT C C Rewriting points on interfaces as points on - side of interface: DO 104, I1=OPOI,NPOIE-NQPOI,NQPOI IF (IRAM(I1+6).LT.0) THEN IRAM(I1+7)=IRAM(I1+5) IRAM(I1+8)=IRAM(I1+4) ELSE IRAM(I1+4)=IRAM(I1+8) IRAM(I1+5)=IRAM(I1+7) ENDIF IRAM(NPOI+6)=0 104 CONTINUE C C Recording gridpoints to the array of points: IF (NPOI+NGPS*NQPOI.GT.OCON) CALL MSEROR(1) DO 106, I1=1,NGPS CALL MSCP(I1,X1,X2,X3) RAM(NPOI+1)=X1 RAM(NPOI+2)=X2 RAM(NPOI+3)=X3 IRAM(NPOI+4)=IRAM(2*(I1-1)+2) IRAM(NPOI+5)=IRAM(2*(I1-1)+1) IRAM(NPOI+6)=0 IRAM(NPOI+7)=IRAM(2*(I1-1)+1) IRAM(NPOI+8)=IRAM(2*(I1-1)+2) NPOI=NPOI+NQPOI 106 CONTINUE C Loop over gridfaces: DO 200, I1=1,NFAC C Forming array with connections for the gridface: NJCN=0 C Connections inside the gridface: DO 110, I2=IRAM(OFAC+I1-1),IRAM(OFAC+I1)-3,3 C Connections of the gridface I1: IF (IRAM(I2+3).NE.0) THEN NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+1) JCN(2,NJCN)=IRAM(I2+2) JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=IRAM(I2+2)+ISHIFT JCN(2,NJCN)=IRAM(I2+1)+ISHIFT JCN(3,NJCN)=1 ELSE NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) IF (IRAM(I2+1).LT.0) THEN JCN(1,NJCN)=-IRAM(I2+1) ELSE JCN(1,NJCN)=IRAM(I2+1)+ISHIFT ENDIF IF (IRAM(I2+2).LT.0) THEN JCN(2,NJCN)=-IRAM(I2+2) ELSE JCN(2,NJCN)=IRAM(I2+2)+ISHIFT ENDIF JCN(3,NJCN)=1 ENDIF 110 CONTINUE C Connections along the gridface: C Indices of gridlegs of the gridface: CALL MSGFAC(I1,IL1,IL2,IL3,IL4) C CALL MSGLEG(IL1,IP1,IP2) NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=NPOIEP+IP1*NQPOI DO 112, I2=IRAM(OLEG+IL1-1)+NQPOI,IRAM(OLEG+IL1),NQPOI ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4) ICBNEW=IRAM(I2 -NQPOI+4) IF (ICBOLD.NE.ICBNEW) THEN JCN(2,NJCN)=I2+ISHIFT JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2 ELSE JCN(2,NJCN)=I2 JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2+ISHIFT ENDIF 112 CONTINUE JCN(2,NJCN)=NPOIEP+IP2*NQPOI JCN(3,NJCN)=1 C CALL MSGLEG(IL2,IP1,IP2) NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=NPOIEP+IP1*NQPOI DO 114, I2=IRAM(OLEG+IL2-1)+NQPOI,IRAM(OLEG+IL2),NQPOI ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4) ICBNEW=IRAM(I2 -NQPOI+4) IF (ICBOLD.NE.ICBNEW) THEN JCN(2,NJCN)=I2+ISHIFT JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2 ELSE JCN(2,NJCN)=I2 JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2+ISHIFT ENDIF 114 CONTINUE JCN(2,NJCN)=NPOIEP+IP2*NQPOI JCN(3,NJCN)=1 C CALL MSGLEG(IL3,IP2,IP1) NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=NPOIEP+IP1*NQPOI DO 116, I2=IRAM(OLEG+IL3),IRAM(OLEG+IL3-1)+NQPOI,-NQPOI ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4) ICBNEW=IRAM(I2 -NQPOI+4) IF (ICBOLD.NE.ICBNEW) THEN JCN(2,NJCN)=I2+ISHIFT JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2 ELSE JCN(2,NJCN)=I2 JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2+ISHIFT ENDIF 116 CONTINUE JCN(2,NJCN)=NPOIEP+IP2*NQPOI JCN(3,NJCN)=1 C CALL MSGLEG(IL4,IP2,IP1) NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=NPOIEP+IP1*NQPOI DO 118, I2=IRAM(OLEG+IL4),IRAM(OLEG+IL4-1)+NQPOI,-NQPOI ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4) ICBNEW=IRAM(I2 -NQPOI+4) IF (ICBOLD.NE.ICBNEW) THEN JCN(2,NJCN)=I2+ISHIFT JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2 ELSE JCN(2,NJCN)=I2 JCN(3,NJCN)=1 NJCN=NJCN+1 IF (NJCN.GT.MJCN) CALL MSEROR(4) JCN(1,NJCN)=I2+ISHIFT ENDIF 118 CONTINUE JCN(2,NJCN)=NPOIEP+IP2*NQPOI JCN(3,NJCN)=1 C Array with connections is done. C C Forming polygons along the gridface: 124 CONTINUE DO 130, I2=1,NJCN IF (JCN(3,I2).NE.0) THEN C Initiating polygon: NLJPOL=2 JPOL(1,1)=JCN(1,I2) JPOL(2,1)=JCN(2,I2) JCN(3,I2)=0 126 CONTINUE DO 128, I3=1,NJCN IF ((JCN(3,I3).NE.0).AND. * (JCN(1,I3).EQ.JPOL(NLJPOL,1))) THEN NLJPOL=NLJPOL+1 IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7) JPOL(NLJPOL,1)=JCN(2,I3) JCN(3,I3)=0 IF (JPOL(NLJPOL,1).EQ.JPOL(1,1)) THEN C Polygon is closed - recording the polygon to IRAM: NLJPOL=NLJPOL-1 IF (NPOL-NLJPOL-1.LE.NCON) CALL MSEROR(1) NPOL=NPOL-1 IRAM(NPOL)=NLJPOL DO 127, I4=1,NLJPOL NPOL=NPOL-1 IRAM(NPOL)=JPOL(I4,1) 127 CONTINUE GOTO 124 ELSE GOTO 126 ENDIF ENDIF 128 CONTINUE C MODSRF-14 CALL ERROR ('MODSRF-14: Polygon not closed.') C This error should not appear. Contact the author. ENDIF 130 CONTINUE C End of the loop along gridfaces of the 2-D slice. 200 CONTINUE C C Preparing the normal to the 2-D slice: F(2)=0. F(3)=0. F(4)=0. IF (N1.EQ.1) F(2)=1. IF (N2.EQ.1) F(3)=1. IF (N3.EQ.1) F(4)=1. C ENDIF C C C Storing points on structural interfaces and polygons: WRITE(*,'(A)') *'+MODSRF: Writing. ' C C Storing points: OPEN(LU,FILE=VRTX,FORM='FORMATTED') WRITE(LU,'(A)') '/' NPOINT=(NPOI-OPOI)/NQPOI C Length of the names of points: LEN1=LENGTH(TEXTP) LEN2=0 202 CONTINUE IF (NPOINT.GE.10**LEN2) THEN LEN2=LEN2+1 GOTO 202 ENDIF LENG=LEN1+LEN2 C Loop over points: DO 218, I1=1,NPOINT C Address in RAM just before the current point: J1=OPOI+(I1-1)*NQPOI C Name of the point: DO 204, I2=0,LEN2-1 TEXTP(LENG-I2:LENG-I2)= * CHAR(ICHAR('0')+MOD(I1,10**(I2+1))/10**I2) 204 CONTINUE C Preparing the indices of complex blocks: IF (IRAM(J1+6).LT.0) THEN ICBPOS=IRAM(J1+8) ICBNEG=IRAM(J1+4) ELSE ICBPOS=IRAM(J1+4) ICBNEG=IRAM(J1+8) ENDIF C Preparing the normal to the interface at the point: IF (NCUB.NE.0) CALL SRFC2(IRAM(J1+6),RAM(J1+1),F) C The quantities at the point: J2=0 C Loop over the quantities to be stored: 206 CONTINUE IF ((TEXTC(J2+1).NE.' ').AND.(J2+1.LE.69)) THEN J2=J2+1 IF (TEXTC(J2).EQ.'x1') THEN C First coordinate of the point: IF (POWER(J2).EQ.1.) THEN Z(J2)=RAM(J1+1) ELSE Z(J2)=RAM(J1+1)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'x2') THEN C Second coordinate of the point: IF (POWER(J2).EQ.1.) THEN Z(J2)=RAM(J1+2) ELSE Z(J2)=RAM(J1+2)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'x3') THEN C Third coordinate of the point: IF (POWER(J2).EQ.1.) THEN Z(J2)=RAM(J1+3) ELSE Z(J2)=RAM(J1+3)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'norm1') THEN C First component of the normal: IF (POWER(J2).EQ.1.) THEN Z(J2)=F(2) ELSE Z(J2)=F(2)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'norm2') THEN C Second component of the normal: IF (POWER(J2).EQ.1.) THEN Z(J2)=F(3) ELSE Z(J2)=F(3)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'norm3') THEN C Third component of the normal: IF (POWER(J2).EQ.1.) THEN Z(J2)=F(4) ELSE Z(J2)=F(4)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'isrf') THEN Z(J2)=FLOAT(IABS(IRAM(J1+6)))**POWER(J2) ELSEIF (TEXTC(J2).EQ.'+icb') THEN Z(J2)=FLOAT(ICBPOS)**POWER(J2) ELSEIF (TEXTC(J2).EQ.'-icb') THEN Z(J2)=FLOAT(ICBNEG)**POWER(J2) ELSEIF (TEXTC(J2).EQ.'+isb') THEN IF (IRAM(J1+6).LT.0) THEN Z(J2)=FLOAT(IRAM(J1+7))**POWER(J2) ELSE Z(J2)=FLOAT(IRAM(J1+5))**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'-isb') THEN IF (IRAM(J1+6).LT.0) THEN Z(J2)=FLOAT(IRAM(J1+5))**POWER(J2) ELSE Z(J2)=FLOAT(IRAM(J1+7))**POWER(J2) ENDIF ELSEIF ((TEXTC(J2)(2:2).EQ.'v') * .OR.(TEXTC(J2)(2:2).EQ.'d') * .OR.(TEXTC(J2)(2:3).EQ.'qp') * .OR.(TEXTC(J2)(2:3).EQ.'qs')) THEN VP(1)=0. VS(1)=0. RHO=0. QP=0. QS=0. IF (TEXTC(J2)(1:1).EQ.'+') THEN IF (ICBPOS.NE.0) CALL PARM2(ICBPOS,Z(1),VP,VS,RHO,QP,QS) ELSE IF (ICBNEG.NE.0) CALL PARM2(ICBNEG,Z(1),VP,VS,RHO,QP,QS) ENDIF IF (TEXTC(J2)(2:3).EQ.'vp') THEN Z(J2)=VP(1)**POWER(J2) ELSEIF (TEXTC(J2)(2:3).EQ.'vs') THEN Z(J2)=VS(1)**POWER(J2) ELSEIF (TEXTC(J2)(2:3).EQ.'de') THEN Z(J2)=RHO**POWER(J2) ELSEIF (TEXTC(J2)(2:3).EQ.'qp') THEN Z(J2)=QP**POWER(J2) ELSEIF (TEXTC(J2)(2:3).EQ.'qs') THEN Z(J2)=QS**POWER(J2) ENDIF ELSEIF ((TEXTC(J2)(2:2).EQ.'a') * .OR.(TEXTC(J2)(2:2).EQ.'q')) THEN DO 208, I3=1,21 A(1,I3)=0. Q(1,I3)=0. 208 CONTINUE IF (TEXTC(J2)(1:1).EQ.'+') THEN IF (ICBPOS.NE.0) CALL PARM3(ICBPOS,Z(1),A,RHO,Q) ELSE IF (ICBNEG.NE.0) CALL PARM3(ICBNEG,Z(1),A,RHO,Q) ENDIF IF (TEXTC(J2)(2:4).EQ.'a11') THEN Z(J2)=A(1, 1)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a12') THEN Z(J2)=A(1, 2)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a22') THEN Z(J2)=A(1, 3)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a13') THEN Z(J2)=A(1, 4)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a23') THEN Z(J2)=A(1, 5)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a33') THEN Z(J2)=A(1, 6)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a14') THEN Z(J2)=A(1, 7)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a24') THEN Z(J2)=A(1, 8)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a34') THEN Z(J2)=A(1, 9)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a44') THEN Z(J2)=A(1,10)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a15') THEN Z(J2)=A(1,11)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a25') THEN Z(J2)=A(1,12)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a35') THEN Z(J2)=A(1,13)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a45') THEN Z(J2)=A(1,14)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a55') THEN Z(J2)=A(1,15)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a16') THEN Z(J2)=A(1,16)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a26') THEN Z(J2)=A(1,17)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a36') THEN Z(J2)=A(1,18)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a46') THEN Z(J2)=A(1,19)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a56') THEN Z(J2)=A(1,20)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'a66') THEN Z(J2)=A(1,21)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q11') THEN Z(J2)=Q(1, 1)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q12') THEN Z(J2)=Q(1, 2)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q22') THEN Z(J2)=Q(1, 3)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q13') THEN Z(J2)=Q(1, 4)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q23') THEN Z(J2)=Q(1, 5)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q33') THEN Z(J2)=Q(1, 6)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q14') THEN Z(J2)=Q(1, 7)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q24') THEN Z(J2)=Q(1, 8)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q34') THEN Z(J2)=Q(1, 9)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q44') THEN Z(J2)=Q(1,10)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q15') THEN Z(J2)=Q(1,11)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q25') THEN Z(J2)=Q(1,12)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q35') THEN Z(J2)=Q(1,13)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q45') THEN Z(J2)=Q(1,14)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q55') THEN Z(J2)=Q(1,15)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q16') THEN Z(J2)=Q(1,16)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q26') THEN Z(J2)=Q(1,17)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q36') THEN Z(J2)=Q(1,18)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q46') THEN Z(J2)=Q(1,19)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q56') THEN Z(J2)=Q(1,20)**POWER(J2) ELSEIF (TEXTC(J2)(2:4).EQ.'q66') THEN Z(J2)=Q(1,21)**POWER(J2) ENDIF ELSEIF (TEXTC(J2).EQ.'srf') THEN CALL SRFC2(IVALUE(J2),RAM(J1+1),FF) IF (POWER(J2).EQ.1.) THEN Z(J2)=FF(1) ELSE Z(J2)=FF(1)**POWER(J2) ENDIF ENDIF GOTO 206 ENDIF C End of the loop for other quantities. C Writing the quantities at the point: C Setting output format for the array: FORMAT='(3A,00(F00.0,1X),A)' FORMAT(6:6)=CHAR(ICHAR('0')+MOD(J2,10)) FORMAT(5:5)=CHAR(ICHAR('0')+J2/10) OUTMIN=0. OUTMAX=0. DO 214, I2=1,J2 IF (OUTMIN.GT.Z(I2)) OUTMIN=Z(I2) IF (OUTMAX.LT.Z(I2)) OUTMAX=Z(I2) 214 CONTINUE CALL FORM1(OUTMIN,OUTMAX,FORMAT(8:15)) FORMAT(14:17)= '1X),' C Output format is set. WRITE(LU,FORMAT) '''',TEXTP(1:LENG),''' ',(Z(I2),I2=1,J2),'/' 218 CONTINUE C End of the loop over all points. WRITE(LU,'(A)') '/' CLOSE(LU) C C Storing polygons: IF(PLGN.NE.' ') OPEN(LU1,FILE=PLGN,FORM='FORMATTED') IF(PLGNS.NE.' ')OPEN(LU2,FILE=PLGNS,FORM='FORMATTED') FORMA1='(00(I8,1X),A)' FORMA2='(00(3A,1X),A)' I1=MRAM C Loop over polygons: 220 CONTINUE IF (I1.GT.NPOL) THEN J1=IRAM(I1) I3=MOD(J1,10) FORMA1(3:3)=CHAR(ICHAR('0')+I3) FORMA2(3:3)=CHAR(ICHAR('0')+I3) I3=J1/10 FORMA1(2:2)=CHAR(ICHAR('0')+I3) FORMA2(2:2)=CHAR(ICHAR('0')+I3) J2=IRAM(I1-1) IF ((IRAM(J2-NQPOI+4).NE.0).OR.(IRAM(J2-NQPOI+8).NE.0).OR. * (LFREE)) THEN IF (PLGN.NE.' ') THEN WRITE(LU1,FORMA1) * ((IRAM(I2)-OPOI)/NQPOI,I2=I1-1,I1-J1,-1),'/' ENDIF IF (PLGNS.NE.' ') THEN IF (J1.GT.MLJPOL) THEN C MODSRF-15 CALL ERROR ('MODSRF-15: Disorder in polygons.') C This error should not appear. Contact the author. ENDIF DO 224, I2=I1-1,I1-J1,-1 C Index of the point (from 1 to NPOINT): J2=(IRAM(I2)-OPOI)/NQPOI C Name of the point: I4=I1-I2 DO 222, I3=0,LEN2-1 TEXTS(I4)(1:LEN1)=TEXTP(1:LEN1) TEXTS(I4)(LENG-I3:LENG-I3)= * CHAR(ICHAR('0')+MOD(J2,10**(I3+1))/10**I3) 222 CONTINUE 224 CONTINUE WRITE(LU2,FORMA2) * ('''',TEXTS(I2)(1:LENG),'''',I2=1,J1),'/' ENDIF ENDIF I1=I1-(J1+1) GOTO 220 ENDIF C End of the loop over polygons. IF (PLGN.NE.' ') THEN WRITE(LU1,'(A)') '/' CLOSE(LU1) ENDIF IF (PLGNS.NE.' ') THEN WRITE(LU2,'(A)') '/' CLOSE(LU2) ENDIF C WRITE(*,'(A)') *'+MODSRF: Done. ' STOP END C C======================================================================= C SUBROUTINE FCTMS(X,Y,T) C C----------------------------------------------------------------------- C REAL X,Y(*),T(*) C C Subroutine evaluating the right-hand sides of the interface tracing C equations. I.E. subroutine computing vector T tangent to the C interface. C C Input: C X... Value of the independent variable. C Y... Array containing two coordinates of a C point of the interface, determined by means of numerical C integration. C Output: C Y... Array containing two coordinates of the C interface, corrected by means of the linearization in C the direction of the gradient (perpendicular to the C interface). C T... Array containing derivatives of the coordinates C Y with respect to X (vector tangent to the interface). C----------------------------------------------------------------------- C C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: REAL S1,S2,S3,F(10),AUX C C....................................................................... C CALL SRFC2(ISRFO,Y,F) IF (ISRFO.LT.0) THEN S1=F(2) S2=F(3) S3=F(4) ELSE S1=-F(2) S2=-F(3) S3=-F(4) ENDIF IF (IFACE.LE.NFAC1) THEN F(2)=0. S1=0. T(1)=0. AUX=SQRT(S2*S2+S3*S3) T(2)=-S3/AUX T(3)= S2/AUX ELSEIF (IFACE.LE.NFAC12) THEN F(3)=0. S2=0. T(2)=0. AUX=SQRT(S1*S1+S3*S3) T(3)=-S1/AUX T(1)= S3/AUX ELSE F(4)=0. S3=0. T(3)=0. AUX=SQRT(S1*S1+S2*S2) T(1)=-S2/AUX T(2)= S1/AUX ENDIF C C Correction of the isoline AUX=F(1)/AUX/AUX Y(1)=Y(1)-F(2)*AUX Y(2)=Y(2)-F(3)*AUX Y(3)=Y(3)-F(4)*AUX C RETURN END C C======================================================================= C SUBROUTINE OUTMS(X,Y,DERY,IHLF,NDIM,PRMT) C C----------------------------------------------------------------------- C INTEGER IHLF,NDIM REAL X,Y(NDIM),DERY(NDIM),PRMT(*) C C Subroutine to test, whether an interface between two complex blocks C was crossed between the old point (stored from the previous invocation C of OUTMS), and the new point. The new point along the interface C being traced was computed by RKGS. C If the interface is crossed, the indices of simple blocks, C complex blocks, and surfaces ISBO1,ICBO1,ISBO2,ICBO2,ISRFN1,ISRFN2, C ISBN1,ICBN1,ISBN2,ICBN2 are computed and stored in common block MSC. C C Input: C X... Value of the independent variable in the new point. C Y... Array containing the coordinates of the new C point of the interface, determined by means of numerical C integration. C DERY... Array containing derivatives of the coordinates C Y with respect to X (vector tangent to the interface). C IHLF... Number of bisections of the initial increment PRMT(3). C NDIM... Dimension of arrays, should be 3. C PRMT(1)... Lower bound of the interval for the independent C variable. C PRMT(2)... Upper bound of the interval for the independent C variable. C PRMT(3)... Initial increment of the independent variable. C PRMT(4)... Upper error bound. C Output: C If any interface was crossed C and the new point lies within the considered gridface: C No output C If any interface was crossed C and the new point lies outside the considered gridface: C PRMT(5)=2. C If an interface was crossed C and the point of intersection lies within the considered gridface: C PRMT(5)=1. C Y... Array containing the coordinates of the point C of the intersection of the traced and the crossed C interface. C ISBO1,ICBO1,ISBO2,ICBO2,ISRFN1,ISRFN2,ISBN1,ICBN1,ISBN2,ICBN2... C The indices of simple blocks, complex blocks, and surfaces C are computed and stored in common block MSC. C----------------------------------------------------------------------- C C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... EXTERNAL ISIDE,MSLPIF INTEGER ISIDE LOGICAL MSLPIF C Auxiliary storage locations: REAL XTMP1,XTMP(3),DTMP(3),DUMMY(1),ERRSRF,XINT1,DINT(3) REAL X01,X02,X1(3),X2(3) INTEGER IDUMMY,IY1(8),IY2(8),ISUR(2) C....................................................................... IF (IHLF.GT.10) THEN C MODSRF-16 CALL ERROR ('MODSRF-16: Wrong IHLF.') C This error should not appear. Contact the author. ENDIF IF (X.GT.PRMT(1)) THEN X01=0. X02=0. C Calling CDE along side 1: ERRSRF=PRMT(4) XTMP1=X XTMP(1)=Y(1) XTMP(2)=Y(2) XTMP(3)=Y(3) DTMP(1)=DERY(1) DTMP(2)=DERY(2) DTMP(3)=DERY(3) IY1(4)=ISBO1 IY1(5)=ICBO1 IY1(6)=0 CALL CDE(ISRFO,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY1,ERRSRF, * X0,X0,YOLD,DYOLD,X,Y,DERY,XTMP1,XTMP,DTMP,XINT1,X1,DINT) IF (IY1(6).NE.0) THEN C Interface crossed: X01=XTMP1 X1(1)=XTMP(1) X1(2)=XTMP(2) X1(3)=XTMP(3) ENDIF C C Calling CDE along side 2: XTMP1=X XTMP(1)=Y(1) XTMP(2)=Y(2) XTMP(3)=Y(3) DTMP(1)=DERY(1) DTMP(2)=DERY(2) DTMP(3)=DERY(3) IY2(4)=ISBO2 IY2(5)=ICBO2 IY2(6)=0 CALL CDE(ISRFO,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY2,ERRSRF, * X0,X0,YOLD,DYOLD,X,Y,DERY,XTMP1,XTMP,DTMP,XINT1,X2,DINT) IF (IY2(6).NE.0) THEN C Interface crossed: X02=XTMP1 X2(1)=XTMP(1) X2(2)=XTMP(2) X2(3)=XTMP(3) ENDIF C IF ((X01.NE.0.).OR.(X02.NE.0.)) THEN C Interface crossed. C IF ((X01.NE.0.).AND.(X02.NE.0.)) THEN C Choosing the nearer of the two crossed interfaces: IF (X01.LT.X02) THEN X02=0. ELSEIF (X02.LT.X01) THEN X01=0. ENDIF ENDIF C C Recording indices of blocks computed by CDE: IF (X01.NE.0.) THEN ISBO1= IY1(4) ISRFN1=IY1(6) ISBN1= IY1(7) ICBN1= IY1(8) Y(1)=X1(1) Y(2)=X1(2) Y(3)=X1(3) ENDIF IF (X02.NE.0.) THEN ISBO2= IY2(4) ISRFN2=IY2(6) ISBN2= IY2(7) ICBN2= IY2(8) Y(1)=X2(1) Y(2)=X2(2) Y(3)=X2(3) ENDIF C Computing remaining indices of blocks: IF (X01.EQ.0.) THEN ISUR(1)=ISRFO ISUR(2)=ISRFN2 CALL BLOCKS(X2,2,ISUR,ISBO2,IDUMMY,ISBO1,ICBO1) ISRFN1=ISRFN2 IF (ISIDE(ISRFO,ISBN2).LT.2) THEN CALL BLOCKS(X2,2,ISUR,ISBN2,IDUMMY,ISBN1,ICBN1) ELSE ISBN1=ISBN2 ICBN1=ICBN2 ENDIF ENDIF IF (X02.EQ.0.) THEN ISUR(1)=ISRFO ISUR(2)=ISRFN1 CALL BLOCKS(X1,2,ISUR,ISBO1,IDUMMY,ISBO2,ICBO2) ISRFN2=ISRFN1 IF (ISIDE(ISRFO,ISBN1).LT.2) THEN CALL BLOCKS(X1,2,ISUR,ISBN1,IDUMMY,ISBN2,ICBN2) ELSE ISBN2=ISBN1 ICBN2=ICBN1 ENDIF ENDIF c IF (MSLPIF(Y,IFACE)) THEN C Point of intersection lies within the gridface. PRMT(5)=1. c ELSE cC Point of intersection lies outside the gridface. c PRMT(5)=2. c ENDIF ELSE C Interface not crossed. IF (.NOT.MSLPIF(Y,IFACE)) THEN C New point Y is outside the gridface, C RKGS is to be terminated. PRMT(5)=2. ENDIF ENDIF ENDIF C C Storing the old point: X0=X YOLD(1)=Y(1) YOLD(2)=Y(2) YOLD(3)=Y(3) DYOLD(1)=DERY(1) DYOLD(2)=DERY(2) DYOLD(3)=DERY(3) RETURN END C C C======================================================================= C SUBROUTINE MSGLEG(ILEG,IPOIN1,IPOIN2) C C----------------------------------------------------------------------- C INTEGER ILEG,IPOIN1,IPOIN2 C Input: C ILEG... Index of gridleg. C Output: C IPOIN1,IPOIN2... Indices of gridpoints forming the gridleg. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,ILEG1 C....................................................................... C ILEG1=ILEG-1 IF (ILEG.LE.NLEG1) THEN I31=ILEG1 / N11N2 I21=(ILEG1 - I31*N11N2) / N11 I1=ILEG - I31*N11N2 - I21*N11 IPOIN1=I31*N1N2+I21*N1+I1 IPOIN2=IPOIN1+1 ELSEIF (ILEG.LE.NLEG12) THEN I31=(ILEG1-NLEG1) / N1N21 I21=((ILEG1-NLEG1) - I31*N1N21) / N1 I1=(ILEG-NLEG1) - I31*N1*N21 - I21*N1 IPOIN1=I31*N1N2+I21*N1+I1 IPOIN2=IPOIN1+N1 ELSEIF (ILEG.LE.NLEG) THEN I31=(ILEG1-NLEG12) / N1N2 I21=((ILEG1-NLEG12) - I31*N1N2) / N1 I1=(ILEG-NLEG12) - I31*N1N2 - I21*N1 IPOIN1=I31*N1N2+I21*N1+I1 IPOIN2=IPOIN1+N1*N2 ELSE C MODSRF-17 CALL ERROR ('MODSRF-17: Wrong ILEG.') C This error should not appear. Contact the author. ENDIF RETURN END C C C======================================================================= C SUBROUTINE MSGFAC(IFAC,ILEG1,ILEG2,ILEG3,ILEG4) C C----------------------------------------------------------------------- C INTEGER IFAC,ILEG1,ILEG2,ILEG3,ILEG4 C IFAC... Index of gridface. C ILEG1,ILEG2,ILEG3,ILEG4... Indices of gridlegs forming C the gridface. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,IFAC1 C....................................................................... C IFAC1=IFAC-1 IF (IFAC.LE.NFAC1) THEN I31= IFAC1 / N1N21 I21=(IFAC1 - I31*N1N21) / N1 I1 = IFAC - I31*N1N21 - I21*N1 ILEG1=NLEG1 + I31*N1N21+I21*N1+I1 ILEG4=NLEG12 + I31*N1N2+I21*N1+I1 ILEG2=ILEG4+N1 ILEG3=ILEG1+N1N21 ELSEIF (IFAC.LE.NFAC12) THEN I31= (IFAC1-NFAC1) / N11N2 I21=((IFAC1-NFAC1) - I31*N11N2) / N11 I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11 ILEG1=NLEG12 + I31*N1N2+I21*N1+I1 ILEG4= I31*N11N2+I21*N11+I1 ILEG2=ILEG4+N11N2 ILEG3=ILEG1+1 ELSEIF (IFAC.LE.NFAC) THEN I31= (IFAC1-NFAC12) / N11N21 I21=((IFAC1-NFAC12) - I31*N11N21) / N11 I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11 ILEG1= I31*N11N2+I21*N11+I1 ILEG4=NLEG1 + I31*N1N21+I21*N1+I1 ILEG2=ILEG4+1 ILEG3=ILEG1+N11 ELSE C MODSRF-18 CALL ERROR ('MODSRF-18: Wrong IFAC.') C This error should not appear. Contact the author. ENDIF RETURN END C C C======================================================================= C SUBROUTINE MSGCUB(ICUB,IFAC1,IFAC2,IFAC3,IFAC4,IFAC5,IFAC6) C C----------------------------------------------------------------------- C INTEGER ICUB,IFAC1,IFAC2,IFAC3,IFAC4,IFAC5,IFAC6 C ICUB... Index of gridcube. C IFAC1,...,IFAC6... Indices of gridfaces forming the gridcube. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,ICUB1 C....................................................................... C IF ((ICUB.LT.1).OR.(ICUB.GT.NCUB)) THEN C MODSRF-19 CALL ERROR ('MODSRF-19: Wrong ICUB.') C This error should not appear. Contact the author. ENDIF ICUB1=ICUB-1 I31= ICUB1 / N11N21 I21=(ICUB1 - I31*N11N21) / N11 I1= ICUB - I31*N11N21 - I21*N11 IFAC1= I31*N1N21 +I21*N1 +I1 IFAC2=NFAC1 +I31*N11N2 +I21*N11+I1 IFAC3=NFAC12+I31*N11N21+I21*N11+I1 IFAC4=IFAC1+1 IFAC5=IFAC2+N11 IFAC6=IFAC3+N11N21 RETURN END C C C======================================================================= C LOGICAL FUNCTION MSLPIF(XX,IFAC) C C----------------------------------------------------------------------- C Subroutine for decision, whether the point with coordinates XX C lies inside gridface IFAC. C INTEGER IFAC REAL XX(3) C Input: C XX... Coordinates of the point. C IFAC... Index of the gridface. C Output: C MSLPIF... .TRUE. when the point lies inside the gridface. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,IFAC1,IP1,IP2 INTEGER ILEG1,ILEG2 REAL XA1,XA2,XA3,XB1,XB2,XB3,XC1,XC2,XC3 C....................................................................... C IFAC1=IFAC-1 IF (IFAC.LE.NFAC1) THEN I31= IFAC1 / N1N21 I21=(IFAC1 - I31*N1N21) / N1 I1 = IFAC - I31*N1N21 - I21*N1 ILEG1=NLEG1 + I31*N1N21+I21*N1+I1 ILEG2=NLEG12 + I31*N1N2+I21*N1+I1 + N1 ELSEIF (IFAC.LE.NFAC12) THEN I31= (IFAC1-NFAC1) / N11N2 I21=((IFAC1-NFAC1) - I31*N11N2) / N11 I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11 ILEG1=NLEG12 + I31*N1N2+I21*N1+I1 ILEG2= I31*N11N2+I21*N11+I1 + N11N2 ELSEIF (IFAC.LE.NFAC) THEN I31= (IFAC1-NFAC12) / N11N21 I21=((IFAC1-NFAC12) - I31*N11N21) / N11 I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11 ILEG1= I31*N11N2+I21*N11+I1 ILEG2=NLEG1 + I31*N1N21+I21*N1+I1 + 1 ELSE C MODSRF-20 CALL ERROR ('MODSRF-20: Wrong IFAC.') C This error should not appear. Contact the author. ENDIF C MSLPIF=.TRUE. C CALL MSGLEG(ILEG1,IP1,IP2) CALL MSCP(IP1,XA1,XA2,XA3) CALL MSCP(IP2,XB1,XB2,XB3) CALL MSGLEG(ILEG2,IP1,IP2) CALL MSCP(IP2,XC1,XC2,XC3) IF ((XX(1).LT.AMIN1(XA1,XB1,XC1)).OR. * (XX(1).GT.AMAX1(XA1,XB1,XC1))) MSLPIF=.FALSE. IF ((XX(2).LT.AMIN1(XA2,XB2,XC2)).OR. * (XX(2).GT.AMAX1(XA2,XB2,XC2))) MSLPIF=.FALSE. IF ((XX(3).LT.AMIN1(XA3,XB3,XC3)).OR. * (XX(3).GT.AMAX1(XA3,XB3,XC3))) MSLPIF=.FALSE. RETURN END C C C======================================================================= C LOGICAL FUNCTION MSLPIL(XX,ILEG) C C----------------------------------------------------------------------- C Subroutine for decision, whether the point with coordinates XX C lies inside gridleg ILEG. C INTEGER ILEG REAL XX(3) C Input: C XX... Coordinates of the point. C ILEG... Index of the gridleg. C Output: C MSLPIL... .TRUE. when the point lies inside the gridleg. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER IP1,IP2 REAL XA1,XA2,XA3,XB1,XB2,XB3 C....................................................................... C CALL MSGLEG(ILEG,IP1,IP2) CALL MSCP(IP1,XA1,XA2,XA3) CALL MSCP(IP2,XB1,XB2,XB3) MSLPIL=.TRUE. IF ((XX(1).LT.AMIN1(XA1,XB1)).OR. * (XX(1).GT.AMAX1(XA1,XB1))) MSLPIL=.FALSE. IF ((XX(2).LT.AMIN1(XA2,XB2)).OR. * (XX(2).GT.AMAX1(XA2,XB2))) MSLPIL=.FALSE. IF ((XX(3).LT.AMIN1(XA3,XB3)).OR. * (XX(3).GT.AMAX1(XA3,XB3))) MSLPIL=.FALSE. RETURN END C C C======================================================================= C SUBROUTINE MSGCGP(ICUB,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8) C C----------------------------------------------------------------------- C Subroutine for debugging purposes, gives indices of points C of cube ICUB. C C INTEGER ICUB,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8 C ICUB... Index of gridcube. C IP1,...,IP8... Indices of gridpoints forming the gridcube. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,ICUB1 C....................................................................... C IF ((ICUB.LT.1).OR.(ICUB.GT.NCUB)) THEN C MODSRF-21 CALL ERROR ('MODSRF-21: Wrong ICUB.') C This error should not appear. Contact the author. ENDIF ICUB1=ICUB-1 I31= ICUB1 / N11N21 I21=(ICUB1 - I31*N11N21) / N11 I1= ICUB - I31*N11N21 - I21*N11 IP1=I31*N1N2+I21*N1+I1 IP2=IP1+N1 IP3=IP2+N1N2 IP4=IP1+N1N2 IP5=IP1+1 IP6=IP5+N1 IP7=IP6+N1N2 IP8=IP5+N1N2 RETURN END C C C======================================================================= C SUBROUTINE MSXFAC(IFAC) C C----------------------------------------------------------------------- C Subroutine for debugging purposes, gives coordinates of points C of face IFAC. C INTEGER IFAC,ILEG1,ILEG2,ILEG3,ILEG4 C IFAC... Index of gridface. C ILEG1,ILEG2,ILEG3,ILEG4... Indices of gridlegs forming C the gridface. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I1,IFAC1,IP1,IP2 REAL X1,X2,X3 C....................................................................... C IFAC1=IFAC-1 IF (IFAC.LE.NFAC1) THEN I31= IFAC1 / N1N21 I21=(IFAC1 - I31*N1N21) / N1 I1 = IFAC - I31*N1N21 - I21*N1 ILEG1=NLEG1 + I31*N1N21+I21*N1+I1 ILEG4=NLEG12 + I31*N1N2+I21*N1+I1 ILEG2=ILEG4+N1 ILEG3=ILEG1+N1N21 ELSEIF (IFAC.LE.NFAC12) THEN I31= (IFAC1-NFAC1) / N11N2 I21=((IFAC1-NFAC1) - I31*N11N2) / N11 I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11 ILEG1=NLEG12 + I31*N1N2+I21*N1+I1 ILEG4= I31*N11N2+I21*N11+I1 ILEG2=ILEG4+N11N2 ILEG3=ILEG1+1 ELSEIF (IFAC.LE.NFAC) THEN I31= (IFAC1-NFAC12) / N11N21 I21=((IFAC1-NFAC12) - I31*N11N21) / N11 I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11 ILEG1= I31*N11N2+I21*N11+I1 ILEG4=NLEG1 + I31*N1N21+I21*N1+I1 ILEG2=ILEG4+1 ILEG3=ILEG1+N11 ELSE C MODSRF-22 CALL ERROR ('MODSRF-22: Wrong IFAC.') C This error should not appear. Contact the author. ENDIF CALL MSGLEG(ILEG1,IP1,IP2) CALL MSCP(IP1,X1,X2,X3) CALL MSCP(IP2,X1,X2,X3) C CALL MSGLEG(ILEG2,IP1,IP2) CALL MSCP(IP1,X1,X2,X3) CALL MSCP(IP2,X1,X2,X3) C CALL MSGLEG(ILEG3,IP1,IP2) CALL MSCP(IP1,X1,X2,X3) CALL MSCP(IP2,X1,X2,X3) C CALL MSGLEG(ILEG4,IP1,IP2) CALL MSCP(IP1,X1,X2,X3) CALL MSCP(IP2,X1,X2,X3) C RETURN END C C C======================================================================= C SUBROUTINE MSCP(IPOIN,X1,X2,X3) C C----------------------------------------------------------------------- C INTEGER IPOIN REAL X1,X2,X3 C IPOIN... Index of a point. C X1,X2,X3... Coordinates of the point. C....................................................................... C Common block /RAMC/ to store the information about points on C structural interfaces and common block /MSC/ to store auxiliary C quantities: INCLUDE 'modsrf.inc' C modsrf.inc. C ........................... C Auxiliary storage locations: INTEGER I31,I21,I11,IPOIN1 C....................................................................... C IPOIN1=IPOIN-1 I31= IPOIN1 / N1N2 I21=(IPOIN1 - I31*N1N2) / N1 I11= IPOIN - I31*N1N2 - I21*N1 - 1 X1=O1+FLOAT(I11)*D1 X2=O2+FLOAT(I21)*D2 X3=O3+FLOAT(I31)*D3 RETURN END C C======================================================================= C SUBROUTINE MSJCN(JCN,NJCN,MJCN) C C----------------------------------------------------------------------- C Subroutine for marking sure and wrong connections in array JCN. C INTEGER MJCN,JCN(4,MJCN),NJCN C JCN(1 to 4,i)... Connection i of the face: address of the first C point, address of the second point, index of the C interface, status of the connection. C Auxiliary storage locations: INTEGER I4,I3,I2 C....................................................................... C 10 CONTINUE I4=0 DO 60, I2=1,NJCN C Loop over unsure connections: IF (JCN(4,I2).EQ.0) THEN DO 20, I3=1,NJCN C Loop over other connections with the same starting point: IF ((I3.NE.I2).AND.(JCN(1,I3).EQ.JCN(1,I2))) THEN IF (JCN(4,I3).EQ.1) THEN C Another connection is sure, connection I2 is wrong. JCN(4,I2)=-1 GOTO 50 ELSEIF (JCN(4,I3).EQ.0) THEN C Another unsure connection, connection I2 is unsure. GOTO 30 ENDIF ENDIF 20 CONTINUE C No other connection, connection I2 is sure. JCN(4,I2)=1 GOTO 50 30 CONTINUE DO 40, I3=1,NJCN C Loop over other connections with the same end point: IF ((I3.NE.I2).AND.(JCN(2,I3).EQ.JCN(2,I2))) THEN IF (JCN(4,I3).EQ.1) THEN C Another connection is sure, connection I2 is wrong. JCN(4,I2)=-1 GOTO 50 ELSEIF (JCN(4,I3).EQ.0) THEN C Another unsure connection, connection I2 is unsure. GOTO 50 ENDIF ENDIF 40 CONTINUE C No other connection, connection I2 is sure. JCN(4,I2)=1 50 CONTINUE C Noting the change in the status of the connection: IF (JCN(4,I2).NE.0) I4=1 ENDIF 60 CONTINUE C Change in the status of a connection implies C the repetition of the loop: IF (I4.EQ.1) GOTO 10 RETURN END C C======================================================================= C SUBROUTINE MSEROR(IERR) C C----------------------------------------------------------------------- C INTEGER IERR C IERR... Index of the error. C....................................................................... IF (IERR.EQ.1) THEN C MODSRF-23 CALL ERROR('MODSRF-23: Small array (I)RAM.') C Try to enlarge the dimension MRAM in the file C ram.inc. C The memory requirements of this program are C roughly 8*N1*N2*N3+12*NPTS, where N1*N2*N3 C is the dimension of the input grid, and NPTS is the number of C intersections of gridlegs with structural interfaces. ELSEIF (IERR.EQ.2) THEN C MODSRF-24 CALL ERROR('MODSRF-24: Small array JPT.') C Try to enlarge the dimension MJPT at the beginning of this file. C If this does not help, contact the author. ELSEIF (IERR.EQ.3) THEN C MODSRF-25 CALL ERROR('MODSRF-25: Small array IPTE.') C Try to enlarge the dimension MIPTE at the beginning of this C file. If this does not help, contact the author. ELSEIF (IERR.EQ.4) THEN C MODSRF-26 CALL ERROR('MODSRF-26: Small array JCN.') C Try to enlarge the dimension MJCN at the beginning of this file. C If this does not help, contact the author. ELSEIF (IERR.EQ.6) THEN C MODSRF-27 CALL ERROR('MODSRF-27: Small array JPOL.') C Try to enlarge the dimension MJPOL at the beginning of this C file. If this does not help, contact the author. ELSEIF (IERR.EQ.7) THEN C MODSRF-28 CALL ERROR('MODSRF-28: Small array JPOL.') C Try to enlarge the dimension MLJPOL at the beginning of this C file. If this does not help, contact the author. ELSEIF (IERR.EQ.8) THEN C MODSRF-29 CALL ERROR('MODSRF-29: Interface not found.') C This might happen when the whole gridleg coincides with an C interface. If possible, slightly change the origin of C the grid. If this does not help, contact the author. ELSE C MODSRF-30 CALL ERROR('MODSRF-30: Wrong IERR.') C This subroutine was invocated with wrong value of IERR. C This error should not appear, please contact the author. ENDIF END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'model.for' C model.for INCLUDE 'metric.for' C metric.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'parm.for' C parm.for INCLUDE 'val.for' C val.for INCLUDE 'fit.for' C fit.for INCLUDE 'means.for' C means.for INCLUDE 'rkgs.for' C rkgs.for C C======================================================================= C