C SUBROUTINE FILE 'MODSEC.FOR' TO DETERMINE THE INTERFACES AND VELOCITY C ISOLINES IN 2-D SECTIONS OF A 3-D SEISMIC MODEL. C C BY LUDEK KLIMES, IVAN PSENCIK C C THIS FILE CONSISTS OF THE FOLLOWING EXTERNAL PROCEDURES: C SECTB...BLOCK DATA SUBROUTINE DEFINING COMMON BLOCKS /SECTC/ AND C /VALUES/ TO STORE THE INPUT DATA. C FUNC... SUBROUTINE DESIGNED TO EVALUATE THE VALUE, FIRST AND C SECOND DERIVATIVES OF THE GIVEN FUNCTION WITH RESPECT TO C THE PLOT COORDINATES. C DISC... SUBROUTINE DESIGNED TO DETERMINE THE POINT OF INTERSECTION C OF THE GIVEN LINE SEGMENT WITH THE BOUNDARY OF THE COMPLEX C BLOCK. IT MAY ALSO BE USED TO DETERMINE THE INDEX OF THE C BLOCK IN WHICH THE GIVEN POINT IS SITUATED. C ISOL... SUBROUTINE DESIGNED TO FIND THE POINT OF INTERSECTION C OF THE GIVEN LINE SEGMENT WITH AN ISOLINE. C ISOLA...AUXILIARY SUBROUTINE TO THE SUBROUTINE ISOL. C SECT1...SUBROUTINE DESIGNED TO READ THE INPUT DATA FOR THE C PLOTTED SECTION OF THE MODEL AND TO STORE THEM IN THE C MEMORY. C SECT2...AUXILIARY SUBROUTINE TO THE SUBROUTINES DISC AND ISOL. C THE SUBROUTINE TRANSFORMS THE PLOT COORDINATES TO THE C MODEL COORDINATES. C SECT3...AUXILIARY SUBROUTINE TO THE SUBROUTINES DISC AND ISOL. C THE SUBROUTINE TRANSFORMS THE MODEL COORDINATES TO THE C PLOT COORDINATES. C MODSEC..SUBROUTINE DEMONSTRATING THE FUNCTION OF THE SUBROUTINES C DISC AND ISOL (AND, PARTIALLY, ALSO FUNC). IT EMPLOYES C THE SUBROUTINES WHEN SKETCHING A SECTION OF THE MODEL BY C MEANS OF EXTENDED ASCII CHARACTERS ONTO THE SCREEN. THE C SUBROUTINE IS INTENDED TO BE JUST AN EXAMPLE HOW TO USE C THE SUBROUTINES FUNC, DISC AND ISOL OF THIS FILE. C CONT1...SUBROUTINE DESIGNED TO INITIALIZE ARRAYS CONTAINING THE C POINTS OF INTERSECTION OF ISOLINES WITH VERTICAL LINES C LIMITING THE REGIONS OF NUMERICALY TRACING THE ISOLINES. C IT IS CALLED ONCE BEFORE TRACING ISOLINES IN A NEW COLUMN C OF THE SECTION. C CONT2...SUBROUTINE DESIGNED TO TRACE AN ISOLINE BY MEANS OF C NUMERICAL INTEGRATION WITHIN THE GIVEN COLUMN. C FCTI... SUBROUTINE EVALUATING THE RIGHT-HAND SIDES OF THE ISOLINE C TRACING EQUATIONS. C OUTI... SUBROUTINE DESIGNED TO CHECK FOR THE INTERSECTIONS OF THE C ISOLINE WITH STRUCTURAL INTERFACES OR BOUNDARIES OF THE C COLUMN IN WHICH THE ISOLINE IS TRACED. C C TO DEMONSTRATE THE USAGE OF THIS SUBROUTINE FILE, REMOVE THE C ASTERISKS FROM THE FIRST COLUMN OF THE THREE FOLLOWING STATEMENTS AND C COMPILE: * CALL MODSEC * STOP * END C C THE OUTPUT INTERFACES AND ISOLINES MAY BE WRITTEN EITHER IN THE FORM C OF SHORT LINES (INTERSECTIONS OF INTERFACES AND ISOSURFACES WITH GIVEN C SECTIONS, STORED BY PARTS), OR IN THE FORM OF POINTS (INTERSECTIONS C OF INTERFACES AND ISOSURFACES WITH GIVEN LINES). C C FILENAME OF THE INPUT DATA FOR THE SPECIFICATION OF THE PLOTTED C SECTIONS OF THE MODEL, READ FROM THE * INTERACTIVE INPUT DEVICE: C (1) 'MODSEC' C 'MODSEC'... NAME OF THE FILE CONTAINING INPUT DATA FOR THE C SPECIFICATION OF THE PLOTTED SECTIONS. C DEFAULT: 'SEC.DAT' C C INPUT DATA FOR THE SPECIFICATION OF THE PLOTTED SECTION OF THE MODEL, C READ FROM THE ABOVE MENTIONED FILE 'MODSEC': C THE DATA ARE READ IN BY THE LIST DIRECTED INPUT (FREE FORMAT). IN C THE LIST OF INPUT DATA BELOW, EACH NUMBERED PARAGRAPH INDICATES C THE BEGINNING OF A NEW INPUT OPERATION (NEW READ STATEMENT). C THE VARIABLE NAMES ENCLOSED IN APOSTROPHES CORRESPOND TO CHARACTER C STRINGS. IF THE FIRST LETTER OF THE SYMBOLIC NAME OF THE INPUT C VARIABLE IS I-N, THE CORRESPONDING VALUE IN INPUT DATA MUST BE OF C THE TYPE INTEGER. OTHERWISE, THE INPUT PARAMETER IS OF THE TYPE C REAL. C (1) 'MODEL','SECTS' C 'MODEL'...STRING CONTAINING THE NAME OF THE FILE WITH THE INPUT C DATA FOR THE MODEL. THE DATA WILL BE READ IN BY THE C SUBROUTINE MODEL1. ONLY THE FIRST 80 CHARACTERS OF THE C STRING ARE SIGNIFICANT. C 'SECTS'...THE STRING CONTAINING THE NAME OF THE OUTPUT FILE WITH C THE GENERATED MODEL SECTIONS. ONLY THE FIRST 80 C CHARACTERS OF THE STRING ARE SIGNIFICANT. C DEFAULT: 'MODEL'='MODEL.DAT', 'SECTS'='SEC.OUT'. C (2) IPS,NC1,NC2,STEP,ERR C IPS... POSITIVE TO PLOT P-WAVE VELOCITY ISOLINES, C NEGATIVE TO PLOT S-WAVE VELOCITY ISOLINES. C NC1... NUMBER OF COLUMNS IN EACH SECTION. IN CARTESIAN C COORDINATES EACH SECTION HAS THE SHAPE OF RHOMBOID. C IN THE NORMALIZED SECTION COORDINATES, EACH SECTION IS A C SQUARE (0,1)*(0,1). EACH SECTION IS DIVIDED INTO NC1 C COLUMNS OF THE SAME WIDTH, PARALLEL WITH THE SECOND C NORMALIZED SECTION COORDINATE AXIS. C NC2... NUMBER OF STEPS IN THE DIRECTION OF THE SECOND NORMALIZED C SECTION COORDINATE AXIS (I.E. ALONG THE COLUMNS), WHEN C SEARCHING FOR THE INTERFACES OR ISOLINES. C STEP... STEP.GT.0: THE OUTPUT INTERFACES AND ISOLINES ARE WRITTEN C IN THE FORM OF SHORT LINES (INTERSECTIONS OF INTERFACES C AND ISOSURFACES WITH GIVEN SECTIONS, STORED BY PARTS). C THE FILE WITH LINES HAS THE NAME GIVEN BY INPUT C PARAMETER 'SECTS' ON THE INPUT LINE (1). C STEP... THE RELATIVE STEP OF THE NUMERICAL INTEGRATION C ALONG INTERFACES AND ISOLINES. MEASURED IN NORMALIZED C PLOT COORDINATES IN WHICH A SECTION IS REPRESENTED BY A C UNIT SQUARE (0,1)*(0,1). THE INITIAL POINTS OF THE C ISOLINES ARE DETERMINED ALONG THE AXES OF INDIVIDUAL C COLUMNS (SEE NC1 ABOVE). THEN THE ISOLINES ARE TRACED C BY MEANS OF NUMERICAL INTEGRATION. C STEP.EQ.0: NO OUTPUT FILE WITH THE SECTIONS IS GENERATED. C INSTEAD, THE SECTIONS ARE ROUGHLY DISPLAYED ON THE C SCREEN IN A TEXT MODE. THIS OPTION IS NOT ASSUMED TO BE C USED. C STEP.LT.0: THE OUTPUT INTERFACES AND ISOLINES ARE WRITTEN C IN THE FORM OF POINTS. THE POINTS ARE DETERMINED AS THE C POINT OF INTERSECTION OF INTERFACES AND ISOSURFACES WITH C THE GIVEN LINES. THE LINES ARE DEFINED AS THE AXES OF C INDIVIDUAL COLUMNS (SEE NC1 ABOVE). C THE FILE WITH POINTS HAS THE NAME GIVEN BY INPUT C PARAMETER 'SECTS' ON THE INPUT LINE (1). C THE ABSOLUTE VALUE OF STEP HAS NO SIGNIFICANCE FOR THIS C OPTION. C ERR... MAXIMUM RELATIVE ERROR IN THE DIRECTION OF COLUMNS WHEN C DETERMINING THE POSITIONS OF THE INTERFACES OR ISOLINES C AND, SIMULTANEOUSLY, THE UPPER ERROR BOUND OF THE C NUMERICAL INTEGRATION WHEN TRACING THE ISOLINES. MEASURED C IN NORMALIZED PLOT COORDINATES IN WHICH A SECTION IS C REPRESENTED BY A UNIT SQUARE (0,1)*(0,1). C DEFAULT VALUES: IPS=1, NC1=4, NC2=4, STEP=0.3/NC1, ERR=0.001. C (3) VALUES OF ISOLINES TO BE PLOTTED TERMINATED BY A SLASH. C (4) ANY TIMES THE FOLLOWING DATA (4.1): C (4.1) C10,C20,C30,C11,C21,C31,C12,C22,C32,C13,C23,C33,NREPET,/ C C10,C20,C30... CARTESIAN COORDINATES CORRESPONDING TO THE POINT OF C PLOT COORDINATES (0,0) (I.E. TO THE ORIGIN OF THE PLOT C CORDINATES). IF THE MODEL COORDINATES ARE CURVILINEAR, C THE MAPPING TO THE CARTESIAN COORDINATES IS GIVEN BY MEANS C OF THE SUBROUTINE CARTES. C C11,C21,C31... CARTESIAN COORDINATES CORRESPONDING TO THE POINT OF C PLOT COORDINATES (1,0) MINUS C10,C20,C30, RESPECTIVELY. C C12,C22,C32... CARTESIAN COORDINATES CORRESPONDING TO THE POINT OF C PLOT COORDINATES (0,1) MINUS C10,C20,C30, RESPECTIVELY. C C13,C23,C33... NEED NOT BE SPECIFIED FOR NREPET=0 (DEFAULT). C FOR NREPET POSITIVE, C13,C23,C33 IS THE VECTORIAL DISTANCE C BETWEEN THE ORIGINS OF THE TWO CONSECUTIVE PARALLEL C SECTIONS. C NREPET..IF POSITIVE, NREPET ADDITIONAL SECTIONS, PARALLEL WITH THE C GIVEN SECTION, WILL BE GENERATED. ORIGIN (C10,C20,C30) OF C EACH PARALLEL SECTION IS ORIGIN (C10,C20,C30) OF THE C PREVIOUS SECTION PLUS (C13,C23,C33). C /... OBLIGATORY SLASH AT THE END OF LINE TO FACILITATE FUTURE C EXTENSIONS. C DEFAULT: NREPET=0. C (5) / (A SLASH) C FOR AN EXAMPLE REFER TO THE SAMPLE INPUT DATA 'SEC.DAT'. C C INPUT FILE 'MODEL' IS DESCRIBED WITHIN THE SOURCE FILE 'MODEL.FOR'. C C OUTPUT FILE 'SECTS' WITH LINES (FOR STEP.GT.0): C (1) NONE TO SEVERAL STRINGS TERMINATED BY / (A SLASH). C (2) V1,V2,...,VN,/ C LIST OF ISOLINE VALUES TERMINATED BY A SLASH. C (3) FOR EACH MODEL SECTION (3.1), (3.2), AND (3.3): C (3.1) C10,C20,C30,C11,C21,C31,C12,C22,C32 C TRANSFORMATION MATRIX FROM NORMALIZED SECTION COORDINATES C1,C2 TO C THE CARTESIAN COORDINATES X1,X2,X3: C X1=C10+C11*C1+C12*C2, C X2=C20+C21*C1+C22*C2, C X3=C30+C31*C1+C32*C2. C (3.2) FOR EACH ISOLINE ELEMENT (3.2.1), (3.2.2), AND (3.2.3): C (3.2.1) LINE,IV C LINE... POSITIVE: INDEX OF THE COMPLEX BLOCK. THE ISOLINE IS A C VELOCITY ISOLINE. C NEGATIVE: MINUS THE INDEX OF A SURFACE. THE ISOLINE IS C A PART OF THE STRUCTURAL INTERFACE. C IV... INDEX OF THE VELOCITY VALUE CORRESPONDING TO THE ISOLINE C FOR NLINE POSITIVE, C 0 FOR NLINE NEGATIVE. C (3.2.2) FOR EACH POINT OF THE ELEMENT (3.2.2.1): C (3.2.2.1) C1,C2 C NORMALIZED SECTION COORDINATES OF A POINT OF THE ISOLINE. C (3.2.3) / (A SLASH). C (3.3) / (A SLASH). C (4) / (A SLASH) OR END OF FILE. C C OUTPUT FILE 'SECTS' WITH POINTS (FOR STEP.LT.0): C (1) FOR EACH POINT OF INTERSECTION OF THE VERTICAL LINE (REPRESENTING C THE AXIS OF THE VERTICAL COLUMN) THE FOLLOWING LINE: C (1.1) ISRF,X1,X2,X3 C ISRF... INDEX OF THE SURFACE COVERING THE STRUCTURAL INTERFACE AT C THE POINT OF INTERSECTION. C X1,X2,X3... CARTESIAN COORDINATES OF THE POINT OF INTERSECTION. C C STORAGE IN THE MEMORY: C THE INPUT DATA (2) TO (4) ARE STORED IN COMMON BLOCKS /SECTC/ AND C /VALUES/ DEFINED IN THE FOLLOWING SUBROUTINE: C ------------------------------------------------------------------ BLOCK DATA SECTB REAL C10,C20,C30,C11,C21,C31,C12,C22,C32 COMMON/SECTC/ C10,C20,C30,C11,C21,C31,C12,C22,C32 SAVE /SECTC/ INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV SAVE /VALUES/ REAL ZL(100),ZM(100),ZR(100),COLL,COLM,COLR INTEGER INTL,INTM,INTR COMMON/COLUMN/ZL,ZM,ZR,COLL,COLM,COLR,INTL,INTM,INTR SAVE /COLUMN/ INTEGER IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD COMMON/CONTC/ IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD END C ------------------------------------------------------------------ C C10,C20,C30,C11,C21,C31,C12,C22,C32... INPUT DATA (2). C IPS... INPUT DATA (3). C VALUE...INPUT DATA (4). C NV... NUMBER OF NON-ZERO ISOLINE VALUES IN THE INPUT DATA (4). C IV... INDEX IN THE ARRAY VALUE OF THE ISOLINE CURRENTLY BEING C PLOTTED. C ZERO, IF THERE IS NO ISOLINE IN THE GIVEN INTERVAL. C OUTPUT OF THE SUBROUTINE ISOL. IF THE INTERVAL WHERE AN C ISOLINE IS LOOKED FOR HAS BEEN CHANGED SINCE THE LAST C INVOCATION OF THE SUBROUTINE ISOL, IV HAS TO BE SET TO C ZERO. C ZL,ZM,ZR... THE SECOND NORMALIZED SECTION COORDINATES OF THE C POINTS OF INTERSECTION OF THE ISOLINES WITH LINES COLL, C COLM, AND COLR, RESPECTIVELY. C COLL... THE FIRST NORMALIZED SECTION COORDINATE OF THE LEFT-HAND C BOUNDARY LINE OF THE LEFT-HAND REGION FOR ISOLINE TRACING. C COLM... THE FIRST NORMALIZED SECTION COORDINATE OF THE MIDDLE C BOUNDARY LINE BETWEEN THE TWO REGIONS FOR ISOLINE TRACING. C THE INITIAL POINT OF THE ISOLINE IS ASSUMED TO BE SITUATED C AT THIS LINE. C COLR... THE FIRST NORMALIZED SECTION COORDINATE OF THE RIGHT-HAND C BOUNDARY LINE OF THE RIGHT-HAND REGION FOR ISOLINE C TRACING. C INTL,INTM,INTR... NUMBERS OF THE POINTS OF INTERSECTION OF THE C ISOLINES WITH LINES COLL, COLM, AND COLR, RESPECTIVELY. C IFUN... EITHER: C MINUS THE INDEX OF THE FUNCTION DESCRIBING A SURFACE C COVERING A STRUCTURAL INTERFACE, OR C -101, -102, -103, -104, -105 OR -106 FOR THE BOUNDARIES C OF THE MODEL, OR C THE INDEX OF THE COMPLEX BLOCK IN WHICH AN ISOLINE IS C PLOTTED. C NBACK...1 OR 2, INDEX OF THE DIRECTION IN WHICH THE ISOLINE IS TO C BE FOLLOWED. C ISB1O,ICB1O... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK IN WHICH THE INITIAL POINT OF THE ISOLINE IS C SITUATED. C ISB2O,ICB2O... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK, TOUCHING THE COMPLEX BLOCK ICB1O FROM THE OTHER C SIDE OF THE SURFACE ISRF=-IFUN AT THE INITIAL POINT. C NEED NOT BE DEFINED FOR A VELOCITY ISOLINE (IFUN C POSITIVE). C NBOD... NUMBER OF POINTS ALONG THE CALCULATED PART OF THE ISOLINE C (THE PART SITUATED WITHIN THE GIVEN COLUMN OF THE C SECTION). C COMMON BLOCK /SECTC/ IS INCLUDED IN EXTERNAL PROCEDURES FUNC, C SECT1, SECT2 AND SECT3. C COMMON BLOCK /VALUES/ IS INCLUDED IN EXTERNAL PROCEDURES FUNC, C ISOL, ISOLA, SECT1 AND MODSEC. C THE INDEX OF THE LAST ALLOCATED NUMERIC STORAGE UNIT OF ARRAY C VALUE IS NAMED MV AND IS GIVEN BY THE SIXTH AND SEVENTH STATEMENT C OF THE BLOCK DATA SUBROUTINE SECTB. IF THE VALUE OF MV IS C CHANGED, IT MUST BE ADJUSTED IN ALL SUBROUTINES WHICH INCLUDE THE C COMMON BLOCK /VALUES/. C C DATE: 1994, JANUARY 26 C CODED BY LUDEK KLIMES C C======================================================================= C SUBROUTINE FUNC(IFUNC,C,F) INTEGER IFUNC REAL C(2),F(3) C C THIS SUBROUTINE EVALUATES THE VALUE, FIRST AND SECOND DERIVATIVES OF C THE GIVEN FUNCTION WITH RESPECT TO THE PLOT COORDINATES. C C INPUT: C IFUNC...EITHER: C MINUS THE INDEX OF THE FUNCTION DESCRIBING A SURFACE C COVERING A STRUCTURAL INTERFACE, OR C -101, -102, -103, -104, -105 OR -106 FOR THE BOUNDARIES C OF THE MODEL, OR C THE INDEX OF THE COMPLEX BLOCK IN WHICH AN ISOLINE IS C PLOTTED. C C... ARRAY OF TWO PLOT COORDINATES OF THE GIVEN POINT. C C OUTPUT: C F... ARRAY CONTAINING FUNCTION VALUE AND THE FIRST PARTIAL C DERIVATIVES OF THE EVALUATED FUNCTION IN THE ORDER F, F1, C F2. C C COMMON BLOCK /MODELC/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /SECTC/: REAL C10,C20,C30,C11,C21,C31,C12,C22,C32 COMMON/SECTC/ C10,C20,C30,C11,C21,C31,C12,C22,C32 C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL VELOC,CARTES,SRFC2,PARM2 C VELOC...FILE 'MODEL.FOR'. C CARTES,KOOR... FILE 'METRIC.FOR'. C SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C C DATE: 1992, NOVEMBER 2 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I REAL AUX1,AUX2,AUX3 REAL COOR(3),CART(3),PDER(9) REAL UP(10),US(10),QP,QS,VD(10) C C....................................................................... C C CARTESIAN COORDINATES CART(1)=C10+C11*C(1)+C12*C(2) CART(2)=C20+C21*C(1)+C22*C(2) CART(3)=C30+C31*C(1)+C32*C(2) C C MODEL COORDINATES CALL CARTES(COOR,.FALSE.,CART,PDER) C C EVALUATING THE FUNCTION IN THE MODEL COORDINATES IF(IFUNC.GT.0) THEN CALL PARM2(IFUNC,COOR,UP,US,AUX3,QP,QS) CALL VELOC(IPS,UP,US,QP,QS,AUX1,AUX2,VD,AUX3) ELSE IF(IFUNC.GE.-100) THEN CALL SRFC2(-IFUNC,COOR,VD) ELSE I=(-IFUNC-99)/2 VD(1)=COOR(I)-BOUNDM(-IFUNC-100) VD(2)=0. VD(3)=0. VD(4)=0. VD(I+1)=1. END IF C C GRADIENT IN CARTESIAN COORDINATES AUX1=PDER(1)*VD(2)+PDER(2)*VD(3)+PDER(3)*VD(4) AUX2=PDER(4)*VD(2)+PDER(5)*VD(3)+PDER(6)*VD(4) AUX3=PDER(7)*VD(2)+PDER(8)*VD(3)+PDER(9)*VD(4) C C GRADIENT IN PLOT COORDINATES F(1)=VD(1) F(2)=C11*AUX1+C21*AUX2+C31*AUX3 F(3)=C12*AUX1+C22*AUX2+C32*AUX3 RETURN END C C======================================================================= C SUBROUTINE DISC(X1,C1,DC1,X2,C2,DC2,ERR,C1MIN,C1MAX,C2MIN,C2MAX, * LINE,ISB1,ICB1,ISRF,ISB2,ICB2) REAL X1,C1(2),DC1(2),X2,C2(2),DC2(2),ERR,C1MIN,C1MAX,C2MIN,C2MAX INTEGER LINE,ISB1,ICB1,ISRF,ISB2,ICB2 C C THIS SUBROUTINE DETERMINES THE POINT OF INTERSECTION OF THE GIVEN LINE C SEGMENT WITH THE BOUNDARY OF THE COMPLEX BLOCK. IT MAY ALSO BE USED C TO DETERMINE THE INDEX OF THE BLOCK IN WHICH THE GIVEN POINT IS C SITUATED. C C INPUT: C X1... INDEPENDENT VARIABLE CORRESPONDING TO THE FIRST GIVEN C POINT. C C1... ARRAY OF PLOT COORDINATES CORRESPONDING TO THE FIRST GIVEN C POINT. C DC1... ARRAY CONTAINING THE DERIVATIVES OF PLOT COORDINATES WITH C RESPECT TO THE INDEPENDENT VARIABLE AT THE FIRST GIVEN C POINT. C X2... INDEPENDENT VARIABLE CORRESPONDING TO THE SECOND GIVEN C POINT. C C2... ARRAY OF PLOT COORDINATES CORRESPONDING TO THE SECOND GIVEN C POINT. C DC2... ARRAY CONTAINING THE DERIVATIVES OF PLOT COORDINATES WITH C RESPECT TO THE INDEPENDENT VARIABLE AT THE SECOND GIVEN C POINT. C ERR... MAXIMUM ERROR IN INDEPENDENT VARIABLE FOR THE C DETERMINATION OF THE POINT OF INTERSECTION. C C1MIN,C1MAX,C2MIN,C2MAX... BOUNDARIES OF THE REGION IN WHICH THE C LINE SHOULD BE SITUATED. C LINE... IF THE LINE IS A SURFACE, INDEX OF THE SURFACE. C IN THIS CASE, INPUT VALUE OF ISB1 AND ICB1 MAY C CORRESPOND TO A MATERIAL BLOCK ON ITS EITHER SIDE. C IF THE LINE IS SITUATED INSIDE A BLOCK, LINE=0. C ISB1,ICB1... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK IN WHICH THE POINT C1 IS SITUATED, ICB1=0 IF THE C INDICES ARE YET UNKNOWN. C C OUTPUT: C X2,C2,DC2... INDEPENDENT VARIABLE, ARRAY OF PLOT COORDINATES AND C ARRAY CONTAINING THEIR DERIVATIVES, CORRESPONDING TO THE C ENDPOINT OF THE LINE ELEMENT. IF BOTH TWO GIVEN POINTS OF C THE LINE ARE SITUATED IN THE SAME COMPLEX BLOCK, THE C ENDPOINT OF THE LINE ELEMENT COINCIDES WITH THE SECOND C GIVEN POINT. OTHERWISE, THE ENDPOINT IS THE POINT OF C INTERSECTION OF THE LINE WITH THE BOUNDARY OF THE COMPLEX C BLOCK. C ISB1,ICB1... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK IN WHICH THE END OF LINE (AT C2) IS SITUATED. C ISB1=0 AND ICB1=0 IF THE POINT C1 IS SITUATED IN A FREE C SPACE. ISB1=0 AND ICB1=0 FOR THE POINT C1 BEING SITUATED C OUTSIDE THE COMPUTATIONAL VOLUME, TOO. C NOTE: IF THE POINT C1 IS SITUATED IN A FREE SPACE OR C OUTSIDE THE COMPUTATIONAL VOLUME (OUTPUT ICB1=0), C INDICES ISRF, ISB2 AND ICB2 ARE NOT DEFINED. C ISRF... INDEX OF THE SURFACE AT WHICH THE ENDPOINT OF THE LINE C ELEMENT IS SITUATED, SUPPLEMENTED BY A SIGN '+' OR '-' FOR C THE ENDPOINT SITUATED AT THE POSITIVE OR NEGATIVE SIDE OF C THE SURFACE, RESPECTIVELY. C ZERO INSIDE THE COMPLEX BLOCK. C ISRF=201 OR 202 IF CROSSING THE GIVEN BOUNDARIES C1MIN OR C C1MAX, RESPECTIVELY. C ISB2,ICB2... INDICES OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK AT POINT C2 FROM THE OTHER SIDE THAN THE LINE. C IF THE POINT C2 IS SITUATED AT THE BOUNDARY OF THE COMPLEX C BLOCK ICB1 (I.E. FOR ISRF.NE.0), ISB2,ICB2 ARE THE C INDICES OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK, TOUCHING THE COMPLEX BLOCK ICB1 FROM THE OTHER SIDE C OF THE SURFACE ISRF AT THE ENDPOINT OF THE LINE ELEMENT. C ISB2=0 AND ICB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF C ISRF. ISB2=0 AND ICB2=0 FOR THE SURFACE ISRF BEING THE C BOUNDARY OF THE COMPUTATIONAL VOLUME. C C COMMON BLOCK /MODELC/: INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: INTEGER NSRFC EXTERNAL BLOCK,SRFC2,CDE,NSRFC,SECT2,SECT3 C NSRFC,BLOCK... FILE 'MODEL.FOR'. C CARTES... FILE 'METRIC.FOR'. C SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C CDE,CROSS,HIVD2... FILE 'MEANS.FOR'. C SECT2,SECT3... THIS FILE. C C DATE: 1994, JANUARY 26 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I,IY(8) REAL FAUX(10),CBOUND(1) REAL Y1(3),D1(3),Y2(3),D2(3),XA,YA(3),DA(3),XB,YB(3),DB(3) C C....................................................................... C C DETERMINATION OF THE COMPLEX BLOCK IN WHICH THE FIRST POINT LIES CALL SECT2(C1,DC1,Y1,D1) IF(ICB1.EQ.0) THEN C CHECK FOR CROSSING THE END SURFACES * DO 11 I=1,NEND * CALL SRFC2(IABS(KEND(I)),Y1,FAUX) * IF(FAUX(1)*FLOAT(KEND(I)).LE.0.) THEN * GO TO 90 * END IF * 11 CONTINUE C CHECK FOR CROSSING THE COORDINATE BOUNDARIES OF THE C COMPUTATIONAL VOLUME: DO 12 I=1,3 IF(Y1(I).LT.BOUNDM(I+I-1)) THEN GO TO 90 END IF IF(Y1(I).GT.BOUNDM(I+I)) THEN GO TO 90 END IF 12 CONTINUE C SEARCHING FOR THE COMPLEX BLOCK IN WHICH THE FIRST POINT LIES: CALL BLOCK(Y1,0,0,ISRF,ISB1,ICB1,FAUX) IF(ICB1.EQ.0) THEN GO TO 90 END IF END IF C C INTERSECTION OF THE GIVEN LINE SEGMENT WITH THE BOUNDARIES C1MIN, C C1MAX OF THE REGION IN WHICH THE LINE SHOULD BE SITUATED ISRF=0 20 CONTINUE IF(ISRF.NE.201.AND.C1MIN.GT.C2(1)) THEN I=1 ISRF=201 FAUX(1)=C1MIN ELSE IF(ISRF.NE.202.AND.C2(1).GT.C1MAX) THEN I=1 ISRF=202 FAUX(1)=C1MAX ELSE IF(ISRF.NE.203.AND.C2MIN.GT.C2(2)) THEN I=2 ISRF=203 FAUX(1)=C2MIN ELSE IF(ISRF.NE.204.AND.C2(2).GT.C2MAX) THEN I=2 ISRF=204 FAUX(1)=C2MAX ELSE GO TO 30 END IF XA=X2 YA(1)=C2(1) YA(2)=C2(2) DA(1)=DC2(1) DA(2)=DC2(2) CALL CROSS(SECT2,101,I,I,2,ERR,X1,C1,DC1,X2,C2,DC2,XA,YA,DA, * XB,YB,DB,FAUX) C HERE SECT2 JUST FILLS THE PLACE RESERVED FOR EXTERNAL PROCEDURE. X2=XA C2(1)=YA(1) C2(2)=YA(2) DC2(1)=DA(1) DC2(2)=DA(2) GO TO 20 C C INTERSECTION OF THE GIVEN LINE SEGMENT WITH THE BOUNDARY OF THE C COMPLEX BLOCK 30 CONTINUE CALL SECT2(C2,DC2,Y2,D2) XA=X2 DO 31 I=1,3 YA(I)=Y2(I) DA(I)=D2(I) 31 CONTINUE IY(4)=ISB1 IY(5)=ICB1 IY(6)=ISRF CALL CDE(IABS(LINE),0,I,0,I,CBOUND, * 1,3,3,IY,ERR,X1,X1,Y1,D1,X2,Y2,D2,XA,YA,DA,XB,YB,DB) IF(IY(6).EQ.ISRF) THEN ISB2=IY(4) ICB2=IY(5) ELSE IF(IABS(IY(6)).LE.NSRFC()) THEN ISB2=IY(7) ICB2=IY(8) ELSE ISB2=0 ICB2=0 END IF ISB1=IY(4) ISRF=IY(6) X2=XA CALL SECT3(YA,DA,C2,DC2) RETURN C 90 CONTINUE ISB1=0 ISRF=0 ISB2=0 ICB2=0 RETURN END C C======================================================================= C SUBROUTINE ISOL(X1,C1,DC1,X2,C2,DC2,ERR,ICB1) REAL X1,C1(2),DC1(2),X2,C2(2),DC2(2),ERR INTEGER ICB1 C C THIS SUBROUTINE FINDS THE POINT OF INTERSECTION OF THE GIVEN LINE C SEGMENT WITH AN ISOLINE. C C INPUT: C X1... INDEPENDENT VARIABLE CORRESPONDING TO THE FIRST GIVEN C POINT. C C1... ARRAY OF PLOT COORDINATES CORRESPONDING TO THE FIRST GIVEN C POINT. C DC1... ARRAY CONTAINING THE DERIVATIVES OF PLOT COORDINATES WITH C RESPECT TO THE INDEPENDENT VARIABLE AT THE FIRST GIVEN C POINT. C X2... INDEPENDENT VARIABLE CORRESPONDING TO THE SECOND GIVEN C POINT. C C2... ARRAY OF PLOT COORDINATES CORRESPONDING TO THE SECOND C GIVEN POINT. C DC2... ARRAY CONTAINING THE DERIVATIVES OF PLOT COORDINATES WITH C RESPECT TO THE INDEPENDENT VARIABLE AT THE SECOND GIVEN C POINT. C ERR... MAXIMUM ERROR IN INDEPENDENT VARIABLE FOR THE C DETERMINATION OF THE POINT OF INTERSECTION. C ICB1... INDEX OF THE COMPLEX BLOCK IN WHICH THE POINTS C1 AND C2 C ARE SITUATED. C C OUTPUT: C X1,C1,DC1... INDEPENDENT VARIABLE, ARRAY OF PLOT COORDINATES AND C ARRAY CONTAINING THEIR DERIVATIVES, CORRESPONDING TO THE C POINT OF INTERSECTION OF THE LINE WITH THE ISOLINE CLOSEST C TO THE FIRST GIVEN POINT. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C INPUT: C IPS,VALUE,NV... GIVEN VALUES, SEE THE DESCRIPTION OF THE BLOCK C DATA SECTB. C IV... IF THE INPUT VALUES OF X1,C1,DC1 ARE THE OUTPUT ONES FROM C THE LAST INVOCATION OF THIS SUBROUTINE, IV SHOULD BE C UNCHANGED SINCE THAT INVOCATION. C IF THE INTERVAL X1,C1,DC1 TO X2,C2,DC2 HAS BEEN CHANGED C SINCE THE LAST INVOCATION OF THIS SUBROUTINE, IV HAS TO C BE SET TO ZERO. C OUTPUT: C IV... INDEX IN THE ARRAY VALUE OF THE ISOLINE CLOSEST TO THE C FIRST GIVEN POINT. C ZERO, IF THERE IS NO ISOLINE IN THE GIVEN INTERVAL. C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK, EXCEPT IV, ARE C ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL CROSS,FUNC,ISOLA,SECT2,SECT3 C VELOC... FILE 'MODEL.FOR'. C CARTES,KOOR... FILE 'METRIC.FOR'. C SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C CROSS,HIVD2... FILE 'MEANS.FOR'. C FUNC,ISOLA,SECT2,SECT3... THIS FILE. C C DATE: 1992, NOVEMBER 2 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I REAL F1(3),F2(10) REAL Y1(3),D1(3),Y2(3),D2(3),XA,YA(3),DA(3),XB,YB(3),DB(3) C C....................................................................... C CALL FUNC(ICB1,C1,F1) CALL FUNC(ICB1,C2,F2) IF(F1(1).LT.F2(1)) THEN IF(IV.EQ.0) THEN DO 11 IV=1,NV IF(VALUE(IV).GT.F1(1)) THEN GO TO 12 END IF 11 CONTINUE IV=NV+1 12 CONTINUE ELSE IV=IV+1 END IF IF(IV.GT.NV) THEN IV=0 ELSE IF(VALUE(IV).GT.F2(1)) THEN IV=0 END IF ELSE IF(IV.EQ.0) THEN DO 13 IV=NV,1,-1 IF(VALUE(IV).LT.F1(1)) THEN GO TO 14 END IF 13 CONTINUE IV=0 14 CONTINUE ELSE IV=IV-1 END IF IF(IV.GT.0) THEN IF(VALUE(IV).LT.F2(1)) THEN IV=0 END IF END IF END IF C IF(IV.GT.0) THEN C THERE IS AN ISOLINE NO. IV IN THE GIVEN INTERVAL CALL SECT2(C1,DC1,Y1,D1) CALL SECT2(C2,DC2,Y2,D2) XA=X2 DO 21 I=1,3 YA(I)=Y2(I) DA(I)=D2(I) 21 CONTINUE CALL ISOLA(ICB1,YA,F2) CALL CROSS(ISOLA,ICB1,1,3,3,ERR,X1,Y1,D1,X2,Y2,D2,XA,YA,DA, * XB,YB,DB,F2) X1=XB CALL SECT3(YA,DA,C1,DC1) END IF C RETURN END C C======================================================================= C SUBROUTINE ISOLA(ICB1,COOR,F) INTEGER ICB1 REAL COOR(3),F(10) C C AUXILIARY SUBROUTINE TO THE SUBROUTINE ISOL. C THIS SUBROUTINE EVALUATES THE VALUE, FIRST AND SECOND DERIVATIVES OF C THE GIVEN FUNCTION WITH RESPECT TO THE MODEL COORDINATES. IT IS C INTENDED TO BE CALLED BY THE SUBROUTINE CROSS OF THE FILE 'MEANS.FOR' C IN ORDER TO FIND THE POINT OF INTERSECTION OF THE GIVEN LINE SEGMENT C WITH AN ISOLINE. NOTE THAT THE ISOLINE IS THE ZERO ISOLINE OF THE C FUNCTION EVALUATED BY THIS SUBROUTINE. C C INPUT: C ICB1... INDEX OF THE COMPLEX BLOCK IN WHICH THE POINTS C1 AND C2 C ARE SITUATED. C COOR... ARRAY OF THREE MODEL COORDINATES OF THE GIVEN POINT. C C OUTPUT: C F... ARRAY CONTAINING FUNCTION VALUE, THE FIRST AND SECOND C PARTIAL DERIVATIVES OF THE EVALUATED FUNCTION IN THE C ORDER F, F1, F2, F3, F11, F12, F22, F13, F23, F33. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: C EXTERNAL VELOC,PARM2 C VELOC... FILE 'MODEL.FOR'. C PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR'). C C DATE: 1991, JUNE 25 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: REAL AUX1,AUX2,AUX3 REAL UP(10),US(10),QP,QS C C....................................................................... C C EVALUATING THE FUNCTION IN THE MODEL COORDINATES CALL PARM2(ICB1,COOR,UP,US,AUX3,QP,QS) CALL VELOC(IPS,UP,US,QP,QS,AUX1,AUX2,F,AUX3) F(1)=F(1)-VALUE(IV) RETURN END C C======================================================================= C SUBROUTINE SECT1(LU1,LU2,ISECT,NC1,NC2,STEP,ERR) INTEGER LU1,LU2,ISECT,NC1,NC2 REAL STEP,ERR C C THIS SUBROUTINE READS THE DATA FOR THE PLOTTED SECTION OF THE MODEL. C C INPUT: C LU1... LOGICAL NUMBER OF THE EXTERNAL UNIT CONNECTED TO THE FILE C WITH THE MAIN INPUT DATA. C IF ZERO, MAIN INPUT DATA ARE READ FROM THE * EXTERNAL C INTERACTIVE DEVICE. C LU2... LOGICAL UNIT NUMBER OF THE EXTERNAL OUTPUT DEVICE TO STORE C THE LINES OR POINTS SITUATED AT STRUCTURAL INTERFACES OR C VELOCITY ISOSURFACES. C THIS LOGICAL UNIT IS USED FOR READING THE INPUT DATA, C THEN CONNECTED TO THE OUTPUT FILE OPENED IN THIS ROUTINE. C NOTE THAT THE OUTPUT FILE IS NOT OPENED IF STEP=0 IN THE C INPUT DATA FILE. C ISECT...FOR THE FIRST INVOCATION ISECT=0, C OTHERWISE THE UNCHANGED OUTPUT FROM THE PREVIOUS C INVOCATION. C C OUTPUT: C ISECT...ISECT=0 IF THERE REMAINS NO MODEL SECTION TO COMPUTE. C THEN THE OUTPUT PARAMETERS NC1,NC2,STEP,ERR ARE C UNCHANGED INPUT. C OTHERWISE THE INPUT VALUE INCREASED BY ONE. C THEN THE OUTPUT PARAMETERS NC1,NC2,STEP,ERR ARE READ C FROM THE INPUT DATA FILE. C NC1... NUMBER OF VERTICAL COLUMNS IN THE PRINT OF THE MODEL C SECTION. C NC2... NUMBER OF STEPS IN THE VERTICAL DIRECTION WHEN LOOKING FOR C THE INTERFACES OR ISOLINES. C STEP... RELATIVE STEP OF THE NUMERICAL INTEGRATION ALONG C INTERFACES AND ISOLINES. C ERR... RELATIVE ERROR IN THE VERTICAL DIRECTION WHEN DETERMINING C THE POSITIONS OF THE INTERFACES OR ISOLINES AND, C SIMULTANEOUSLY, THE UPPER ERROR BOUND OF THE NUMERICAL C INTEGRATION. C C COMMON BLOCK /MODELT/: CHARACTER*80 TEXTM COMMON/MODELT/TEXTM C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /SECTC/: REAL C10,C20,C30,C11,C21,C31,C12,C22,C32 COMMON/SECTC/ C10,C20,C30,C11,C21,C31,C12,C22,C32 C THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE DEFINED IN THIS C SUBROUTINE. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE DEFINED IN THIS C SUBROUTINE. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL SECTB,MODEL1 C SECTB.. BLOCK DATA SUBROUTINE OF THIS FILE. C MODEL1 AND SUBSEQUENT ROUTINES... FILE 'MODEL.FOR' AND SUBSEQUENT C FILES (LIKE 'METRIC.FOR', 'SRFC.FOR', 'PARM.FOR', C 'VAL.FOR', AND 'FIT.FOR'). C C DATE: 1994, JANUARY 22 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUTOMATIC GENERATION OF PARALLEL SECTIONS: INTEGER NREPET REAL C13,C23,C33 SAVE C13,C23,C33,NREPET C NREPET..NUMBER OF SECTIONS PARALLEL WITH THE GIVEN SECTION. C C13,C23,C33... OFFSET OF THE ORIGINS OF THE CONSECUTIVE PARALLEL C SECTIONS. C C AUXILIARY STORAGE LOCATIONS: CHARACTER*80 MODEL,SECTS REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8,AUX9 C C....................................................................... C IF(ISECT.EQ.0) THEN C C (1) OPENING DATA FILES AND READING THE INPUT DATA: MODEL='MODEL.DAT' SECTS='SEC.OUT' IF(LU1.EQ.0) THEN READ(*,*) MODEL,SECTS ELSE READ(LU1,*) MODEL,SECTS END IF OPEN(LU2,FILE=MODEL,STATUS='OLD') CALL MODEL1(LU2) CLOSE(LU2) C C (2) IPS=1 NC1=4 NC2=4 STEP=0.3/FLOAT(NC1) ERR=0.001 IF(LU1.EQ.0) THEN READ(*,*) IPS,NC1,NC2,STEP,ERR ELSE READ(LU1,*) IPS,NC1,NC2,STEP,ERR END IF C C (3) DO 11 NV=1,MV VALUE(NV)=0. 11 CONTINUE IF(LU1.EQ.0) THEN READ(*,*) VALUE ELSE READ(LU1,*) VALUE END IF DO 12 NV=1,MV IF(VALUE(NV).EQ.0.) THEN GO TO 13 END IF 12 CONTINUE NV=MV+1 13 CONTINUE NV=NV-1 IV=0 C C OUTPUT FILE: IF(STEP.NE.0.) THEN OPEN(LU2,FILE=SECTS) IF(STEP.GT.0.) THEN WRITE(LU2,'(3A)') * '''LINES SITUATED AT INTERFACES OR VELOCITY ISOSURFACES:''' ELSE WRITE(LU2,'(3A)') * '''POINTS SITUATED AT INTERFACES OR VELOCITY ISOSURFACES:''' END IF WRITE(LU2,'(3A)') '''',TEXTM(1:78),'''' WRITE(LU2,'(A)') '/' IF(STEP.GT.0.) THEN C TERMINAL LINE OF FILE, OVERWRITTEN USING BACKSPACE IN THE C CASE OF OUTPUT WRITE(LU2,'(A)') '/' END IF END IF C C INITIALIZATION: NREPET=0 C END IF C ISECT=ISECT+1 IF(NREPET.LE.0) THEN C C (4) AUX1=C10 AUX2=C20 AUX3=C30 AUX4=C11 AUX5=C21 AUX6=C31 AUX7=C12 AUX8=C22 AUX9=C32 IF(LU1.EQ.0) THEN READ(*,*) C10,C20,C30,C11,C21,C31,C12,C22,C32,C13,C23,C33, * NREPET ELSE READ(LU1,*) C10,C20,C30,C11,C21,C31,C12,C22,C32,C13,C23,C33, * NREPET END IF IF(AUX1.EQ.C10.AND.AUX4.EQ.C11.AND.AUX7.EQ.C12.AND. * AUX2.EQ.C20.AND.AUX5.EQ.C21.AND.AUX8.EQ.C22.AND. * AUX3.EQ.C30.AND.AUX6.EQ.C31.AND.AUX9.EQ.C32) THEN C END OF COMPUTATION ISECT=0 END IF C ELSE C C PARALLEL SECTION WITH THE PREVIOUS ONE: C10=C10+C13 C20=C20+C23 C30=C30+C33 NREPET=NREPET-1 C END IF RETURN END C C======================================================================= C SUBROUTINE SECT2(C,DC,COOR,DCOOR) REAL C(2),DC(2),COOR(3),DCOOR(3) C C AUXILIARY SUBROUTINE TO THE SUBROUTINES DISC AND ISOL. C THIS SUBROUTINE TRANSFORMS THE PLOT COORDINATES TO THE MODEL C COORDINATES. C C INPUT: C C... ARRAY OF TWO PLOT COORDINATES OF THE GIVEN POINT. C C OUTPUT: C COOR... ARRAY CONTAINING THE MODEL COORDINATES X1, X2, X3 OF THE C GIVEN POINT. C C COMMON BLOCK /SECTC/: REAL C10,C20,C30,C11,C21,C31,C12,C22,C32 COMMON/SECTC/ C10,C20,C30,C11,C21,C31,C12,C22,C32 C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL CARTES C CARTES,KOOR... FILE 'METRIC.FOR'. C C DATE: 1991, OCTOBER 7 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: REAL CART(3),PDER(9),AUX1,AUX2,AUX3 C C....................................................................... C C CARTESIAN COORDINATES: CART(1)=C10+C11*C(1)+C12*C(2) CART(2)=C20+C21*C(1)+C22*C(2) CART(3)=C30+C31*C(1)+C32*C(2) AUX1=C11*DC(1)+C12*DC(2) AUX2=C21*DC(1)+C22*DC(2) AUX3=C31*DC(1)+C32*DC(2) C C MODEL COORDINATES: CALL CARTES(COOR,.FALSE.,CART,PDER) DCOOR(1)=PDER(1)*AUX1+PDER(4)*AUX2+PDER(7)*AUX3 DCOOR(2)=PDER(2)*AUX1+PDER(5)*AUX2+PDER(8)*AUX3 DCOOR(3)=PDER(3)*AUX1+PDER(6)*AUX2+PDER(9)*AUX3 C RETURN END C C======================================================================= C SUBROUTINE SECT3(COOR,DCOOR,C,DC) REAL COOR(3),DCOOR(3),C(2),DC(2) C C AUXILIARY SUBROUTINE TO THE SUBROUTINES DISC AND ISOL. C THIS SUBROUTINE TRANSFORMS THE MODEL COORDINATES TO THE PLOT C COORDINATES. C C INPUT: C COOR... ARRAY CONTAINING THE MODEL COORDINATES X1, X2, X3 OF THE C GIVEN POINT. C C OUTPUT: C C... ARRAY OF TWO PLOT COORDINATES OF THE GIVEN POINT. C C COMMON BLOCK /SECTC/: REAL C10,C20,C30,C11,C21,C31,C12,C22,C32 COMMON/SECTC/ C10,C20,C30,C11,C21,C31,C12,C22,C32 C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL CARTES C CARTES,KOOR... FILE 'METRIC.FOR'. C C DATE: 1991, OCTOBER 7 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: REAL CART(3),AUX1,AUX2 REAL B11,B12,B22,BDET,DAUX1,DAUX2,DAUX3,PDER(9) C C....................................................................... C C CARTESIAN COORDINATES: CALL CARTES(COOR,.TRUE.,CART,PDER) DAUX1=PDER(1)*DCOOR(1)+PDER(4)*DCOOR(2)+PDER(7)*DCOOR(3) DAUX2=PDER(2)*DCOOR(1)+PDER(5)*DCOOR(2)+PDER(8)*DCOOR(3) DAUX3=PDER(3)*DCOOR(1)+PDER(6)*DCOOR(2)+PDER(9)*DCOOR(3) C C PLOT COORDINATES: B11=C11*C11+C21*C21+C31*C31 B12=C11*C12+C21*C22+C31*C32 B22=C12*C12+C22*C22+C32*C32 BDET=B11*B22-B12*B12 AUX1=(CART(1)-C10)*C11+(CART(2)-C20)*C21+(CART(3)-C30)*C31 AUX2=(CART(1)-C10)*C12+(CART(2)-C20)*C22+(CART(3)-C30)*C32 C(1)=( B22*AUX1-B12*AUX2)/BDET C(2)=(-B12*AUX1+B11*AUX2)/BDET AUX1=DAUX1*C11+DAUX2*C21+DAUX3*C31 AUX2=DAUX1*C12+DAUX2*C22+DAUX3*C32 DC(1)=( B22*AUX1-B12*AUX2)/BDET DC(2)=(-B12*AUX1+B11*AUX2)/BDET C RETURN END C C======================================================================= C SUBROUTINE MODSEC C C SUBROUTINE DEMONSTRATING THE FUNCTION OF THE SUBROUTINES DISC AND ISOL C (AND, PARTIALLY, ALSO FUNC). IT EMPLOYES THE SUBROUTINES WHEN C SKETCHING A SECTION OF THE MODEL BY MEANS OF EXTENDED ASCII C CHARACTERS ONTO THE SCREEN. THIS SUBROUTINE IS INTENDED TO BE JUST AN C EXAMPLE HOW TO USE THE SUBROUTINES FUNC, DISC AND ISOL OF THIS FILE. C C NO ARGUMENT INPUT NOR OUTPUT. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK, EXCEPT IV, ARE C ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: INTEGER NSRFC EXTERNAL NSRFC,DISC,ISOL,SECT1,CONT1,CONT2 C C DATE: 1994, JANUARY 26 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C CHARACTER*80 FILE INTEGER LU1,LU2,ISECT,NC1,NC2 REAL STEP,ERR C C FILE... NAME OF THE MAIN INPUT DATA FILE. C LU1... LOGICAL UNIT NUMBER OF THE INPUT DATA FILE. C LU2... LOGICAL UNIT NUMBER OF THE OUTPUT FILE WITH GENERATED C MODEL SECTIONS. C ISECT...SEQUENTIAL NUMBER OF THE CURRENTLY EVALUATED MODEL C SECTION. C NC1... NUMBER OF VERTICAL COLUMNS IN THE MODEL SECTION. C NC2... 1/NC2 IS THE RELATIVE STEP IN THE VERTICAL DIRECTION WHEN C LOOKING FOR THE INTERFACES OR ISOLINES. C STEP... STEP OF THE NUMERICAL INTEGRATION WHEN COMPUTING THE C INTERFACES OR ISOLINES. C ERR... RELATIVE ERROR IN THE VERTICAL DIRECTION WHEN DETERMINING C THE POSITIONS OF THE INTERFACES OR ISOLINES. C INTEGER MSRFC PARAMETER (MSRFC=100) INTEGER KVALUE(MV),KSRFC(MSRFC) INTEGER ISB1,ICB1,ISRF,ISB2,ICB2,ISB0,ICB0,I REAL X1,C1(2),DC1(2),X2,C2(2),DC2(2),C0,DIRECT REAL COOR(3),DCOOR(3) C C KVALUE(I)... THE LAST COLUMN INTERSECTED BY THE I-TH ISOLINE. C KSRFC(I)... THE LAST COLUMN INTERSECTED BY THE I-TH SURFACE. C ISB1,ICB1... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK IN WHICH THE POINT C1 IS SITUATED. C ISRF... INDEX OF THE SURFACE AT WHICH THE ENDPOINT OF THE LINE C ELEMENT IS SITUATED, SUPPLEMENTED BY A SIGN '+' OR '-' FOR C THE ENDPOINT SITUATED AT THE POSITIVE OR NEGATIVE SIDE OF C THE SURFACE, RESPECTIVELY. C ZERO INSIDE THE COMPLEX BLOCK. C ISB2,ICB2... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK, TOUCHING THE COMPLEX BLOCK ICB1 FROM THE OTHER SIDE C OF THE SURFACE ISRF AT THE ENDPOINT OF THE LINE ELEMENT. C ISB2=0 AND ICB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF C ISRF. C UNDEFINED INSIDE THE COMPLEX BLOCK, DEFINED ONLY AT THE C POINT OF INTERSECTION OF THE LINE WITH THE BOUNDARY OF THE C COMPLEX BLOCK. C ISB0,ICB0... LAST VALUES OF ISB1,ICB1 CORRESPONDING TO DIRECT=+1, C WHEN DIRECT=-1. C I... AUXILIARY LOOP VARRIABLE. C X1,X2...INDEPENDENT VARIABLES ALONG THE LINES, IN THIS CASE C COINCIDING WITH C1(2) AND C2(2), RESPECTIVELY. C C1,C2...VECTORS OF THE RELATIVE POSITION WITHIN THE MODEL SECTION. C DC1,DC2... DERIVATIVES OF THE ARRAYS C1,C2 WITH RESPECT TO THE C INDEPENDENT VARIABLE. C C0... C DIRECT..DIRECTION IN WHICH WE ARE SEARCHING FOR THE INTERFACES OR C ISOLINES, TAKES THE VALUES +1 OR -1. C INTEGER IFUNC,IBACK INTEGER MC1,MC2,IC1,IC2 PARAMETER(MC1=75,MC2=20) CHARACTER*(MC1) PICT(MC2) C C IFUNC...EITHER: C MINUS THE INDEX OF THE FUNCTION DESCRIBING A SURFACE C COVERING A STRUCTURAL INTERFACE, OR C -101, -102, -103, -104, -105 OR -106 FOR THE BOUNDARIES C OF THE MODEL, OR C THE INDEX OF THE COMPLEX BLOCK IN WHICH AN ISOLINE IS C PLOTTED. C IBACK...LOOP VARIABLE. C MC1... MAXIMUM NUMBER OF VERTICAL COLUMNS IN THE PRINT OF THE C MODEL SECTION. C MC2... NUMBER OF HORIZONTAL ROWS IN THE PRINT OF THE MODEL C SECTION. C IC1,IC2...INDICES. C PICT...CHARACTER ARRAY CONTAINING THE PRINT OF THE MODEL. C C....................................................................... C LU1=11 LU2=12 C C OPENING MAIN INPUT DATA FILE: IF(LU1.NE.0) THEN WRITE(*,'(A)') ' ENTER THE NAME OF THE MAIN INPUT DATA FILE:' FILE='SEC.DAT' READ(*,*) FILE OPEN(LU1,FILE=FILE,STATUS='OLD') WRITE(*,'(A)') '+ ' END IF C C LOOP FOR SECTIONS OF THE MODEL ISECT=0 10 CONTINUE C C READING THE INPUT DATA CALL SECT1(LU1,LU2,ISECT,NC1,NC2,STEP,ERR) IF(ISECT.EQ.0) THEN C .............................................................. C END OF CALCULATION IF(STEP.NE.0.) THEN IF(STEP.NE.0.) THEN IF(STEP.GT.0.) THEN BACKSPACE(LU2) END IF WRITE(LU2,'(A)') '/' END IF CLOSE(LU2) END IF C .............................................................. RETURN END IF C DC1(1)=0. DC2(1)=0. C ................................................................ C STARTING WITH A NEW SECTION: IF(STEP.EQ.0.) THEN IF(NC1.GT.MC1) THEN PAUSE 'ERROR: CHARACTER ARRAY PICT TOO SMALL' STOP END IF DO 11 IC2=1,MC2 PICT(IC2)=' ' 11 CONTINUE IF(ISECT.GT.1) THEN WRITE(*,'(A)') ' ' END IF ELSE DO 16 I=1,MIN0(NSRFC(),MSRFC) KSRFC(I)=-1 16 CONTINUE DO 17 I=1,NV KVALUE(I)=-1 17 CONTINUE END IF C ................................................................ C C LOOP FOR THE VERTICAL LINES (COLUMNS IN THE PRINT OF THE MODEL) DO 80 IC1=1,NC1 DIRECT=1. C1(1)=(FLOAT(IC1)-.5)/FLOAT(NC1) C2(1)=C1(1) C1(2)=0. ISB1=0 ICB1=0 C .............................................................. C NEW COLUMN: IF(STEP.GT.0.) THEN CALL CONT1(IC1,NC1) END IF C .............................................................. C LOOP FOR THE INTERVALS OF THE SIZE 1/NC2 ALONG A VERTICAL LINE C WHEN LOOKING FOR THE INTERFACES OR ISOLINES 20 CONTINUE C2(2)=AMIN1(C1(2)+DIRECT/FLOAT(NC2),1.) DC1(2)=DIRECT DC2(2)=DIRECT X1=DIRECT*C1(2) X2=DIRECT*C2(2) IC2=INT(FLOAT(NC2)*(C1(2)+C2(2))/2.+1.) WRITE(*,'(9(A,I4,5X))') * '+SECTION',ISECT,'COLUMN',IC1,'ELEVATION',IC2 CALL DISC(X1,C1,DC1,X2,C2,DC2,ERR,0.,1.,0.,1., * 0,ISB1,ICB1,ISRF,ISB2,ICB2) IF(ICB1.EQ.0) THEN C POINT C1 IS SITUATED WITHIN A FREE SPACE ICB0=0 C1(2)=C1(2)+1./FLOAT(NC2) DIRECT=-1. ELSE C POINT C1 IS SITUATED WITHIN A MATERIAL COMPLEX BLOCK IF(ICB0.EQ.0) THEN ISB0=ISB1 ICB0=ICB1 C0=C1(2) END IF IV=0 30 CONTINUE CALL ISOL(X1,C1,DC1,X2,C2,DC2,ERR,ICB1) IF(IV.NE.0) THEN C POINT OF INTERSECTION OF THE VERTICAL LINE WITH AN C ISOLINE C ........................................................ IF(STEP.EQ.0.) THEN IC2= * MAX0(0,MIN0(MC2,INT((1.-C1(2)+ERR/2.)*FLOAT(MC2)+1.))) IF(PICT(IC2)(IC1:IC1).EQ.' ') THEN C PLOTTING ISOLINE IF(C1(2).GT.1.-(FLOAT(IC2)-.67)/FLOAT(MC2)) THEN PICT(IC2)(IC1:IC1)='~' ELSE IF(C1(2).GT.1.-(FLOAT(IC2)-.33)/FLOAT(MC2))THEN PICT(IC2)(IC1:IC1)='-' ELSE PICT(IC2)(IC1:IC1)='_' END IF END IF ELSE CALL SECT2(C1,DC1,COOR,DCOOR) IF(STEP.LT.0.) THEN IBACK=1 ELSE IBACK=2 END IF DO 31 IBACK=1,IBACK IF(STEP.GT.0.) THEN BACKSPACE(LU2) END IF IF(KVALUE(IV).LT.IC1-1.OR.STEP.LT.0.) THEN C WRITING THE REFERENCE COORDINATES: IF(IPS.GT.0) THEN WRITE(LU2,'(3(A,I4),A,F6.3,A,3F12.6,A)') * '''SECT',ISECT,', BLOC',ICB1,', ISOL',IV, * ', VP =', * VALUE(IV),'''',COOR(1),COOR(2),COOR(3),' /' ELSE WRITE(LU2,'(3(A,I4),A,F6.3,A,3F12.6,A)') * '''SECT',ISECT,', BLOC',ICB1,', ISOL',IV, * ', VS =', * VALUE(IV),'''',COOR(1),COOR(2),COOR(3),' /' END IF ELSE C WRITING / IN PLACE OF REFERENCE COORDINATES: IF(IPS.GT.0) THEN WRITE(LU2,'(3(A,I4),A,F6.3,A)') * '''SECT',ISECT,', BLOC',ICB1,', ISOL',IV, * ', VP =',VALUE(IV),''' /' ELSE WRITE(LU2,'(3(A,I4),A,F6.3,A)') * '''SECT',ISECT,', BLOC',ICB1,', ISOL',IV, * ', VS =',VALUE(IV),''' /' END IF END IF KVALUE(IV)=IC1 IF (STEP.GT.0.) THEN IFUNC=ICB1 CALL CONT2(LU2,IFUNC,IBACK,ISB1,ICB1,0,0,STEP,ERR, * C1) END IF 31 CONTINUE END IF C ........................................................ GO TO 30 END IF IF(ISRF.NE.0) THEN C CROSSING AN INTERFACE C ........................................................ IF(STEP.EQ.0.) THEN IC2=MAX0(0,MIN0(MC2,INT((1.-C2(2))*FLOAT(MC2)+1.))) IF(C2(2).GT.1.-(FLOAT(IC2)-.5)/FLOAT(MC2)) THEN PICT(IC2)(IC1:IC1)=CHAR(223) ELSE PICT(IC2)(IC1:IC1)=CHAR(220) END IF ELSE CALL SECT2(C2,DC2,COOR,DCOOR) IF(STEP.LT.0.) THEN IBACK=1 ELSE IBACK=2 END IF DO 36 IBACK=1,IBACK IF(STEP.GT.0.) THEN BACKSPACE(LU2) END IF IFUNC=IABS(ISRF) IF(KSRFC(MIN0(IFUNC,MSRFC)).LT.IC1-1.OR.STEP.LT.0.) * THEN C WRITING THE REFERENCE COORDINATES: WRITE(LU2,'(2(A,I4),A,3F12.6,A)') * '''SECT',ISECT,', SURF',IFUNC,'''', * COOR(1),COOR(2),COOR(3),' /' ELSE C WRITING / IN PLACE OF REFERENCE COORDINATES: WRITE(LU2,'(2(A,I4),A)') * '''SECT',ISECT,', SURF',IFUNC,''' /' END IF IF(IFUNC.LT.MSRFC) THEN KSRFC(IFUNC)=IC1 END IF IF(STEP.GT.0.) THEN IFUNC=-IFUNC C IFUNC=101,...,106 FOR THE COMPUTATIONAL VOLUME C BOUDARY. CALL CONT2(LU2,IFUNC,IBACK,ISB1,ICB1,ISB2,ICB2, * STEP,ERR,C2) END IF 36 CONTINUE END IF C ........................................................ IF(ICB2.EQ.0) THEN C CROSSING FREE SURFACE IF(DIRECT.LT.0.) THEN ISB2=ISB0 ICB2=ICB0 C2(2)=C0 ELSE C2(2)=C1(2)+1./FLOAT(NC2) END IF DIRECT=-DIRECT END IF END IF ISB1=ISB2 ICB1=ICB2 C1(2)=C2(2) END IF IF(C1(2).LT.0.) THEN PAUSE * 'ERROR: NEGATIVE PLOT COORDINATE - CONTACT THE AUTHOR' STOP END IF IF(C1(2).LT.1.-ERR) GO TO 20 80 CONTINUE C C ................................................................ C SECTION FINISHED: IF(STEP.EQ.0.) THEN WRITE(*,'(1X,999A1)') CHAR(218),(CHAR(196),I=1,NC1),CHAR(191) DO 90 IC2=1,MC2 WRITE(*,'(1X,A1,A,A1)') CHAR(179),PICT(IC2)(1:NC1),CHAR(179) 90 CONTINUE WRITE(*,'(1X,999A1)') CHAR(192),(CHAR(196),I=1,NC1),CHAR(217) END IF C ................................................................ C GO TO 10 C END OF LOOP FOR SECTIONS OF THE MODEL C END C C======================================================================= C SUBROUTINE CONT1(IC1,NC1) INTEGER IC1,NC1 C C SUBROUTINE DESIGNED TO INITIALIZE ARRAYS CONTAINING THE POINTS OF C INTERSECTION OF ISOLINES WITH VERTICAL LINES LIMITING THE REGIONS OF C NUMERICALY TRACING THE ISOLINES. IT IS CALLED ONCE BEFORE TRACING C ISOLINES IN A NEW COLUMN OF THE SECTION. C C INPUT: C IC1... INDEX OF THE GIVEN COLUMN IN THE MODEL SECTION. C NC1... NUMBER OF COLUMNS IN THE MODEL SECTION. C C NO OUTPUT. C C COMMON BLOCK /COLUMN/: REAL ZL(100),ZM(100),ZR(100),COLL,COLM,COLR INTEGER INTL,INTM,INTR COMMON/COLUMN/ZL,ZM,ZR,COLL,COLM,COLR,INTL,INTM,INTR C C DATE: 1994, JANUARY 26 C CODED BY IVAN PSENCIK, AND LUDEK KLIMES C C----------------------------------------------------------------------- C INTEGER I C C....................................................................... C COLL=AMAX1((FLOAT(IC1)-1.5)/FLOAT(NC1),0.) COLM= (FLOAT(IC1)-0.5)/FLOAT(NC1) COLR=AMIN1((FLOAT(IC1)+0.5)/FLOAT(NC1),1.) IF(IC1.EQ.1) THEN INTL=0 ELSE INTL=INTR DO 1 I=1,INTL ZL(I)=ZR(I) 1 CONTINUE END IF INTM=0 INTR=0 RETURN END C C======================================================================= C SUBROUTINE CONT2(LU2,IFUNC,IBACK,ISB1,ICB1,ISB2,ICB2,STEP,ERR,YP) INTEGER LU2,IFUNC,IBACK,ISB1,ICB1,ISB2,ICB2 REAL STEP,ERR,YP(2) C C SUBROUTINE DESIGNED TO TRACE AN ISOLINE BY MEANS OF NUMERICAL C INTEGRATION, WITHIN THE GIVEN COLUMN. C C INPUT: C LU2... LOGICAL UNIT NUMBER OF THE OUTPUT FILE WITH GENERATED C ISOLINES. C STEP... STEP OF THE NUMERICAL INTEGRATION WHEN COMPUTING THE C INTERFACES OR ISOLINES. C ERR... UPPER ERROR BOUND OF POINTS OF THE ISOLINE DETERMINED BY C MEANS OF NUMERICAL INTEGRATION. C IFUNC...EITHER: C MINUS THE INDEX OF THE FUNCTION DESCRIBING A SURFACE C COVERING A STRUCTURAL INTERFACE, OR C -101, -102, -103, -104, -105 OR -106 FOR THE BOUNDARIES C OF THE MODEL, OR C THE INDEX OF THE COMPLEX BLOCK IN WHICH AN ISOLINE IS C PLOTTED. C IBACK...1 OR 2, INDEX OF THE DIRECTION IN WHICH THE ISOLINE IS TO C BE FOLLOWED. C ISB1,ICB1... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK IN WHICH THE INITIAL POINT YP OF THE ISOLINE IS C SITUATED. C ISB2,ICB2... INDEX OF THE SIMPLE BLOCK AND INDEX OF THE COMPLEX C BLOCK, TOUCHING THE COMPLEX BLOCK ICB1 FROM THE OTHER SIDE C OF THE SURFACE ISRF=-IFUNC AT THE INITIAL POINT YP. C NEED NOT BE DEFINED FOR A VELOCITY ISOLINE (IFUNC C POSITIVE). C STEP... STEP OF THE NUMERICAL INTEGRATION WHEN COMPUTING THE C INTERFACES OR ISOLINES. C ERR... UPPER ERROR BOUND OF POINTS OF THE ISOLINE DETERMINED BY C MEANS OF NUMERICAL INTEGRATION. C YP... ARRAY CONTAINING TWO NORMALIZED SECTION COORDINATES OF THE C INITIAL POINT OF THE ISOLINE TO BE CALCULATED. C C NO OUTPUT. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /CONTC/: INTEGER IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD COMMON/CONTC/ IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL RKGS * EXTERNAL HPCG EXTERNAL FCTI,OUTI C C DATE: 1994, JANUARY 26 C CODED BY IVAN PSENCIK, AND LUDEK KLIMES C C----------------------------------------------------------------------- C INTEGER IHLF,I REAL YPOC(2),DERY(2),AUX(16,2),PRMT(6) C C....................................................................... C IFUN=IFUNC PRMT(1)=0. PRMT(2)=999999. PRMT(3)=STEP PRMT(4)=ERR PRMT(6)=FLOAT(LU2)+.5 C NBACK=IBACK ISB1O=ISB1 ICB1O=ICB1 ISB2O=ISB2 ICB2O=ICB2 NBOD=0 YPOC(1)=YP(1) YPOC(2)=YP(2) C FOR RKGS: DERY(1)=.5 DERY(2)=.5 CALL RKGS(PRMT,YPOC,DERY,2,IHLF,FCTI,OUTI,AUX) C FOR HPCG (PRMT(4)=13.444*ERR): * DERY(1)=.03719 * DERY(2)=.03719 * CALL HPCG(PRMT,YPOC,DERY,2,IHLF,FCTI,OUTI,AUX) C C WRITING DATA DESCRIBING ISOLINES INTO THE FILE LU2 IF(NBOD.GT.1)THEN IF(PRMT(5).GT.-100..AND.PRMT(5).LT.100.) THEN C ISOLINE TERMINATES AT AN INTERFACE C I=INT(SIGN(ABS(PRMT(5))+.5,PRMT(5))) I=INT( ABS(PRMT(5)) ) WRITE(LU2,'(A,I3,A)') '/ (TERMINATING AT SURFACE',I,')' ELSE WRITE(LU2,'(A)') '/' END IF C TERMINAL LINE OF FILE, OVERWRITTEN USING BACKSPACE IN THE CASE C OF OUTPUT WRITE(LU2,'(A)') '/' END IF C RETURN END C C======================================================================= C SUBROUTINE FCTI(X,Y,DERY) REAL X,Y(*),DERY(*) C C SUBROUTINE EVALUATING THE RIGHT-HAND SIDES OF THE ISOLINE TRACING C EQUATIONS. C C INPUT: C X... VALUE OF THE INDEPENDENT VARIABLE ALONG THE ISOLINE. C Y... ARRAY CONTAINING TWO NORMALIZED SECTION COORDINATES OF A C POINT OF THE ISOLINE, DETERMINED BY MEANS OF NUMERICAL C INTEGRATION. C C OUTPUT: C Y... ARRAY CONTAINING TWO NORMALIZED SECTION COORDINATES OF THE C OF THE ISOLINE, CORRECTED BY MEANS OF THE LINEARISATION IN C THE DIRECTION OF THE GRADIENT (PERPENDICULAR TO THE C ISOLINE). C DERY... ARRAY CONTAINING DERIVATIVES OF THE NORMALIZED COORDINATES C Y WITH RESPECT TO X. C C COMMON BLOCK /VALUES/: INTEGER IPS,MV,NV,IV PARAMETER (MV=128) REAL VALUE(MV) COMMON/VALUES/ IPS,VALUE,NV,IV C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /CONTC/: INTEGER IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD COMMON/CONTC/ IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD C C DATE: 1994, JANUARY 26 C CODED BY IVAN PSENCIK, AND LUDEK KLIMES C C----------------------------------------------------------------------- C REAL S(3),S1,S2,AUX C C....................................................................... C CALL FUNC(IFUN,Y,S) S1=S(2) S2=S(3) AUX=SQRT(S1*S1+S2*S2) DERY(1)=S2/AUX DERY(2)=-S1/AUX IF (NBACK.EQ.2) THEN DERY(1)=-DERY(1) DERY(2)=-DERY(2) END IF C C CORRECTION OF THE ISOLINE IF(IFUN.GT.0) THEN AUX=(S(1)-VALUE(IV))/AUX/AUX ELSE AUX=S(1)/AUX/AUX END IF Y(1)=Y(1)-S1*AUX Y(2)=Y(2)-S2*AUX C RETURN END C C======================================================================= C SUBROUTINE OUTI(X,Y,DERY,IHLF,NDIM,PRMT) INTEGER IHLF,NDIM REAL X,Y(NDIM),DERY(NDIM),PRMT(*) C C SUBROUTINE DESIGNED TO CHECK FOR THE INTERSECTIONS OF THE ISOLINE WITH C STRUCTURAL INTERFACES OR BOUNDARIES OF THE COLUMN IN WHICH THE ISOLINE C IS TRACED. C C INPUT: C X... VALUE OF THE INDEPENDENT VARIABLE ALONG THE ISOLINE. C Y... ARRAY CONTAINING TWO NORMALIZED SECTION COORDINATES OF A C POINT OF THE ISOLINE. C DERY... ARRAY CONTAINING DERIVATIVES OF THE NORMALIZED COORDINATES C Y WITH RESPECT TO X. C IHLF... NUMBER OF BISECTIONS OF THE INITIAL INCREMENT OF THE C INDEPENDENT VARIABLE. C NDIM... NUMBER OF ORDINARY DIFFERENTIAL EQUATIONS. C PRMT... ARRAY CONTAINING PARAMETERS OF THE INTEGRATION. C PRMT(6) CONTAINS THE LOGICAL UNIT NUMBER OF THE OUTPUT C FILE WITH GENERATED ISOLINES, INCREASED BY 0.5. C C OUTPUT: C X,Y,DERY... VALUES CORRESPONDING TO THE POINT OF INTERSECTION OF C THE ISOLINE WITH THE BOUNDARY OF THE COMPLEX BLOCK OR THE C COMPUTATIONAL REGION (COLUMN OF THE SECTION), IF THE C ISOLINE INTERSECTS THE BOUNDARY. UNALTERED IF THE ISOLINE C DOES NOT INTERSECT THE BOUNDARY (I.E. IF PRMT(5) REMAINS C UNCHANGED). C PRMT(5)=301,302... THE ISOLINE HAS ALREADY BEEN DETERMINED AND C WILL NOT BE TRACED AGAIN. C 201,202,203,204... THE ISOLINE HAS INTERSECTED THE C BOUNDARY OF THE COMPUTATIONAL REGION (COLUMN OF THE C SECTION). C ISRF... INDEX OF THE INTERSECTED SURFACE LIMITING THE C COMPLEX BLOCK IF A STRUCTURAL INTERFACE HAS BEEN C CROSSED. C OTHERWISE UNALTERED. C C COMMON BLOCK /COLUMN/: REAL ZL(100),ZM(100),ZR(100),COLL,COLM,COLR INTEGER INTL,INTM,INTR COMMON/COLUMN/ZL,ZM,ZR,COLL,COLM,COLR,INTL,INTM,INTR C C COMMON BLOCK /CONTC/: INTEGER IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD COMMON/CONTC/ IFUN,NBACK,ISB1O,ICB1O,ISB2O,ICB2O,NBOD C C DATE: 1994, JANUARY 26 C CODED BY IVAN PSENCIK, AND LUDEK KLIMES C C----------------------------------------------------------------------- C LOGICAL LEFT INTEGER LINE,ISB1,ICB1,ISRF,ISRFO,I,LU2 REAL ERR,ERROR,BNDL,BNDR REAL XOLD,YOLD(2),DOLD(2) REAL COOR(3),DCOOR(3) SAVE LEFT,BNDL,BNDR,XOLD,YOLD,DOLD C C....................................................................... C ERR=PRMT(4) ERROR=PRMT(4) * ERROR=10.*PRMT(4) C IF(NBOD.EQ.0) THEN IF(DERY(1).LT.0.) THEN C THE NEXT POINT OF AN ISOLINE OR A DISCONTINUITY WILL BE C SITUATED TO THE LEFT FROM THE FIRST POINT. THE FIRST POINT IS C TESTED WHETHER AN ISOLINE OR A DISCONTINUITY HAS BEEN C CONSTRUCTED THROUGH IT LEFT=.TRUE. BNDL=COLL BNDR=COLM DO 6 I=1,INTL IF(ABS(ZL(I)-Y(2)).LT.ERROR) THEN PRMT(5)=301. RETURN END IF 6 CONTINUE ELSE C THE NEXT POINT OF AN ISOLINE OR A DISCONTINUITY WILL BE C SITUATED TO THE RIGHT FROM THE FIRST POINT. THE FIRST POINT IS C TESTED WHETHER AN ISOLINE OR A DISCONTINUITY HAS BEEN C CONSTRUCTED THROUGH IT LEFT=.FALSE. BNDL=COLM BNDR=COLR DO 7 I=1,INTM IF(ABS(ZM(I)-Y(2)).LT.ERROR) THEN PRMT(5)=302. RETURN END IF 7 CONTINUE END IF ELSE C C CHECK FOR CROSSING THE BOUNDARY OF THE COMPLEX BLOCK ISB1=ISB1O ICB1=ICB1O IF(IFUN.GE.0) THEN LINE=0 ELSE LINE=-IFUN END IF CALL DISC(XOLD,YOLD,DOLD,X,Y,DERY,ERR,BNDL,BNDR,0.,1., * LINE,ISB1,ICB1,ISRF,ISB1O,ICB1O) IF(LINE.NE.0)THEN ISRFO=ISRF ISB1=ISB2O ICB1=ICB2O CALL DISC(XOLD,YOLD,DOLD,X,Y,DERY,ERR,BNDL,BNDR,0.,1., * LINE,ISB1,ICB1,ISRF,ISB2O,ICB2O) IF(ISRF.EQ.0) THEN ISRF=ISRFO END IF END IF C C TEST WHETHER THE ISOLINE OR DISCONTINUITY INTERSECTS MIDDLE C VERTICAL GRID LINE COLM: IF(LEFT)THEN IF(ISRF.EQ.202)THEN C STORING THE POINT OF INTERSECTION OF THE ISOLINE OR C INTERFACE WITH THE MIDDLE VERTICAL GRID LINE INTL=INTL+1 ZL(INTL)=Y(2) END IF ELSE IF(ISRF.EQ.201)THEN C STORING THE POINT OF INTERSECTION OF THE ISOLINE OR C INTERFACE WITH THE MIDDLE VERTICAL GRID LINE INTM=INTM+1 ZM(INTM)=Y(2) END IF C C TEST WHETHER THE ISOLINE OR DISCONTINUITY INTERSECTS RIGHT C VERTICAL GRID LINE COLR: IF(ISRF.EQ.202)THEN C STORING THE POINT OF INTERSECTION OF THE ISOLINE OR C INTERFACE WITH THE RIGHT VERTICAL GRID LINE INTR=INTR+1 ZR(INTR)=Y(2) END IF END IF C C END OF CHECKS FOR CROSSING BOUNDARIES PRMT(5)=FLOAT(ISRF) END IF C XOLD=X YOLD(1)=Y(1) YOLD(2)=Y(2) DOLD(1)=DERY(1) DOLD(2)=DERY(2) C NBOD=NBOD+1 CALL SECT2(Y,DERY,COOR,DCOOR) LU2=INT(PRMT(6)) WRITE(LU2,'(3F12.6,A)') COOR(1),COOR(2),COOR(3),' /' IF(NBOD.GT.1000)THEN PAUSE 'WARNING: MORE THEN 1000 POINTS OF THE ISOLINE' END IF RETURN END C C======================================================================= C