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
