C Program 'sec.for' to determine the interfaces and velocity isolines in
C 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     SEC...  Main program.  It just invokes subroutine MODSEC.
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 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               points of intersection of interfaces and isosurfaces
C               with the given lines.  The lines are defined as the axes
C               of 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               STEP=-1 (strictly STEP.GT.-1.5):
C                 File with points will contain only the coordinates of
C                 points.
C               STEP=-2 (strictly STEP.LE.-1.5):
C                 File with points will contain the coordinates of
C                 points and the gradients of functions describing the
C                 interface (the gradients are normal to the interface).
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     'LINES SITUATED AT INTERFACES OR VELOCITY ISOSURFACES:'
C     'character*78 string describing the model'
C     /
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) None to several strings terminated by / (a slash).
C     In fact, 2 strings and a slash are generated:
C     'POINTS SITUATED AT INTERFACES OR VELOCITY ISOSURFACES:'
C     'character*78 string describing the model'
C     /
C (2) Written only for STEP.LE.-1.5:
C     D11,D21,D31,D12,D22,D32,D13,D23,D33
C             Vectorial discretization steps:
C             D11=C11/NC1, D21=C11/NC1, D31=C11/NC1,
C             D12=C12/NC2, D22=C12/NC2, D32=C12/NC2,
C             D13=C13    , D23=C13    , D33=C13    .
C (3) For each point of intersection of the vertical line (representing
C     the axis of the vertical column) the following line:
C (3.1) 'SECTnnnnSURFisrf',X1,X2,X3,/          (for STEP.GT.-1.5)
C       'SECTnnnnSURFisrf',X1,X2,X3,F1,F2,F3,/ (for STEP.LE.-1.5)
C     'SECTnnnnSURFisrf'... Character*16 string composed of
C             substring 'SECT', index nnnn of the section, substring
C             'SURF', and of index isrf of the surface covering the
C             structural interface at the point of intersection.
C             The indices are expressed in the format (I4).
C             String 'SECTnnnnSURFisrf' is enclosed in apostrophes.
C     X1,X2,X3... Cartesian coordinates of the point of intersection.
C (4) / (a slash) or end of file.  In fact, a slash is generated.
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
      INCLUDE 'sec.inc'
      END
C     ------------------------------------------------------------------
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C=======================================================================
C
      PROGRAM SEC
C
      CALL MODSEC
      STOP
      END
C
C To convert this program into subroutine package, enter the asterisks
C to the first column of the four preceding FORTRAN statements.
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/:
      INCLUDE 'model.inc'
C None of the storage locations of the common block are altered.
C
C Common blocks /SECTC/ and /VALUES/:
      INCLUDE 'sec.inc'
C None of the storage locations of the common blocks 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
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     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/:
      INCLUDE 'model.inc'
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: 1995, October 29
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(.TRUE.,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=XB
      CALL SECT3(YB,DB,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/:
      INCLUDE 'sec.inc'
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/:
      INCLUDE 'sec.inc'
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/:
      INCLUDE 'model.inc'
C None of the storage locations of the common block are altered.
C
C Common blocks /SECTC/ and /VALUES/:
      INCLUDE 'sec.inc'
C The storage locations of the common blocks 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: 1996, September 30
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 D11,D21,D31,D12,D22,D32
      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
        ELSE
          IF(STEP.LE.-1.5) THEN
            IF(NREPET.LT.0) THEN
              PAUSE 'Error: NREPET negative'
              STOP
            END IF
            D11=C11/FLOAT(NC1)
            D21=C21/FLOAT(NC1)
            D31=C31/FLOAT(NC1)
            D12=C12/FLOAT(NC2)
            D22=C22/FLOAT(NC2)
            D32=C32/FLOAT(NC2)
            WRITE(LU2,'(  9(F10.6,1X)  )')
     *                          D11,D21,D31,D12,D22,D32,C13,C23,C33
          END IF
        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/:
      INCLUDE 'sec.inc'
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/:
      INCLUDE 'sec.inc'
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/:
      INCLUDE 'sec.inc'
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: 1995, March 31
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),VD(10)
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     VD...   Temporary storage to evaluate normal to the interface.
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       ................................................................
        WRITE(*,'(9(A,I4,5X))') '+Section',ISECT
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
*-*         WRITE(*,'(9(A,I4,5X))')
*-*  *                     '+SECTION',ISECT,'COLUMN',IC1
            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.AND.IABS(ISRF).LE.100) 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:
                      IF(STEP.GT.-1.5) THEN
                        WRITE(LU2,'(2(A,I4),A,3F12.6,A)')
     *                  '''SECT',ISECT,', SURF',IFUNC,'''',
     *                                      COOR(1),COOR(2),COOR(3),' /'
                      ELSE
                        CALL SRFC2(IFUNC,COOR,VD)
***                     WRITE(LU2,'(2(A,I4),A,6F12.6,A)')
***  *                  '''SECT',ISECT,', SURF',IFUNC,'''',
                        WRITE(LU2,'(        A,6F12.6,A)')
     *                  ''' ''',
     *                    COOR(1),COOR(2),COOR(3),VD(2),VD(3),VD(4),' /'
                      END IF
                    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/:
      INCLUDE 'sec.inc'
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 blocks /VALUES/ and /CONTC/:
      INCLUDE 'sec.inc'
C None of the storage locations of common block /VALUES/ are altered.
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 blocks /VALUES/ and /CONTC/:
      INCLUDE 'sec.inc'
C None of the storage locations of common /VALUES/ block are altered.
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 blocks /COLUMN/ and /CONTC/:
      INCLUDE 'sec.inc'
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
      INCLUDE 'model.for'
      INCLUDE 'metric.for'
      INCLUDE 'srfc.for'
      INCLUDE 'parm.for'
      INCLUDE 'val.for'
      INCLUDE 'fit.for'
      INCLUDE 'means.for'
      INCLUDE 'rkgs.for'
C
C=======================================================================
C