C 
C=======================================================================
C
      PROGRAM MODSRF
C for coverage of structural interfaces by polygons and for computation
C of 2-D slices through a model.
C
C Version: 5.30
C Date: 1999, May 23
C
C----------------------------------------------------------------------
C
C Structural interfaces:
C Program MODSRF covers structural interfaces by polygons composed
C of points with known cartesian coordinates. Structural interfaces
C are in MODEL package defined by implicit functions. This is very
C useful for computation purposes, but not very useful for
C visualization of models. For visualization purposes the interfaces
C must be expressed in explicit form, e.g. as sets of polygons
C composed of points with known coordinates.
C Method: The model is decomposed into a set of cubes given by input
C parameters Oi, Ni, Di. Points of intersections of cube edges with
C structural interfaces are computed together with points of
C intersection of structural edges with faces of the cubes. Points
C which belong to the same interface are then connected into
C polygons for each cube.
C This mode is started if all N1, N2, N3 are greater than 1.
C
C 2-D sections:
C If one of the values N1, N2, N3 equal 1, the 2-D sections are
C generated. Above mentioned points of intersection are computed
C along rectangles given by Oi, Ni, Di. Points which belong to the
C same complex block are then connected into polygons for each
C rectangle.
C
C
C Coded by Petr Bulant
C             Department of Geophysics, Charles University
C             Ke Karlovu 3,  121 16  Praha 2,  Czech Republic
C             E-mail: bulant@seis.karlov.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of input and output data:
C
C Main input data read from the * input device:
C     The data are read in by the list directed input (free format), by
C     a single READ statement.  The data consist of character strings
C     which have to be enclosed in apostrophes.  Note that the
C     interactive * external unit may be piped or redirected to a file.
C     The following items are read:
C (1) 'SEP' /
C     'SEP'..String in apostrophes containing the name of the main
C            input data file with parameters in the SEP format.
C            The parameters are: the filename of the data for the model,
C            dimensions of the grid of N1*N2*N3 gridpoints, the form
C            and names of the output files with points and polygons,
C            etc.
C     /...   It is recommended to terminate these data by a slash.
C     Default: 'SEP'='len-srf.h'.
C
C                                                     
C Input data file 'SEP':
C Data file 'SEP' has the form of the SEP (Stanford Exploration
C Project) parameter file:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
C Data specifying the model by means of the MODEL package:
C     MODEL='string'... Name of the input file with the data specifying
C             the model.  Must be specified.
C             Description of MODEL.
C             Example of data set MODEL.
C             Default: MODEL=' '
C Parameters defining the basic regular rectangular grid:
C     N1=positive integer... Number of gridpoints of the basic grid
C             along the X1 axis.  Default: N1=1
C     N2=positive integer... Number of gridpoints of the basic grid
C             along the X2 axis.  Default: N2=1
C     N3=positive integer... Number of gridpoints of the basic grid
C             along the X3 axis.  Default: N3=1
C       The values of N1, N2, N3 also decide about the mode in which
C       MODSRF runs:
C         all of N1,N2,N3 greater 1: computation of polygons along
C                                    structural interfaces.
C         one of N1,N2,N3 equal 1: computation of polygons along
C                                  2-D section through the model.
C     O1=real... X1 coordinate of the origin of the grid. Default: O1=0.
C     O2=real... X2 coordinate of the origin of the grid. Default: O2=0.
C     O3=real... X3 coordinate of the origin of the grid. Default: O3=0.
C     D1=real... Grid spacing along the X1 axis.  Default: D1=1.
C     D2=real... Grid spacing along the X2 axis.  Default: D2=1.
C     D3=real... Grid spacing along the X3 axis.  Default: D3=1.
C       The grid intervals may also be negative.
C Names of the output files:
C     VRTX='string'... Name of the output file with vertices of the
C             polygons. Description of file VRTX.
C             Default: VRTX='vrtx.out'
C     PLGN='string'... Name of the output file describing the polygons.
C             Description of file PLGN.
C             Default: PLGN='plgn.out'
C     PLGNS='string'... Name of the file describing the polygons in
C             terms of the names of the vertices.
C             Description of file PLGNS.
C             Default: PLGNS=' '
C Parameters specifying the quantities to be written into the file VRTX.
C     TEXTP='string' ... First part of names of vertices. The second
C             part of the name of a vertex is formed by number giving
C             its position in the file VRTX.
C             Default: TEXTP='V'
C     COLUMN07 to COLUMN69, POWER07 to POWER69:
C     POWERii=real ... Power of the quantity to be written in column ii.
C     COLUMNii='string' ... String which specifies the quantity to be
C             written to the column ii of the file VRTX. First six
C             columns are reserved for coordinates of the vertices and
C             for the normals. Column zero is reserved for names of the
C             vertices. Following strings are allowed:
C             ' ' (a space) ... Nothing is to be written to the column
C                      ii and to all the following columns.
C             ISRF ... Index of the interface.
C             +ISB ... Index of simple block at positive side of the
C                      interface.
C             -ISB ... Index of simple block at negative side of the
C                      interface.
C             +ICB -ICB ... Index of complex block.
C             +VP  -VP  ... P-wave velocity.
C             +VS  -VS  ... S-wave velocity.
C             +DEN -DEN ... Density.
C             +QP  -QP  ... P wave loss factor.
C             +QS  -QS  ... S wave loss factor.
C             +A11 -A11 +A12 -A12 +A22 -A22 +A13 -A13 +A23 -A23
C             +A33 -A33 +A14 -A14 +A24 -A24 +A34 -A34 +A44 -A44
C             +A15 -A15 +A25 -A25 +A35 -A35 +A45 -A45 +A55 -A55
C             +A16 -A16 +A26 -A26 +A36 -A36 +A46 -A46 +A56 -A56
C             +A66 -A66 ... Reduced (i.e. divided by the density)
C               anisotropic elastic parameters
C               (components of the real part of the symmetric 6*6
C               stiffness matrix divided by the density).
C             +Q11 -Q11 +Q12 -Q12 +Q22 -Q22 +Q13 -Q13 +Q23 -Q23
C             +Q33 -Q33 +Q14 -Q14 +Q24 -Q24 +Q34 -Q34 +Q44 -Q44
C             +Q15 -Q15 +Q25 -Q25 +Q35 -Q35 +Q45 -Q45 +Q55 -Q55
C             +Q16 -Q16 +Q26 -Q26 +Q36 -Q36 +Q46 -Q46 +Q56 -Q56
C             +Q66 -Q66 ... Reduced (i.e. divided by the density)
C               imaginary anisotropic elastic parameters
C               (components of the imaginary part of the
C               symmetric 6*6 stiffness matrix divided by the density).
C             All strings may be entered either in uppercase or in
C             lowercase letters.
C             Default: COLUMN07='ISRF', COLUMN08 to COLUMN69=' ',
C                      POWER07 to POWER69=1.
C
C                                                    
C Output file VRTX with the vertices:
C (1) / ... a slash.
C (2) For each vertex data (2.1):
C (2.1) 'NAME',X1,X2,X3,Z1,Z2,Z3,/
C     'NAME'... Name of the vertex. See parameter TEXTP above.
C     X1,X2,X3... Coordinates of the vertex.
C     Z1,Z2,Z3... Normal to the interface or to the 2-D section at
C                 the vertex.
C     /...    None to several values terminated by a slash. See
C             parameters COLUMN07 to COLUMN69 above.
C (3) / ... a slash at the end of file.
C
C                                                    
C Output file PLGN with the polygons specified in terms of
C indices of vertices:
C (1) For each polygon data (1.1):
C (1.1) I1,I2,...,IN,/
C     I1,I2,...,IN... Indices of N vertices of the polygon.
C             The vertices in file VRTX are indexed by positive integers
C             according to their order.
C     /...    List of vertices is terminated by a slash.
C (2) / ... a slash at the end of file.
C
C                                                   
C Output file PLGNS with the polygons specified in terms of
C names of vertices. This enables vertices from several files
C to be combined in further application (the names of the vertices must
C differ, see the parameter TEXTP above).
C (1) For each polygon data (1.1):
C (1.1) 'VRTX1','VRTX2',...,'VRTXN',/
C     'VRTX1','VRTX2',...,'VRTXN'... Strings containing the names of N
C             vertices of the polygon.  The names correspond to the
C             names in file VRTX.
C     /...    List of vertices is terminated by a slash.
C (2) / ... a slash at the end of file.
C
C
C Subroutines and external functions required:
      EXTERNAL RKGS
      EXTERNAL FCTMS,OUTMS
      LOGICAL MSLPIF,MSLPIL
      EXTERNAL MSLPIF,MSLPIL
      INTEGER LENGTH
      EXTERNAL LENGTH
C
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Input and output data files:
      CHARACTER*80 FSEP,FMOD
      INTEGER LU,LU1,LU2
      PARAMETER (LU=1,LU1=1,LU2=2)
C     Auxiliary storage locations:
      INTEGER I1,I2,I3,I4,I5,I6,I7,J1,J2,J3,ITER,NPOINT,LEN1,LEN2,LENG
      INTEGER IL1,IL2,IL3,IL4,IF1,IF2,IF3,IF4,IF5,IF6
      INTEGER IX,IP1,IP2,ISHIFT
      INTEGER ISBOLD,ICBOLD,ISBNEW,ICBNEW,ISRF,ICBPOS,ICBNEG
      REAL X1,X2,X3,XX(3),Y1,Y2,Y3,YY(3)
      EQUIVALENCE (X1,XX(1)),(X2,XX(2)),(X3,XX(3))
      EQUIVALENCE (Y1,YY(1)),(Y2,YY(2)),(Y3,YY(3))
      REAL F(10)
      INTEGER IDUMMY(1),IY(8)
      REAL ERR,DUMMY(1),XOLD1,XOLD(3),DOLD(3),XNEW1,XNEW(3),DNEW(3),
     *     XTMP1,XTMP(3),DTMP(3),XINT1,XINT(3),DINT(3)
      REAL PRMT(5),AUX(8,3)
      INTEGER MJPT,NJPT
      PARAMETER (MJPT=8)
      INTEGER JPT(6,MJPT)
C     JPT(1 to 6,i) ... Points along the gridface: address of the point,
C       ICB,ISB,ISRF,ISB,ICB. JPT(1,i).LT.0 means that the point is
C       stored in reverse order than it is in IRAM.
      INTEGER MJCN,NJCN,NNJCN(0:7),MJSP,NJSP,
     * MJPOL,NJPOL,MLJPOL,NLJPOL,NIND
      PARAMETER (MJCN=50,MJSP=100,MJPOL=10,MLJPOL=10)
      INTEGER JCN(3,MJCN),JSP(3,MJSP),JPOL(0:MLJPOL,MJPOL),IND(MJCN)
C     JCN(1 to 3,i) ... Connection i of the face: address of the first
C       point, address of the second point, index of the interface.
C     NNJCN(1:7) ... Number of connections for first gridface, second
C     gridface, ... , sixth gridface, between edges.
C     JSP(1 to 3,i) ... Array of starting points as they go along
C       a polygon: previous point, point, following point.
C     JPOL(0,i) ... Number of points in polygon i.
C     JPOL(1 to JPOL(0,i),i) ... Points in the polygon.
      INTEGER MIPTE
      PARAMETER (MIPTE=5)
      INTEGER IPTE(MIPTE),NIPTE
C     IPTE(i) ... address of the point to be used as starting point
C       when searching for the edge. IPTE(i).LT.0 means that the point
C       is to be used in reverse order than it is stored in array IRAM.
      REAL Z(69),POWER(69),VP(10),VS(10),RHO,QP,QS,A(10,21),Q(21),
     *     OUTMIN,OUTMAX
      LOGICAL LCN
      CHARACTER*20 FORMAT,FORMA1,FORMA2,VRTX,PLGN,PLGNS,TEXTP,
     *             TEXTC(69)*4,TEXTS(MLJPOL)
C
C.......................................................................
C
C     Reading main input data:
      FSEP='len-srf.h'
      WRITE(*,'(A)')
     *' MODSRF: Enter the name of SEP parameter file: '
      READ (*,*) FSEP
      WRITE(*,'(A)')
     *'+MODSRF: Reading input data.                                   '

C     Reading all the data from file FSEP to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      IF (D1.LT.0.) THEN
        O1=O1+(N1-1)*D1
        D1=-D1
      ENDIF
      IF (D2.LT.0.) THEN
        O2=O2+(N2-1)*D2
        D2=-D2
      ENDIF
      IF (D3.LT.0.) THEN
        O3=O3+(N3-1)*D3
        D3=-D3
      ENDIF
      IF ((N1.LE.0).OR.(N2.LE.0).OR.(N3.LE.0).OR.
     *    ((D1.EQ.0.).AND.(N1.NE.1)).OR.
     *    ((D2.EQ.0.).AND.(N2.NE.1)).OR.
     *    ((D3.EQ.0.).AND.(N3.NE.1))) THEN
C       MODSRF-01
        CALL ERROR('MODSRF-01: Wrong specification of the grid.')
C       This specification of the grid may cause problems. Please,
C       specify D1,D2,D3 nonzero and N1,N2,N3 greater than 0. Di may
C       equal zero in case that corresponding Ni equals 1.
      ENDIF
C
C     Reading the data for the model:
      CALL RSEP3T('MODEL',FMOD,' ')
      OPEN(LU,FILE=FMOD,STATUS='OLD')
      CALL MODEL1(LU)
      CLOSE(LU)
C
C     Recalling the output filenames:
      CALL RSEP3T('VRTX',VRTX,'vrtx.out')
      CALL RSEP3T('PLGN',PLGN,'plgn.out')
      CALL RSEP3T('PLGNS',PLGNS,' ')
C
C     Recalling the first part of names of points in output file VRTX:
      CALL RSEP3T('TEXTP',TEXTP,'V')
C
C     Recalling the data specifying the quantities to be written into
C     the output file with points at structural interfaces:
      FORMA1='COLUMN00'
      FORMA2='POWER00'
      I1=7
   5  CONTINUE
        FORMA1(8:8)=CHAR(ICHAR('0')+MOD(I1,10))
        FORMA2(7:7)=CHAR(ICHAR('0')+MOD(I1,10))
        FORMA1(7:7)=CHAR(ICHAR('0')+I1/10)
        FORMA2(6:6)=CHAR(ICHAR('0')+I1/10)
        IF (I1.EQ.7) THEN
          CALL RSEP3T(FORMA1,TEXTC(7),'ISRF')
        ELSE
          CALL RSEP3T(FORMA1,TEXTC(I1),' ')
        ENDIF
        CALL RSEP3R(FORMA2,POWER(I1),1.)
        IF (TEXTC(I1).NE.' ') THEN
          CALL LOWER(TEXTC(I1))
          IF ((TEXTC(I1).NE.'isrf').AND.
     *        (TEXTC(I1).NE.'+isb').AND.(TEXTC(I1).NE.'-isb').AND.
     *        (TEXTC(I1).NE.'+icb').AND.(TEXTC(I1).NE.'-icb').AND.
     *        (TEXTC(I1).NE.'+vp') .AND.(TEXTC(I1).NE.'-vp') .AND.
     *        (TEXTC(I1).NE.'+vs') .AND.(TEXTC(I1).NE.'-vs') .AND.
     *        (TEXTC(I1).NE.'+den').AND.(TEXTC(I1).NE.'-den').AND.
     *        (TEXTC(I1).NE.'+qp') .AND.(TEXTC(I1).NE.'-qp') .AND.
     *        (TEXTC(I1).NE.'+qs') .AND.(TEXTC(I1).NE.'-qs') .AND.
     *        (TEXTC(I1).NE.'+a11').AND.(TEXTC(I1).NE.'-a11').AND.
     *        (TEXTC(I1).NE.'+a12').AND.(TEXTC(I1).NE.'-a12').AND.
     *        (TEXTC(I1).NE.'+a22').AND.(TEXTC(I1).NE.'-a22').AND.
     *        (TEXTC(I1).NE.'+a13').AND.(TEXTC(I1).NE.'-a13').AND.
     *        (TEXTC(I1).NE.'+a23').AND.(TEXTC(I1).NE.'-a23').AND.
     *        (TEXTC(I1).NE.'+a33').AND.(TEXTC(I1).NE.'-a33').AND.
     *        (TEXTC(I1).NE.'+a14').AND.(TEXTC(I1).NE.'-a14').AND.
     *        (TEXTC(I1).NE.'+a24').AND.(TEXTC(I1).NE.'-a24').AND.
     *        (TEXTC(I1).NE.'+a34').AND.(TEXTC(I1).NE.'-a34').AND.
     *        (TEXTC(I1).NE.'+a44').AND.(TEXTC(I1).NE.'-a44').AND.
     *        (TEXTC(I1).NE.'+a15').AND.(TEXTC(I1).NE.'-a15').AND.
     *        (TEXTC(I1).NE.'+a25').AND.(TEXTC(I1).NE.'-a25').AND.
     *        (TEXTC(I1).NE.'+a35').AND.(TEXTC(I1).NE.'-a35').AND.
     *        (TEXTC(I1).NE.'+a45').AND.(TEXTC(I1).NE.'-a45').AND.
     *        (TEXTC(I1).NE.'+a55').AND.(TEXTC(I1).NE.'-a55').AND.
     *        (TEXTC(I1).NE.'+a16').AND.(TEXTC(I1).NE.'-a16').AND.
     *        (TEXTC(I1).NE.'+a26').AND.(TEXTC(I1).NE.'-a26').AND.
     *        (TEXTC(I1).NE.'+a36').AND.(TEXTC(I1).NE.'-a36').AND.
     *        (TEXTC(I1).NE.'+a46').AND.(TEXTC(I1).NE.'-a46').AND.
     *        (TEXTC(I1).NE.'+a56').AND.(TEXTC(I1).NE.'-a56').AND.
     *        (TEXTC(I1).NE.'+a66').AND.(TEXTC(I1).NE.'-a66').AND.
     *        (TEXTC(I1).NE.'+q11').AND.(TEXTC(I1).NE.'-q11').AND.
     *        (TEXTC(I1).NE.'+q12').AND.(TEXTC(I1).NE.'-q12').AND.
     *        (TEXTC(I1).NE.'+q22').AND.(TEXTC(I1).NE.'-q22').AND.
     *        (TEXTC(I1).NE.'+q13').AND.(TEXTC(I1).NE.'-q13').AND.
     *        (TEXTC(I1).NE.'+q23').AND.(TEXTC(I1).NE.'-q23').AND.
     *        (TEXTC(I1).NE.'+q33').AND.(TEXTC(I1).NE.'-q33').AND.
     *        (TEXTC(I1).NE.'+q14').AND.(TEXTC(I1).NE.'-q14').AND.
     *        (TEXTC(I1).NE.'+q24').AND.(TEXTC(I1).NE.'-q24').AND.
     *        (TEXTC(I1).NE.'+q34').AND.(TEXTC(I1).NE.'-q34').AND.
     *        (TEXTC(I1).NE.'+q44').AND.(TEXTC(I1).NE.'-q44').AND.
     *        (TEXTC(I1).NE.'+q15').AND.(TEXTC(I1).NE.'-q15').AND.
     *        (TEXTC(I1).NE.'+q25').AND.(TEXTC(I1).NE.'-q25').AND.
     *        (TEXTC(I1).NE.'+q35').AND.(TEXTC(I1).NE.'-q35').AND.
     *        (TEXTC(I1).NE.'+q45').AND.(TEXTC(I1).NE.'-q45').AND.
     *        (TEXTC(I1).NE.'+q55').AND.(TEXTC(I1).NE.'-q55').AND.
     *        (TEXTC(I1).NE.'+q16').AND.(TEXTC(I1).NE.'-q16').AND.
     *        (TEXTC(I1).NE.'+q26').AND.(TEXTC(I1).NE.'-q26').AND.
     *        (TEXTC(I1).NE.'+q36').AND.(TEXTC(I1).NE.'-q36').AND.
     *        (TEXTC(I1).NE.'+q46').AND.(TEXTC(I1).NE.'-q46').AND.
     *        (TEXTC(I1).NE.'+q56').AND.(TEXTC(I1).NE.'-q56').AND.
     *        (TEXTC(I1).NE.'+q66').AND.(TEXTC(I1).NE.'-q66')) THEN
C           MODSRF-02
            CALL ERROR('MODSRF-02: Wrong value of COLUMN.')
C           See allowed values of COLUMNii in the
C           description of file SEP.
          ENDIF
          I1=I1+1
          IF (I1.GT.69) THEN
C           MODSRF-03
            CALL ERROR('MODSRF-03: More than 69 COLUMNs.')
C           Currently up to 69 values of COLUMNii may be specified.
C           See allowed values of COLUMNii in the
C           description of file SEP.
          ENDIF
          GOTO 5
        ENDIF
C     End of the loop over future columns of the output file.
C
C     Preparing numbers of gridpoints, gridlegs, gridfaces
C     and gridcubes:
      N11=N1-1
      N21=N2-1
      N1N2=   N1   * N2
      N11N2= (N1-1)* N2
      N1N21=  N1   *(N2-1)
      N11N21=(N1-1)*(N2-1)
      NGPS=N1*N2*N3
      OLEG=NGPS*2+1
      NLEG1=(N1-1)*N2*N3
      NLEG2=N1*(N2-1)*N3
      NLEG3=N1*N2*(N3-1)
      NLEG12=NLEG1+NLEG2
      NLEG =NLEG12+NLEG3
      OFAC=OLEG+NLEG+1
      NFAC1=N1*(N2-1)*(N3-1)
      NFAC2=(N1-1)*N2*(N3-1)
      NFAC3=(N1-1)*(N2-1)*N3
      NFAC12=NFAC1+NFAC2
      NFAC=NFAC12+NFAC3
      NCUB=(N1-1)*(N2-1)*(N3-1)
      OPOI=OFAC+NFAC
      NPOI=OPOI
      IF (NPOI.GT.MRAM) CALL MSEROR(1)
C
      DO 10, I1=1,OPOI
        IRAM(I1)=0
  10  CONTINUE
C
C
C     Loop along all gridpoints, recording indices of blocks:
      WRITE(*,'(A)')
     *'+MODSRF: Computing indices of blocks.                          '
      DO 13, I3=1,N3
        DO 12, I2=1,N2
          DO 11, I1=1,N1
            IX=(I3-1)*N1N2+(I2-1)*N1+I1
            IX=2*(IX-1)
            X1=O1+FLOAT(I1-1)*D1
            X2=O2+FLOAT(I2-1)*D2
            X3=O3+FLOAT(I3-1)*D3
            CALL BLOCK(XX,0,0,ISRF,ISBNEW,ICBNEW)
            IRAM(IX+1)=ISBNEW
            IRAM(IX+2)=ICBNEW
  11      CONTINUE
  12    CONTINUE
  13  CONTINUE
C
C
C     Loop along all gridlegs,
C     searching for intersections with structural interfaces:
      WRITE(*,'(A)')
     *'+MODSRF: Computing intersections along gridlegs.               '
C     Address of intersection points on gridleg 0:
      IRAM(OLEG)=NPOI
C     Index of coordinate in the direction of gridlegs:
      I4=1
C     Maximum alllowed error in the position of intersection point:
      ERR=0.001*D1
C     Auxiliary quantities:
      DOLD(1)=1.
      DOLD(2)=0.
      DOLD(3)=0.
      DNEW(1)=1.
      DNEW(2)=0.
      DNEW(3)=0.
      DTMP(1)=1.
      DTMP(2)=0.
      DTMP(3)=0.
      DO 29, I1=1,NLEG
        IF (I1.EQ.NLEG1+1) THEN
          I4=2
          ERR=0.001*D2
          DOLD(1)=0.
          DOLD(2)=1.
          DNEW(1)=0.
          DNEW(2)=1.
          DTMP(1)=0.
          DTMP(2)=1.
        ENDIF
        IF (I1.EQ.NLEG12+1) THEN
          I4=3
          ERR=0.001*D3
          DOLD(2)=0.
          DOLD(3)=1.
          DNEW(2)=0.
          DNEW(3)=1.
          DTMP(2)=0.
          DTMP(3)=1.
        ENDIF
        CALL MSGLEG(I1,IP1,IP2)
        ISBOLD=IRAM(2*(IP1-1)+1)
        ISBNEW=IRAM(2*(IP2-1)+1)
        IF (ISBOLD.NE.ISBNEW) THEN
C         The points of the gridleg are in different simple blocks:
          ICBOLD=IRAM(2*(IP1-1)+2)
          ICBNEW=IRAM(2*(IP2-1)+2)
          CALL MSCP(IP1,XOLD(1),XOLD(2),XOLD(3))
          CALL MSCP(IP2,XNEW(1),XNEW(2),XNEW(3))
          XTMP(1)=XNEW(1)
          XTMP(2)=XNEW(2)
          XTMP(3)=XNEW(3)
C         Loop for intersection points along the gridleg:
          ITER=0
  21      CONTINUE
            ITER=ITER+1
            IF (ITER.GT.100) THEN
C             MODSRF-04
              CALL ERROR ('MODSRF-04: More than 100 intersections.')
C             More than 100 points of intersection of the
C             gridline element with the structural interfaces.
C             Check the input data and then contact the author.
            ENDIF
C
C           Determining the interface between points XOLD, XNEW:
            IY(4)=ISBOLD
            IY(5)=ICBOLD
            IY(6)=0
            XTMP(I4)=XNEW(I4)
            XOLD1=XOLD(I4)
            XNEW1=XNEW(I4)
            XTMP1=XNEW1
            CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERR,
     *               XOLD1,XOLD1,XOLD,DOLD,XNEW1,XNEW,DNEW,
     *                     XTMP1,XTMP,DTMP,XINT1,XINT,DINT)
            IF ((ICBOLD.NE.ICBNEW).AND.(IY(6).EQ.0)) THEN
C             The interface was not found between the points with
C             different ICB. This may happen when the NEW point is
C             situated directly at the interface.
C             Repeating the search in reverse direction:
              DOLD(I4)=-DOLD(I4)
              DNEW(I4)=-DNEW(I4)
              DTMP(I4)=-DTMP(I4)
              XTMP1=XOLD1
              XTMP(1)=XOLD(1)
              XTMP(2)=XOLD(2)
              XTMP(3)=XOLD(3)
              IY(4)=ISBNEW
              IY(5)=ICBNEW
              CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERR,
     *               (XTMP1-ERR),XNEW1,XNEW,DNEW,XOLD1,XOLD,DOLD,
     *                     XTMP1,XTMP,DTMP,XINT1,XINT,DINT)
              IF (IY(6).EQ.0) THEN
C               The interface was not found even in reverse direction.
C               This may happen when the whole gridleg coincides with an
C               interface:
                CALL BLOCK(XOLD,0,ISBOLD,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBOLD) CALL MSEROR(8)
                CALL BLOCK(XOLD,0,ISBNEW,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBNEW) CALL MSEROR(8)
                CALL BLOCK(XNEW,0,ISBNEW,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBNEW) CALL MSEROR(8)
                CALL BLOCK(XNEW,0,ISBOLD,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBOLD) CALL MSEROR(8)
C               Whole gridleg coincides with an interface:
                CALL SEPAR(ISBNEW,ISBOLD,IDUMMY,IY(6))
                IF (IDUMMY(1).LT.1) CALL MSEROR(8)
                IY(7)=ISBOLD
                IY(8)=ICBOLD
              ENDIF
              IF (IY(8).NE.ICBOLD) THEN
C               There should be only one interface between the points.
C               
C               MODSRF-06
                CALL ERROR('MODSRF-06: More than one interface.')
C               This error should not appear. Contact the author.
              ENDIF
C             Reverting the points to the original order:
              DOLD(I4)=-DOLD(I4)
              DNEW(I4)=-DNEW(I4)
              DTMP(I4)=-DTMP(I4)
              XTMP1=XNEW1
              XTMP(1)=XNEW(1)
              XTMP(2)=XNEW(2)
              XTMP(3)=XNEW(3)
              ISBOLD=IY(7)
              ICBOLD=IY(8)
              ISBNEW=IY(4)
              ICBNEW=IY(5)
              IY(4)=ISBOLD
              IY(5)=ICBOLD
              IY(6)=-IY(6)
              IY(7)=ISBNEW
              IY(8)=ICBNEW
            ENDIF
            IF (IY(6).NE.0) THEN
C             Structural interface, recording the point of intersection:
C             XINT is the point of intersection with the interface,
C             IY(4) and IY(5) are the indices of the simple block
C             and the complex block before the interface,
C             IY(6) is the index of the interface, IY(7) and IY(8) are
C             the indices of the simple block and the complex block
C             behind the interface:
              IF (.NOT.MSLPIL(XINT,I1)) THEN
C               MODSRF-07
                CALL ERROR('MODSRF-07: Point not on the gridleg.')
C               This error should not appear. Contact the author.
              ENDIF
              IF (NPOI+NQPOI.GT.MRAM) CALL MSEROR(1)
              RAM(NPOI+1)=XINT(1)
              RAM(NPOI+2)=XINT(2)
              RAM(NPOI+3)=XINT(3)
              IRAM(NPOI+4)=IY(5)
              IRAM(NPOI+5)=IY(4)
              IRAM(NPOI+6)=IY(6)
              IRAM(NPOI+7)=IY(7)
              IRAM(NPOI+8)=IY(8)
              NPOI=NPOI+NQPOI
              IF (IY(7).NE.ISBNEW) THEN
C               IY(7) is the index of simple block behind the interface,
C               the search for intersection points must continue.
C               Shifting point XOLD to the point of intersection:
                XOLD(I4)=XINT(I4)
                XOLD1=XOLD(I4)
                ISBOLD=IY(7)
                ICBOLD=IY(8)
                GOTO 21
              ENDIF
            ENDIF
C         End of the loop for intersection points along the gridleg.
        ENDIF
        IRAM(OLEG+I1)=NPOI
  29  CONTINUE
      NPOIL=NPOI
C
C
C     Loop along all gridfaces,
C     connecting the points of intersection of structural interfaces
C     with gridlegs, if necessary computing points of intersection
C     of structural edges with gridfaces.
      WRITE(*,'(A)')
     *'+MODSRF: Connecting points along gridfaces.                     '
      IF (NCUB.NE.0) THEN
        OCON=NPOI+(NPOI-OPOI)/2
      ELSE
        OCON=NPOI + (NPOI-OPOI)*3 + NGPS*NQPOI
      ENDIF
      NCON=OCON
      IF (NCON.GT.MRAM) CALL MSEROR(1)
C     Initialization for RKGS:
      PRMT(1)=0.
      PRMT(2)=(D1+D2+D3)
      PRMT(3)=PRMT(2)/30.
      PRMT(4)=PRMT(3)/100.
      DTMP(1)=1./3.
      DTMP(2)=DTMP(1)
      DTMP(3)=DTMP(1)
C     Address of connections on gridface 0:
      IRAM(OFAC)=NCON
      DO 49, I1=1,NFAC
C       Indices of gridlegs of the gridface:
        CALL MSGFAC(I1,IL1,IL2,IL3,IL4)
C
C       Forming array with intersection points for the gridface:
        NJPT=0
        DO 31, I2=IRAM(OLEG+IL1-1),IRAM(OLEG+IL1)-NQPOI,NQPOI
C         I2 points to the start of records for intersection point.
          NJPT=NJPT+1
          IF (NJPT.GT.MJPT) CALL MSEROR(2)
          JPT(1,NJPT)=I2+NQPOI
          JPT(2,NJPT)=IRAM(I2+4)
          JPT(3,NJPT)=IRAM(I2+5)
          JPT(4,NJPT)=IRAM(I2+6)
          JPT(5,NJPT)=IRAM(I2+7)
          JPT(6,NJPT)=IRAM(I2+8)
  31    CONTINUE
        DO 32, I2=IRAM(OLEG+IL2-1),IRAM(OLEG+IL2)-NQPOI,NQPOI
C         I2 points to the start of records for intersection point.
          NJPT=NJPT+1
          IF (NJPT.GT.MJPT) CALL MSEROR(2)
          JPT(1,NJPT)=I2+NQPOI
          JPT(2,NJPT)=IRAM(I2+4)
          JPT(3,NJPT)=IRAM(I2+5)
          JPT(4,NJPT)=IRAM(I2+6)
          JPT(5,NJPT)=IRAM(I2+7)
          JPT(6,NJPT)=IRAM(I2+8)
  32    CONTINUE
        DO 33, I2=IRAM(OLEG+IL3)-NQPOI,IRAM(OLEG+IL3-1),-NQPOI
C         I2 points to the start of records for intersection point.
          NJPT=NJPT+1
          IF (NJPT.GT.MJPT) CALL MSEROR(2)
          JPT(1,NJPT)=-(I2+NQPOI)
          JPT(2,NJPT)=IRAM(I2+8)
          JPT(3,NJPT)=IRAM(I2+7)
          JPT(4,NJPT)=-IRAM(I2+6)
          JPT(5,NJPT)=IRAM(I2+5)
          JPT(6,NJPT)=IRAM(I2+4)
  33    CONTINUE
        DO 34, I2=IRAM(OLEG+IL4)-NQPOI,IRAM(OLEG+IL4-1),-NQPOI
C         I2 points to the start of records for intersection point.
          NJPT=NJPT+1
          IF (NJPT.GT.MJPT) CALL MSEROR(2)
          JPT(1,NJPT)=-(I2+NQPOI)
          JPT(2,NJPT)=IRAM(I2+8)
          JPT(3,NJPT)=IRAM(I2+7)
          JPT(4,NJPT)=-IRAM(I2+6)
          JPT(5,NJPT)=IRAM(I2+5)
          JPT(6,NJPT)=IRAM(I2+4)
  34    CONTINUE
C
C       Making connections between points of the gridface:
        DO 46, I2=1,NJPT
          DO 36, I3=I2+1,NJPT
            IF ((JPT(2,I2).EQ.JPT(6,I3)).AND.
     *          (JPT(6,I2).EQ.JPT(2,I3))) THEN
C             Same indices of complex blocks:
              IF (IABS(JPT(4,I2)).EQ.IABS(JPT(4,I3))) THEN
C               Same interface - connection found:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                IF (JPT(4,I2).LT.0) THEN
                  IRAM(NCON+1)=IABS(JPT(1,I2))
                  IRAM(NCON+2)=IABS(JPT(1,I3))
                  IRAM(NCON+3)=JPT(4,I2)
                  NCON=NCON+3
                ELSE
                  IRAM(NCON+1)=IABS(JPT(1,I3))
                  IRAM(NCON+2)=IABS(JPT(1,I2))
                  IRAM(NCON+3)=JPT(4,I3)
                  NCON=NCON+3
                ENDIF
                GOTO 45
              ENDIF
            ENDIF
  36      CONTINUE
          DO 37, I3=IRAM(OFAC+I1-1),NCON-3,3
            IF ((IRAM(I3+1).EQ.IABS(JPT(1,I2))).OR.
     *          (IRAM(I3+2).EQ.IABS(JPT(1,I2))))
C             The point I2 is already connected:
     *        GOTO 45
  37      CONTINUE
C         Point I2 is not connected, looking for point of intersection
C         of structural edges with gridface. (Following the intersection
C         of the interface with the gridface to the edge.)
          IFACE=I1
          IPTE(1)=JPT(1,I2)
          NIPTE=1
C         Loop for identification of all edges connected with point I2:
  39      CONTINUE
C           Initiating the search for the edge:
            I3=IABS(IPTE(NIPTE))-NQPOI
            X1=    RAM(I3+1)
            X2=    RAM(I3+2)
            X3=    RAM(I3+3)
            IF (IPTE(NIPTE).LT.0) THEN
              ICBO1=IRAM(I3+8)
              ISBO1=IRAM(I3+7)
              ISRFO=-IRAM(I3+6)
              ISBO2=IRAM(I3+5)
              ICBO2=IRAM(I3+4)
            ELSE
              ICBO1=IRAM(I3+4)
              ISBO1=IRAM(I3+5)
              ISRFO=IRAM(I3+6)
              ISBO2=IRAM(I3+7)
              ICBO2=IRAM(I3+8)
            ENDIF
            NIPTE=NIPTE-1
C           Computing the point of intersection of the edge
C           with the gridface:
            CALL RKGS(PRMT,XX,DTMP,3,I7,FCTMS,OUTMS,AUX)
            IF (PRMT(5).NE.1.) THEN
C             MODSRF-08
              CALL ERROR('MODSRF-08: PRMT(5) not equal to 1.')
C             This error should not appear. Contact the author.
            ENDIF
            IF (.NOT.MSLPIF(XX,IFACE)) THEN
C             MODSRF-09
              CALL ERROR('MODSRF-09: Point outside the gridface.')
C             This error should not appear. Contact the author.
            ENDIF
C           Storing the first point on the edge:
            IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
            RAM(NPOI+1)=X1
            RAM(NPOI+2)=X2
            RAM(NPOI+3)=X3
            IRAM(NPOI+4)=ICBO1
            IRAM(NPOI+5)=ISBO1
            IRAM(NPOI+6)=ISRFO
            IRAM(NPOI+7)=ISBO2
            IRAM(NPOI+8)=ICBO2
            NPOI=NPOI+NQPOI
            J1=NPOI
C           Recording the connection:
            IF (NCON+3.GT.MRAM) CALL MSEROR(1)
            IF (ISRFO.LT.0) THEN
              IRAM(NCON+1)=IABS(IPTE(NIPTE+1))
              IRAM(NCON+2)=NPOI
              IRAM(NCON+3)=ISRFO
            ELSE
              IRAM(NCON+1)=NPOI
              IRAM(NCON+2)=IABS(IPTE(NIPTE+1))
              IRAM(NCON+3)=-ISRFO
            ENDIF
            NCON=NCON+3
C
            IF (ICBO1.NE.ICBN1) THEN
C             Storing the second point on the edge:
              IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
              RAM(NPOI+1)=X1
              RAM(NPOI+2)=X2
              RAM(NPOI+3)=X3
              IRAM(NPOI+4)=ICBN1
              IRAM(NPOI+5)=ISBN1
              IRAM(NPOI+6)=-ISRFN1
              IRAM(NPOI+7)=ISBO1
              IRAM(NPOI+8)=ICBO1
              NPOI=NPOI+NQPOI
              IF (NCUB.EQ.0) THEN
C               'Connection' between first and second point on the edge:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                J2=NPOI-NQPOI
                IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
                J3=NPOI
                IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
                IRAM(NCON+1)=J2
                IRAM(NCON+2)=J3
                IRAM(NCON+3)=0
                NCON=NCON+3
              ENDIF
              DO 40, I3=1,NJPT
                IF ((ICBN1.EQ.JPT(2,I3)).AND.
     *              (ISBN1.EQ.JPT(3,I3)).AND.
     *              (ISBO1.EQ.JPT(5,I3)).AND.
     *              (ICBO1.EQ.JPT(6,I3)).AND.
     *              (IABS(ISRFN1).EQ.IABS(JPT(4,I3)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (JPT(4,I3).LT.0) THEN
                    IRAM(NCON+1)=IABS(JPT(1,I3))
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=JPT(4,I3)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=IABS(JPT(1,I3))
                    IRAM(NCON+3)=-JPT(4,I3)
                  ENDIF
                  NCON=NCON+3
                  GOTO 405
                ENDIF
  40          CONTINUE
              DO 401, I3=1,NIPTE
                I4=IABS(IPTE(I3))-NQPOI
                IF ((ICBN1.EQ.IRAM(I4+8)).AND.
     *              (ISBN1.EQ.IRAM(I4+7)).AND.
     *              (ISBO1.EQ.IRAM(I4+5)).AND.
     *              (ICBO1.EQ.IRAM(I4+4)).AND.
     *              (IABS(ISRFN1).EQ.IABS(IRAM(I4+6)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (IRAM(I4+6).LT.0) THEN
                    IRAM(NCON+1)=I4+NQPOI
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=IRAM(I4+6)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=I4+NQPOI
                    IRAM(NCON+3)=-IRAM(I4+6)
                  ENDIF
                  NCON=NCON+3
                  GOTO 405
                ENDIF
  401         CONTINUE
C             More than one edge, current point will become a starting
C             point of the search for new edge:
              NIPTE=NIPTE+1
              IF (NIPTE.GT.MIPTE) CALL MSEROR(3)
              IPTE(NIPTE)=-NPOI
  405         CONTINUE
            ENDIF
C
            IF (ICBN2.NE.ICBN1) THEN
C             Storing the third point on the edge:
              IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
              RAM(NPOI+1)=X1
              RAM(NPOI+2)=X2
              RAM(NPOI+3)=X3
              IRAM(NPOI+4)=ICBN2
              IRAM(NPOI+5)=ISBN2
              IRAM(NPOI+6)=-ISRFO
              IRAM(NPOI+7)=ISBN1
              IRAM(NPOI+8)=ICBN1
              NPOI=NPOI+NQPOI
              IF (NCUB.EQ.0) THEN
C               'Connection' between second and third point on the edge:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                J2=NPOI-NQPOI
                IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
                J3=NPOI
                IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
                IRAM(NCON+1)=J2
                IRAM(NCON+2)=J3
                IRAM(NCON+3)=0
                NCON=NCON+3
              ENDIF
              DO 42, I3=1,NJPT
                IF ((ICBN2.EQ.JPT(2,I3)).AND.
     *              (ISBN2.EQ.JPT(3,I3)).AND.
     *              (ISBN1.EQ.JPT(5,I3)).AND.
     *              (ICBN1.EQ.JPT(6,I3)).AND.
     *              (IABS(ISRFO).EQ.IABS(JPT(4,I3)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (JPT(4,I3).LT.0) THEN
                    IRAM(NCON+1)=IABS(JPT(1,I3))
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=JPT(4,I3)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=IABS(JPT(1,I3))
                    IRAM(NCON+3)=-JPT(4,I3)
                  ENDIF
                  NCON=NCON+3
                  GOTO 425
                ENDIF
  42          CONTINUE
              DO 421, I3=1,NIPTE
                I4=IABS(IPTE(I3))-NQPOI
                IF ((ICBN2.EQ.IRAM(I4+8)).AND.
     *              (ISBN2.EQ.IRAM(I4+7)).AND.
     *              (ISBN1.EQ.IRAM(I4+5)).AND.
     *              (ICBN1.EQ.IRAM(I4+4)).AND.
     *              (IABS(ISRFO).EQ.IABS(IRAM(I4+6)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (IRAM(I4+6).LT.0) THEN
                    IRAM(NCON+1)=I4+NQPOI
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=IRAM(I4+6)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=I4+NQPOI
                    IRAM(NCON+3)=-IRAM(I4+6)
                  ENDIF
                  NCON=NCON+3
                  GOTO 425
                ENDIF
  421         CONTINUE
C             More than one edge, current point will become a starting
C             point of the search for new edge:
              NIPTE=NIPTE+1
              IF (NIPTE.GT.MIPTE) CALL MSEROR(3)
              IPTE(NIPTE)=-NPOI
  425         CONTINUE
            ENDIF
C
            IF (ICBO2.NE.ICBN2) THEN
C             Storing the fourth point on the edge:
              IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
              RAM(NPOI+1)=X1
              RAM(NPOI+2)=X2
              RAM(NPOI+3)=X3
              IRAM(NPOI+4)=ICBO2
              IRAM(NPOI+5)=ISBO2
              IRAM(NPOI+6)=ISRFN2
              IRAM(NPOI+7)=ISBN2
              IRAM(NPOI+8)=ICBN2
              NPOI=NPOI+NQPOI
              IF (NCUB.EQ.0) THEN
C               'Connection' between third and fourth point on the edge:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                J2=NPOI-NQPOI
                IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
                J3=NPOI
                IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
                IRAM(NCON+1)=J2
                IRAM(NCON+2)=J3
                IRAM(NCON+3)=0
                NCON=NCON+3
              ENDIF
              DO 41, I3=1,NJPT
                IF ((ICBO2.EQ.JPT(2,I3)).AND.
     *              (ISBO2.EQ.JPT(3,I3)).AND.
     *              (ISBN2.EQ.JPT(5,I3)).AND.
     *              (ICBN2.EQ.JPT(6,I3)).AND.
     *              (IABS(ISRFN2).EQ.IABS(JPT(4,I3)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (JPT(4,I3).LT.0) THEN
                    IRAM(NCON+1)=IABS(JPT(1,I3))
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=JPT(4,I3)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=IABS(JPT(1,I3))
                    IRAM(NCON+3)=-JPT(4,I3)
                  ENDIF
                  NCON=NCON+3
                  GOTO 415
                ENDIF
  41          CONTINUE
              DO 411, I3=1,NIPTE
                I4=IABS(IPTE(I3))-NQPOI
                IF ((ICBO2.EQ.IRAM(I4+8)).AND.
     *              (ISBO2.EQ.IRAM(I4+7)).AND.
     *              (ISBN2.EQ.IRAM(I4+5)).AND.
     *              (ICBN2.EQ.IRAM(I4+4)).AND.
     *              (IABS(ISRFN2).EQ.IABS(IRAM(I4+6)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (IRAM(I4+6).LT.0) THEN
                    IRAM(NCON+1)=I4+NQPOI
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=IRAM(I4+6)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=I4+NQPOI
                    IRAM(NCON+3)=-IRAM(I4+6)
                  ENDIF
                  NCON=NCON+3
                  GOTO 415
                ENDIF
  411         CONTINUE
C             More than one edge, current point will become a starting
C             point of the search for new edge:
              NIPTE=NIPTE+1
              IF (NIPTE.GT.MIPTE) CALL MSEROR(3)
              IPTE(NIPTE)=-NPOI
  415         CONTINUE
            ENDIF
C
            IF (NCUB.EQ.0) THEN
C             'Connection' between fourth and first point on the edge:
              IF (NCON+3.GT.MRAM) CALL MSEROR(1)
              J2=NPOI
              IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
              J3=J1
              IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
              IRAM(NCON+1)=J2
              IRAM(NCON+2)=J3
              IRAM(NCON+3)=0
              NCON=NCON+3
            ENDIF
          IF (NIPTE.GE.1)
C           Continuing with the search for next edge:
     *      GOTO 39
C
C         No more edges, the point is connected.
  45      CONTINUE
C         Continuing with next point of the gridface:
  46    CONTINUE
        IRAM(OFAC+I1)=NCON
  49  CONTINUE
      NPOIE=NPOI
C
C
      IF (NCUB.NE.0) THEN
C       Loop along all gridcubes,
C       merging the connections between the points on structural
C       interfaces into polygons approximating the interfaces.
        WRITE(*,'(A)')
     *'+MODSRF: Making polygons approximating the interfaces.          '
        NPOL=MRAM+1
        DO 100, I1=1,NCUB
C         Indices of gridfaces of the gridcube:
          CALL MSGCUB(I1,IF1,IF2,IF3,IF4,IF5,IF6)
C
C         Forming array with connections for the gridcube:
          NJCN=0
          NNJCN(0)=NJCN
          DO 51, I2=IRAM(OFAC+IF1-1),IRAM(OFAC+IF1)-3,3
C           Connections of the gridface IF1:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+1)
            JCN(2,NJCN)=IRAM(I2+2)
            JCN(3,NJCN)=IRAM(I2+3)
  51      CONTINUE
          NNJCN(1)=NJCN
          DO 52, I2=IRAM(OFAC+IF2-1),IRAM(OFAC+IF2)-3,3
C           Connections of the gridface IF2:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+1)
            JCN(2,NJCN)=IRAM(I2+2)
            JCN(3,NJCN)=IRAM(I2+3)
  52      CONTINUE
          NNJCN(2)=NJCN
          DO 53, I2=IRAM(OFAC+IF3-1),IRAM(OFAC+IF3)-3,3
C           Connections of the gridface IF3:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+1)
            JCN(2,NJCN)=IRAM(I2+2)
            JCN(3,NJCN)=IRAM(I2+3)
  53      CONTINUE
          NNJCN(3)=NJCN
          DO 54, I2=IRAM(OFAC+IF4-1),IRAM(OFAC+IF4)-3,3
C           Connections of the gridface IF4:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+2)
            JCN(2,NJCN)=IRAM(I2+1)
            JCN(3,NJCN)=IRAM(I2+3)
  54      CONTINUE
          NNJCN(4)=NJCN
          DO 55, I2=IRAM(OFAC+IF5-1),IRAM(OFAC+IF5)-3,3
C           Connections of the gridface IF5:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+2)
            JCN(2,NJCN)=IRAM(I2+1)
            JCN(3,NJCN)=IRAM(I2+3)
  55      CONTINUE
          NNJCN(5)=NJCN
          DO 56, I2=IRAM(OFAC+IF6-1),IRAM(OFAC+IF6)-3,3
C           Connections of the gridface IF6:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+2)
            JCN(2,NJCN)=IRAM(I2+1)
            JCN(3,NJCN)=IRAM(I2+3)
  56      CONTINUE
          NNJCN(6)=NJCN
          DO 58, I2=1,NNJCN(6)
C           Connections between points on edges:
            IF (JCN(2,I2).GT.NPOIL) THEN
              ICBO1=IRAM(JCN(2,I2)-NQPOI+4)
              ICBO2=IRAM(JCN(2,I2)-NQPOI+8)
              LCN=.FALSE.
              DO 57, I3=1,NNJCN(6)
                IF (I3.NE.I2) THEN
                  IF ((JCN(1,I3).GT.NPOIL).AND.
     *                (IABS(JCN(3,I2)).EQ.IABS(JCN(3,I3)))) THEN
                    ICBN1=IRAM(JCN(1,I3)-NQPOI+4)
                    ICBN2=IRAM(JCN(1,I3)-NQPOI+8)
                    IF (((ICBO1.EQ.ICBN1).AND.(ICBO2.EQ.ICBN2)).OR.
     *                  ((ICBO2.EQ.ICBN1).AND.(ICBO1.EQ.ICBN2))) THEN
                      NJCN=NJCN+1
                      IF (NJCN.GT.MJCN) CALL MSEROR(4)
                      JCN(1,NJCN)=JCN(2,I2)
                      JCN(2,NJCN)=JCN(1,I3)
                      JCN(3,NJCN)=JCN(3,I3)
                      LCN=.TRUE.
                    ENDIF
                  ENDIF
                ENDIF
  57          CONTINUE
              IF (.NOT.LCN) THEN
C               Connection for edge not found, this connection is to be
C               canceled by means of creation of linear 'triangle':
                NJCN=NJCN+1
                IF (NJCN.GT.MJCN) CALL MSEROR(4)
                JCN(1,NJCN)=JCN(2,I2)
                JCN(2,NJCN)=JCN(1,I2)
                JCN(3,NJCN)=JCN(3,I2)
              ENDIF
            ENDIF
  58      CONTINUE
          NNJCN(7)=NJCN
C
C         Forming array with starting points
C         (from point, point, to point)
          NJSP=0
          DO 60, I2=1,NJCN
            DO 59, I3=1,NJCN
              IF ((JCN(2,I2).EQ.JCN(1,I3)).AND.
     *            (JCN(2,I3).NE.JCN(1,I2))) THEN
                NJSP=NJSP+1
                IF (NJSP.GT.MJSP) CALL MSEROR(5)
                JSP(1,NJSP)=JCN(1,I2)
                JSP(2,NJSP)=JCN(2,I2)
                JSP(3,NJSP)=JCN(2,I3)
              ENDIF
  59        CONTINUE
  60      CONTINUE
C
C         Looking for polygons on the faces of the cube.
          NJPOL=0
C         Loop for faces of the cube:
          DO 70, I2=1,6
C           Initiating indices of connections:
            DO 61, I3=1,MJCN
              IND(I3)=0
  61        CONTINUE
            NIND=NNJCN(I2)-NNJCN(I2-1)
C           Loop along connections:
            I3=1
  62        CONTINUE
              IF (NJPOL+1.GT.MJPOL) CALL MSEROR(6)
C             Search for first unused connection:
              IF (I3.LE.NIND-2) THEN
                IF (IND(I3).GE.0) THEN
C                 Unused connection, recording first two points
C                 of the polygon:
                  NLJPOL=2
                  JPOL(1,NJPOL+1)=JCN(1,NNJCN(I2-1)+I3)
                  JPOL(2,NJPOL+1)=JCN(2,NNJCN(I2-1)+I3)
                  IND(I3)=-1
                ELSE
                  I3=I3+1
                  GOTO 62
                ENDIF
C               Search for next point of polygon:
  63            CONTINUE
                DO 69, I4=I3+1,NIND
                  IF ( JCN(1,NNJCN(I2-1)+I4)       .EQ.
     *                 JPOL(NLJPOL,NJPOL+1) ) THEN
C                   Recording next point of polygon:
                    NLJPOL=NLJPOL+1
                    IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7)
                    JPOL(NLJPOL,NJPOL+1)=JCN(2,NNJCN(I2-1)+I4)
                    IND(I4)=-1
                    DO 68, I5=NLJPOL-3,1,-1
C                     Examining whether the polygon is closed:
                      IF (JPOL(I5,NJPOL+1).EQ.JPOL(NLJPOL,NJPOL+1)) THEN
C                       Polygon is closed.
C                       Recording the polygon:
                        NJPOL=NJPOL+1
                        NLJPOL=NLJPOL-I5
                        JPOL(0,NJPOL)=NLJPOL
                        DO 65, I6=1,NLJPOL
                          JPOL(I6,NJPOL)=JPOL(I5+I6,NJPOL)
  65                    CONTINUE
C                       Removing used starting points:
                        DO 67, I6=1,NLJPOL
                          IF (I6.EQ.1) THEN
                            J1=JPOL(NLJPOL,NJPOL)
                          ELSE
                            J1=JPOL(I6-1,NJPOL)
                          ENDIF
                          J2=JPOL(I6,NJPOL)
                          IF (I6.EQ.NJPOL) THEN
                            J3=JPOL(1,NJPOL)
                          ELSE
                            J3=JPOL(I6+1,NJPOL)
                          ENDIF
                          DO 66, I7=1,NJSP
                            IF ((J1.EQ.JSP(1,I7)).AND.(J2.EQ.JSP(2,I7))
     *                          .AND.(J3.EQ.JSP(3,I7))) THEN
                              JSP(1,I7)=0
                              JSP(2,I7)=0
                              JSP(3,I7)=0
                            ENDIF
  66                      CONTINUE
  67                    CONTINUE
C                       Continuing with next polygon:
                        GOTO 62
                      ENDIF
  68                CONTINUE
                    GOTO 63
                  ENDIF
  69            CONTINUE
C               Polygon is opened, searching for other polygon:
                GOTO 62
              ENDIF
C             No more unused connections,
C             no more polygons on this face of the cube.
C           End of the loop along connections.
  70      CONTINUE
C
C         If there are two or more polygons on the faces of the cube,
C         the polygons are to be omitted:
          IF (NJPOL.GE.2) THEN
C           MODSRF-10
            CALL ERROR('MODSRF-10: Two or more polygons on gridface.')
C           This situation cannot be handled by current version of
C           MODSRF. Try to reduce the gridstep (D1,D2,D3).
          ENDIF
C
C         If there is one polygon on the faces of the cube,
C         the polygon is to be recorded as a polygon approximating
C         a structural interface:
          IF (NJPOL.EQ.1) THEN
C           Searching, whether the polygon is to be recorded:
C           Loop for already recorded polygons:
            I2=MRAM+1
  73        CONTINUE
              IF (I2.GT.NPOL) THEN
                J1=IRAM(I2-1)
                I2=I2-J1-1
                IF (J1.EQ.JPOL(0,1)) THEN
                  DO 74, I3=1,J1
                    IF (IRAM(I2+I3).NE.JPOL(I3,1)) GOTO 73
  74              CONTINUE
C                 The polygon is already recorded:
                  GOTO 77
                ELSE
                  GOTO 73
                ENDIF
              ENDIF
C           End of the loop for already recorded polygons.
C           Recording the polygon:
            IF (NPOL-JPOL(0,1)-1.LE.NCON) CALL MSEROR(1)
            NPOL=NPOL-1
            IRAM(NPOL)=JPOL(0,1)
            DO 75, I2=1,JPOL(0,1)
              NPOL=NPOL-1
              IRAM(NPOL)=JPOL(I2,1)
  75        CONTINUE
  77        CONTINUE
          ENDIF
C
C         Recording polygons approximating structural interfaces:
          NJPOL=0
  80      CONTINUE
C         Loop for adding new polygons
C           Initializing a polygon:
            DO 82, I2=1,NJSP
C             Prefering points on edges:
              IF (JSP(1,I2).GT.NPOIL) THEN
                NJPOL=NJPOL+1
                IF (NJPOL.GT.MJPOL) CALL MSEROR(6)
                NLJPOL=3
                JPOL(0,NJPOL)=0
                JPOL(1,NJPOL)=JSP(1,I2)
                JPOL(2,NJPOL)=JSP(2,I2)
                JPOL(3,NJPOL)=JSP(3,I2)
                JSP(1,I2)=0
                JSP(2,I2)=0
                JSP(3,I2)=0
                GOTO  85
              ENDIF
  82        CONTINUE
            DO 84, I2=1,NJSP
C             No points on edges, taking remaining points:
              IF (JSP(1,I2).NE.0) THEN
                NJPOL=NJPOL+1
                IF (NJPOL.GT.MJPOL) CALL MSEROR(6)
                NLJPOL=3
                JPOL(0,NJPOL)=0
                JPOL(1,NJPOL)=JSP(1,I2)
                JPOL(2,NJPOL)=JSP(2,I2)
                JPOL(3,NJPOL)=JSP(3,I2)
                JSP(1,I2)=0
                JSP(2,I2)=0
                JSP(3,I2)=0
                GOTO  85
              ENDIF
  84        CONTINUE
C           No other point to initialize a polygon:
            GOTO 95
C
  85        CONTINUE
C           Loop for adding a points to the polygon:
              DO 93, I2=1,NJSP
                IF ((JSP(1,I2).EQ.JPOL(NLJPOL-1,NJPOL)).AND.
     *              (JSP(2,I2).EQ.JPOL(NLJPOL  ,NJPOL))) THEN
                  DO 87, I3=2,NLJPOL-2
                    IF (JSP(3,I2).EQ.JPOL(I3,NJPOL)) THEN
C                     This point would cause overlooping of the polygon:
                      JSP(1,I2)=0
                      JSP(2,I2)=0
                      JSP(3,I2)=0
                      DO 86, I4=1,NJSP
                        IF ((JSP(1,I4).EQ.JPOL(NLJPOL,NJPOL)).AND.
     *                      (JSP(2,I4).EQ.JPOL(I3    ,NJPOL)).AND.
     *                      (JSP(3,I4).EQ.JPOL(I3+1  ,NJPOL))) THEN
                          JSP(1,I4)=0
                          JSP(2,I4)=0
                          JSP(3,I4)=0
                          GOTO 85
                        ENDIF
  86                  CONTINUE
C                     
C                     MODSRF-11
                      CALL ERROR('MODSRF-11: Disorder in array JSP.')
C                     This error should not appear. Contact the author.
                    ENDIF
  87              CONTINUE
                  IF (JSP(3,I2).EQ.JPOL(1,NJPOL)) THEN
C                   This point would close the polygon:
                    DO 91, I3=I2+1,NJSP
C                     Looking for other possible continuation, which
C                     would make polygon longer:
                      IF ((JSP(1,I3).EQ.JPOL(NLJPOL-1,NJPOL)).AND.
     *                    (JSP(2,I3).EQ.JPOL(NLJPOL  ,NJPOL)).AND.
     *                    (JSP(3,I3).NE.JPOL(1,NJPOL))) THEN
C                       Other possible continuation of the polygon:
                        DO 90, I4=2,NLJPOL-2
                          IF (JSP(3,I3).EQ.JPOL(I4,NJPOL)) THEN
C                           This point would cause overlooping of the
C                           polygon:
                            JSP(1,I3)=0
                            JSP(2,I3)=0
                            JSP(3,I3)=0
                            DO 89, I5=1,NJSP
                              IF (   (JSP(1,I5).EQ.JPOL(NLJPOL,NJPOL))
     *                          .AND.(JSP(2,I5).EQ.JPOL(I4    ,NJPOL))
     *                          .AND.(JSP(3,I5).EQ.JPOL(I4+1  ,NJPOL)))
     *                          THEN
                                JSP(1,I5)=0
                                JSP(2,I5)=0
                                JSP(3,I5)=0
                                GOTO 85
                              ENDIF
  89                        CONTINUE
C                           
C                           MODSRF-12
                            CALL ERROR
     *                      ('MODSRF-12: Disorder in array JSP.')
C                           This error should not appear.
C                           Contact the author.
                          ENDIF
  90                    CONTINUE
C                       Continuation I2 is to be omitted because
C                       continuation I3 makes the polygon longer:
                        JSP(1,I2)=0
                        JSP(2,I2)=0
                        JSP(3,I2)=0
                        IF (JPOL(0,NJPOL).EQ.0) THEN
C                         Removing the second point which closes the
C                         polygon from JSP:
                          DO 88, I4=1,NJSP
                            IF ((JSP(1,I4).EQ.JPOL(NLJPOL,NJPOL)).AND.
     *                          (JSP(2,I4).EQ.JPOL(1,NJPOL)).AND.
     *                          (JSP(3,I4).EQ.JPOL(2,NJPOL))) THEN
                              JSP(1,I4)=0
                              JSP(2,I4)=0
                              JSP(3,I4)=0
                              JPOL(0,NJPOL)=0
                              GOTO 85
                            ENDIF
  88                      CONTINUE
C                         
C                         MODSRF-13
                          CALL ERROR
     *                    ('MODSRF-13: Disorder in array JSP.')
C                         This error should not appear.
C                         Contact the author.
                        ELSE
                          JPOL(0,NJPOL)=0
                          GOTO 85
                        ENDIF
                      ENDIF
  91                CONTINUE
C                   No other continuation, polygon is to be closed:
                    JSP(1,I2)=0
                    JSP(2,I2)=0
                    JSP(3,I2)=0
                    IF (JPOL(0,NJPOL).EQ.0) THEN
C                     First point which closes the polygon.
                      JPOL(0,NJPOL)=NLJPOL
                      GOTO 85
                    ELSE
C                     Second point which closes the polygon.
                      GOTO 80
                    ENDIF
                  ELSE
C                   Adding a point to the polygon:
                    NLJPOL=NLJPOL+1
                    IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7)
                    JPOL(NLJPOL,NJPOL)=JSP(3,I2)
                    JSP(1,I2)=0
                    JSP(2,I2)=0
                    JSP(3,I2)=0
                    GOTO 85
                  ENDIF
                ELSEIF ((JSP(1,I2).EQ.JPOL(NLJPOL,NJPOL)).AND.
     *                  (JSP(2,I2).EQ.JPOL(1,NJPOL)).AND.
     *                  (JSP(3,I2).EQ.JPOL(2,NJPOL))) THEN
C                 The polygon is closed.
                  JSP(1,I2)=0
                  JSP(2,I2)=0
                  JSP(3,I2)=0
                  IF (JPOL(0,NJPOL).EQ.0) THEN
C                   First point which closes the polygon.
                    JPOL(0,NJPOL)=NLJPOL
                    GOTO 85
                  ELSE
C                   Second point which closes the polygon.
                    GOTO 80
                  ENDIF
                ENDIF
  93          CONTINUE
C             MODSRF-14
              CALL ERROR ('MODSRF-14: Polygon opened.')
C             This error should not appear. Contact the author.
C           End of the loop for adding new points to the polygon.
C         End of the loop for adding new polygons.
C
  95      CONTINUE
C         Writing polygons to IRAM:
          DO 97, I2=1,NJPOL
            IF (NPOL-JPOL(0,I2)-1.LE.NCON) CALL MSEROR(1)
            NPOL=NPOL-1
            IRAM(NPOL)=JPOL(0,I2)
            DO 96, I3=1,JPOL(0,I2)
              NPOL=NPOL-1
              IRAM(NPOL)=JPOL(I3,I2)
  96        CONTINUE
  97      CONTINUE
          NJPOL=0
C
C       End of the loop over gridcubes.
 100    CONTINUE
C
C
      ELSE
C       Loop over faces of the 2-D grid, making polygons along the 2-D
C       slice through the model.
        WRITE(*,'(A)')
     *'+MODSRF: Making polygons along the 2-D slice.                   '
        NPOL=MRAM+1
C
C       Rewriting points on interfaces as points on + side of interface:
        ISHIFT=NPOIE-OPOI
        IF (NPOI+ISHIFT.GT.OCON) CALL MSEROR(1)
        DO 102, I1=OPOI,NPOIE-NQPOI,NQPOI
          RAM(NPOI+1)=RAM(I1+1)
          RAM(NPOI+2)=RAM(I1+2)
          RAM(NPOI+3)=RAM(I1+3)
          IF (IRAM(I1+6).LT.0) THEN
            IRAM(NPOI+4)=IRAM(I1+8)
            IRAM(NPOI+5)=IRAM(I1+7)
            IRAM(NPOI+6)=0
            IRAM(NPOI+7)=IRAM(I1+7)
            IRAM(NPOI+8)=IRAM(I1+8)
          ELSE
            IRAM(NPOI+4)=IRAM(I1+4)
            IRAM(NPOI+5)=IRAM(I1+5)
            IRAM(NPOI+6)=0
            IRAM(NPOI+7)=IRAM(I1+5)
            IRAM(NPOI+8)=IRAM(I1+4)
          ENDIF
          NPOI=NPOI+NQPOI
 102    CONTINUE
        NPOILP=NPOIL+ISHIFT
        NPOIEP=NPOIE+ISHIFT
C
C       Rewriting points on interfaces as points on - side of interface:
        DO 104, I1=OPOI,NPOIE-NQPOI,NQPOI
          IF (IRAM(I1+6).LT.0) THEN
            IRAM(I1+7)=IRAM(I1+5)
            IRAM(I1+8)=IRAM(I1+4)
          ELSE
            IRAM(I1+4)=IRAM(I1+8)
            IRAM(I1+5)=IRAM(I1+7)
          ENDIF
          IRAM(NPOI+6)=0
 104    CONTINUE
C
C       Recording gridpoints to the array of points:
        IF (NPOI+NGPS*NQPOI.GT.OCON) CALL MSEROR(1)
        DO 106, I1=1,NGPS
          CALL MSCP(I1,X1,X2,X3)
          RAM(NPOI+1)=X1
          RAM(NPOI+2)=X2
          RAM(NPOI+3)=X3
          IRAM(NPOI+4)=IRAM(2*(I1-1)+2)
          IRAM(NPOI+5)=IRAM(2*(I1-1)+1)
          IRAM(NPOI+6)=0
          IRAM(NPOI+7)=IRAM(2*(I1-1)+1)
          IRAM(NPOI+8)=IRAM(2*(I1-1)+2)
          NPOI=NPOI+NQPOI
 106    CONTINUE
C       Loop over gridfaces:
        DO 200, I1=1,NFAC
C         Forming array with connections for the gridface:
          NJCN=0
C         Connections inside the gridface:
          DO 110, I2=IRAM(OFAC+I1-1),IRAM(OFAC+I1)-3,3
C           Connections of the gridface I1:
            IF (IRAM(I2+3).NE.0) THEN
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN) CALL MSEROR(4)
              JCN(1,NJCN)=IRAM(I2+1)
              JCN(2,NJCN)=IRAM(I2+2)
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN) CALL MSEROR(4)
              JCN(1,NJCN)=IRAM(I2+2)+ISHIFT
              JCN(2,NJCN)=IRAM(I2+1)+ISHIFT
              JCN(3,NJCN)=1
            ELSE
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN) CALL MSEROR(4)
              IF (IRAM(I2+1).LT.0) THEN
                JCN(1,NJCN)=-IRAM(I2+1)
              ELSE
                JCN(1,NJCN)=IRAM(I2+1)+ISHIFT
              ENDIF
              IF (IRAM(I2+2).LT.0) THEN
                JCN(2,NJCN)=-IRAM(I2+2)
              ELSE
                JCN(2,NJCN)=IRAM(I2+2)+ISHIFT
              ENDIF
              JCN(3,NJCN)=1
            ENDIF
 110      CONTINUE
C         Connections along the gridface:
C         Indices of gridlegs of the gridface:
          CALL MSGFAC(I1,IL1,IL2,IL3,IL4)
C
          CALL MSGLEG(IL1,IP1,IP2)
          NJCN=NJCN+1
          IF (NJCN.GT.MJCN)  CALL MSEROR(4)
          JCN(1,NJCN)=NPOIEP+IP1*NQPOI
          DO 112, I2=IRAM(OLEG+IL1-1)+NQPOI,IRAM(OLEG+IL1),NQPOI
            ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4)
            ICBNEW=IRAM(I2         -NQPOI+4)
            IF (ICBOLD.NE.ICBNEW) THEN
              JCN(2,NJCN)=I2+ISHIFT
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2
            ELSE
              JCN(2,NJCN)=I2
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2+ISHIFT
            ENDIF
 112      CONTINUE
          JCN(2,NJCN)=NPOIEP+IP2*NQPOI
          JCN(3,NJCN)=1
C
          CALL MSGLEG(IL2,IP1,IP2)
          NJCN=NJCN+1
          IF (NJCN.GT.MJCN)  CALL MSEROR(4)
          JCN(1,NJCN)=NPOIEP+IP1*NQPOI
          DO 114, I2=IRAM(OLEG+IL2-1)+NQPOI,IRAM(OLEG+IL2),NQPOI
            ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4)
            ICBNEW=IRAM(I2         -NQPOI+4)
            IF (ICBOLD.NE.ICBNEW) THEN
              JCN(2,NJCN)=I2+ISHIFT
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2
            ELSE
              JCN(2,NJCN)=I2
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2+ISHIFT
            ENDIF
 114      CONTINUE
          JCN(2,NJCN)=NPOIEP+IP2*NQPOI
          JCN(3,NJCN)=1
C
          CALL MSGLEG(IL3,IP2,IP1)
          NJCN=NJCN+1
          IF (NJCN.GT.MJCN)  CALL MSEROR(4)
          JCN(1,NJCN)=NPOIEP+IP1*NQPOI
          DO 116, I2=IRAM(OLEG+IL3),IRAM(OLEG+IL3-1)+NQPOI,-NQPOI
            ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4)
            ICBNEW=IRAM(I2         -NQPOI+4)
            IF (ICBOLD.NE.ICBNEW) THEN
              JCN(2,NJCN)=I2+ISHIFT
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2
            ELSE
              JCN(2,NJCN)=I2
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2+ISHIFT
            ENDIF
 116      CONTINUE
          JCN(2,NJCN)=NPOIEP+IP2*NQPOI
          JCN(3,NJCN)=1
C
          CALL MSGLEG(IL4,IP2,IP1)
          NJCN=NJCN+1
          IF (NJCN.GT.MJCN)  CALL MSEROR(4)
          JCN(1,NJCN)=NPOIEP+IP1*NQPOI
          DO 118, I2=IRAM(OLEG+IL4),IRAM(OLEG+IL4-1)+NQPOI,-NQPOI
            ICBOLD=IRAM(JCN(1,NJCN)-NQPOI+4)
            ICBNEW=IRAM(I2         -NQPOI+4)
            IF (ICBOLD.NE.ICBNEW) THEN
              JCN(2,NJCN)=I2+ISHIFT
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2
            ELSE
              JCN(2,NJCN)=I2
              JCN(3,NJCN)=1
              NJCN=NJCN+1
              IF (NJCN.GT.MJCN)  CALL MSEROR(4)
              JCN(1,NJCN)=I2+ISHIFT
            ENDIF
 118      CONTINUE
          JCN(2,NJCN)=NPOIEP+IP2*NQPOI
          JCN(3,NJCN)=1
C         Array with connections is done.
C
C         Forming polygons along the gridface:
 124      CONTINUE
          DO 130, I2=1,NJCN
            IF (JCN(3,I2).NE.0) THEN
C             Initiating polygon:
              NLJPOL=2
              JPOL(1,1)=JCN(1,I2)
              JPOL(2,1)=JCN(2,I2)
              JCN(3,I2)=0
 126          CONTINUE
              DO 128, I3=1,NJCN
                IF ((JCN(3,I3).NE.0).AND.
     *              (JCN(1,I3).EQ.JPOL(NLJPOL,1))) THEN
                  NLJPOL=NLJPOL+1
                  IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7)
                  JPOL(NLJPOL,1)=JCN(2,I3)
                  JCN(3,I3)=0
                  IF (JPOL(NLJPOL,1).EQ.JPOL(1,1)) THEN
C                   Polygon is closed - recording the polygon to IRAM:
                    NLJPOL=NLJPOL-1
                    IF (NPOL-NLJPOL-1.LE.NCON) CALL MSEROR(1)
                    NPOL=NPOL-1
                    IRAM(NPOL)=NLJPOL
                    DO 127, I4=1,NLJPOL
                      NPOL=NPOL-1
                      IRAM(NPOL)=JPOL(I4,1)
 127                CONTINUE
                    GOTO 124
                  ELSE
                    GOTO 126
                  ENDIF
                ENDIF
 128          CONTINUE
C             MODSRF-15
              CALL ERROR ('MODSRF-15: Polygon not closed.')
C             This error should not appear. Contact the author.
            ENDIF
 130      CONTINUE
C       End of the loop along gridfaces of the 2-D slice.
 200    CONTINUE
C
C       Preparing the normal to the 2-D slice:
        F(2)=0.
        F(3)=0.
        F(4)=0.
        IF (N1.EQ.1) F(2)=1.
        IF (N2.EQ.1) F(3)=1.
        IF (N3.EQ.1) F(4)=1.
C
      ENDIF
C
C
C     Storing points on structural interfaces and polygons:
      WRITE(*,'(A)')
     *'+MODSRF: Writing.                                               '
C
C     Storing points:
      OPEN(LU,FILE=VRTX,FORM='FORMATTED')
      WRITE(LU,'(A)') '/'
      NPOINT=(NPOI-OPOI)/NQPOI
C     Length of the names of points:
      LEN1=LENGTH(TEXTP)
      LEN2=0
 202  CONTINUE
        IF (NPOINT.GE.10**LEN2) THEN
          LEN2=LEN2+1
          GOTO 202
        ENDIF
      LENG=LEN1+LEN2
C     Loop over points:
      DO 218, I1=1,NPOINT
C       Address in RAM just before the current point:
        J1=OPOI+(I1-1)*NQPOI
C       Name of the point:
        DO 204, I2=0,LEN2-1
          TEXTP(LENG-I2:LENG-I2)=
     *    CHAR(ICHAR('0')+MOD(I1,10**(I2+1))/10**I2)
 204    CONTINUE
C       Coordinates of the point:
        Z(1)=RAM(J1+1)
        Z(2)=RAM(J1+2)
        Z(3)=RAM(J1+3)
C       Normal to the interface at the point:
        IF (NCUB.NE.0) CALL SRFC2(IRAM(J1+6),RAM(J1+1),F)
        Z(4)=F(2)
        Z(5)=F(3)
        Z(6)=F(4)
C       Other quantities:
        J2=6
 206    CONTINUE
          IF ((TEXTC(J2+1).NE.' ').AND.(J2+1.LE.69)) THEN
            J2=J2+1
            IF (IRAM(J1+6).LT.0) THEN
              ICBPOS=IRAM(J1+8)
              ICBNEG=IRAM(J1+4)
            ELSE
              ICBPOS=IRAM(J1+4)
              ICBNEG=IRAM(J1+8)
            ENDIF
            IF     (TEXTC(J2).EQ.'isrf') THEN
              Z(J2)=FLOAT(IABS(IRAM(J1+6)))**POWER(J2)
            ELSEIF (TEXTC(J2).EQ.'+icb') THEN
              Z(J2)=FLOAT(ICBPOS)**POWER(J2)
            ELSEIF (TEXTC(J2).EQ.'-icb') THEN
              Z(J2)=FLOAT(ICBNEG)**POWER(J2)
            ELSEIF (TEXTC(J2).EQ.'+isb') THEN
              IF (IRAM(J1+6).LT.0) THEN
                Z(J2)=FLOAT(IRAM(J1+7))**POWER(J2)
              ELSE
                Z(J2)=FLOAT(IRAM(J1+5))**POWER(J2)
              ENDIF
            ELSEIF (TEXTC(J2).EQ.'-isb') THEN
              IF (IRAM(J1+6).LT.0) THEN
                Z(J2)=FLOAT(IRAM(J1+5))**POWER(J2)
              ELSE
                Z(J2)=FLOAT(IRAM(J1+7))**POWER(J2)
              ENDIF
            ELSEIF ((TEXTC(J2)(2:2).EQ.'v')
     *          .OR.(TEXTC(J2)(2:2).EQ.'d')
     *          .OR.(TEXTC(J2)(2:3).EQ.'qp')
     *          .OR.(TEXTC(J2)(2:3).EQ.'qs')) THEN
              VP(1)=0.
              VS(1)=0.
              RHO=0.
              QP=0.
              QS=0.
              IF (TEXTC(J2)(1:1).EQ.'+') THEN
                IF (ICBPOS.NE.0) CALL PARM2(ICBPOS,Z(1),VP,VS,RHO,QP,QS)
              ELSE
                IF (ICBNEG.NE.0) CALL PARM2(ICBNEG,Z(1),VP,VS,RHO,QP,QS)
              ENDIF
              IF     (TEXTC(J2)(2:3).EQ.'vp') THEN
                Z(J2)=VP(1)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:3).EQ.'vs') THEN
                Z(J2)=VS(1)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:3).EQ.'de') THEN
                Z(J2)=RHO**POWER(J2)
              ELSEIF (TEXTC(J2)(2:3).EQ.'qp') THEN
                Z(J2)=QP**POWER(J2)
              ELSEIF (TEXTC(J2)(2:3).EQ.'qs') THEN
                Z(J2)=QS**POWER(J2)
              ENDIF
            ELSEIF ((TEXTC(J2)(2:2).EQ.'a')
     *          .OR.(TEXTC(J2)(2:2).EQ.'q')) THEN
              DO 208, I3=1,21
                A(1,I3)=0.
                Q(I3)=0.
 208          CONTINUE
              IF (TEXTC(J2)(1:1).EQ.'+') THEN
                IF (ICBPOS.NE.0) CALL PARM3(ICBPOS,Z(1),A,RHO,Q)
              ELSE
                IF (ICBNEG.NE.0) CALL PARM3(ICBNEG,Z(1),A,RHO,Q)
              ENDIF
              IF     (TEXTC(J2)(2:4).EQ.'a11') THEN
                Z(J2)=A(1, 1)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a12') THEN
                Z(J2)=A(1, 2)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a22') THEN
                Z(J2)=A(1, 3)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a13') THEN
                Z(J2)=A(1, 4)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a23') THEN
                Z(J2)=A(1, 5)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a33') THEN
                Z(J2)=A(1, 6)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a14') THEN
                Z(J2)=A(1, 7)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a24') THEN
                Z(J2)=A(1, 8)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a34') THEN
                Z(J2)=A(1, 9)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a44') THEN
                Z(J2)=A(1,10)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a15') THEN
                Z(J2)=A(1,11)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a25') THEN
                Z(J2)=A(1,12)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a35') THEN
                Z(J2)=A(1,13)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a45') THEN
                Z(J2)=A(1,14)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a55') THEN
                Z(J2)=A(1,15)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a16') THEN
                Z(J2)=A(1,16)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a26') THEN
                Z(J2)=A(1,17)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a36') THEN
                Z(J2)=A(1,18)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a46') THEN
                Z(J2)=A(1,19)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a56') THEN
                Z(J2)=A(1,20)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'a66') THEN
                Z(J2)=A(1,21)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q11') THEN
                Z(J2)=Q( 1)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q12') THEN
                Z(J2)=Q( 2)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q22') THEN
                Z(J2)=Q( 3)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q13') THEN
                Z(J2)=Q( 4)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q23') THEN
                Z(J2)=Q( 5)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q33') THEN
                Z(J2)=Q( 6)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q14') THEN
                Z(J2)=Q( 7)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q24') THEN
                Z(J2)=Q( 8)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q34') THEN
                Z(J2)=Q( 9)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q44') THEN
                Z(J2)=Q(10)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q15') THEN
                Z(J2)=Q(11)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q25') THEN
                Z(J2)=Q(12)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q35') THEN
                Z(J2)=Q(13)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q45') THEN
                Z(J2)=Q(14)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q55') THEN
                Z(J2)=Q(15)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q16') THEN
                Z(J2)=Q(16)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q26') THEN
                Z(J2)=Q(17)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q36') THEN
                Z(J2)=Q(18)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q46') THEN
                Z(J2)=Q(19)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q56') THEN
                Z(J2)=Q(20)**POWER(J2)
              ELSEIF (TEXTC(J2)(2:4).EQ.'q66') THEN
                Z(J2)=Q(21)**POWER(J2)
              ENDIF
            ENDIF
            GOTO 206
          ENDIF
C       End of the loop for other quantities.
C       Writing the quantities at the point:
C       Setting output format for the array:
        FORMAT='(3A,00(F00.0,1X),A)'
        FORMAT(6:6)=CHAR(ICHAR('0')+MOD(J2,10))
        FORMAT(5:5)=CHAR(ICHAR('0')+J2/10)
        OUTMIN=0.
        OUTMAX=0.
        DO 214, I2=1,J2
          IF (OUTMIN.GT.Z(I2)) OUTMIN=Z(I2)
          IF (OUTMAX.LT.Z(I2)) OUTMAX=Z(I2)
 214    CONTINUE
        CALL FORM1(OUTMIN,OUTMAX,FORMAT(8:15))
        FORMAT(14:17)= '1X),'
C       Output format is set.
        WRITE(LU,FORMAT) '''',TEXTP(1:LENG),''' ',(Z(I2),I2=1,J2),'/'
 218  CONTINUE
C     End of the loop over all points.
      WRITE(LU,'(A)') '/'
      CLOSE(LU)
C
C     Storing polygons:
      IF(PLGN.NE.' ') OPEN(LU1,FILE=PLGN,FORM='FORMATTED')
      IF(PLGNS.NE.' ')OPEN(LU2,FILE=PLGNS,FORM='FORMATTED')
      FORMA1='(00(I8,1X),A)'
      FORMA2='(00(3A,1X),A)'
      I1=MRAM
C     Loop over polygons:
 220  CONTINUE
      IF (I1.GT.NPOL) THEN
        J1=IRAM(I1)
        I3=MOD(J1,10)
        FORMA1(3:3)=CHAR(ICHAR('0')+I3)
        FORMA2(3:3)=CHAR(ICHAR('0')+I3)
        I3=J1/10
        FORMA1(2:2)=CHAR(ICHAR('0')+I3)
        FORMA2(2:2)=CHAR(ICHAR('0')+I3)
        J2=IRAM(I1-1)
        IF ((IRAM(J2-NQPOI+4).NE.0).OR.
     *      (IRAM(J2-NQPOI+8).NE.0)) THEN
          IF (PLGN.NE.' ') THEN
            WRITE(LU1,FORMA1)
     *           ((IRAM(I2)-OPOI)/NQPOI,I2=I1-1,I1-J1,-1),'/'
          ENDIF
          IF (PLGNS.NE.' ') THEN
            IF (J1.GT.MLJPOL) THEN
C             MODSRF-16
              CALL ERROR ('MODSRF-16: Disorder in polygons.')
C             This error should not appear. Contact the author.
            ENDIF
            DO 224, I2=I1-1,I1-J1,-1
C             Index of the point (from 1 to NPOINT):
              J2=(IRAM(I2)-OPOI)/NQPOI
C             Name of the point:
              I4=I1-I2
              DO 222, I3=0,LEN2-1
                TEXTS(I4)(1:LEN1)=TEXTP(1:LEN1)
                TEXTS(I4)(LENG-I3:LENG-I3)=
     *          CHAR(ICHAR('0')+MOD(J2,10**(I3+1))/10**I3)
 222          CONTINUE
 224        CONTINUE
            WRITE(LU2,FORMA2)
     *           ('''',TEXTS(I2)(1:LENG),'''',I2=1,J1),'/'
          ENDIF
        ENDIF
        I1=I1-(J1+1)
        GOTO 220
      ENDIF
C     End of the loop over polygons.
      IF (PLGN.NE.' ') THEN
        WRITE(LU1,'(A)') '/'
        CLOSE(LU1)
      ENDIF
      IF (PLGNS.NE.' ') THEN
        WRITE(LU2,'(A)') '/'
        CLOSE(LU2)
      ENDIF
C
      WRITE(*,'(A)')
     *'+MODSRF: Finished                                              '
      END
C
C=======================================================================
C
      SUBROUTINE FCTMS(X,Y,T)
C
C-----------------------------------------------------------------------
C
      REAL X,Y(*),T(*)
C
C Subroutine evaluating the right-hand sides of the interface tracing
C equations. I.E. subroutine computing vector T tangent to the
C interface.
C
C Input:
C     X...    Value of the independent variable.
C     Y...    Array containing two coordinates of a
C             point of the interface, determined by means of numerical
C             integration.
C Output:
C     Y...    Array containing two coordinates of the
C             interface, corrected by means of the linearization in
C             the direction of the gradient (perpendicular to the
C             interface).
C     T...    Array containing derivatives of the coordinates
C             Y with respect to X (vector tangent to the interface).
C-----------------------------------------------------------------------
C
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      REAL S1,S2,S3,F(10),AUX
C
C.......................................................................
C
      CALL SRFC2(ISRFO,Y,F)
      IF (ISRFO.LT.0) THEN
        S1=F(2)
        S2=F(3)
        S3=F(4)
      ELSE
        S1=-F(2)
        S2=-F(3)
        S3=-F(4)
      ENDIF
      IF (IFACE.LE.NFAC1) THEN
        F(2)=0.
        S1=0.
        T(1)=0.
        AUX=SQRT(S2*S2+S3*S3)
        T(2)=-S3/AUX
        T(3)= S2/AUX
      ELSEIF (IFACE.LE.NFAC12) THEN
        F(3)=0.
        S2=0.
        T(2)=0.
        AUX=SQRT(S1*S1+S3*S3)
        T(3)=-S1/AUX
        T(1)= S3/AUX
      ELSE
        F(4)=0.
        S3=0.
        T(3)=0.
        AUX=SQRT(S1*S1+S2*S2)
        T(1)=-S2/AUX
        T(2)= S1/AUX
      ENDIF
C
C     Correction of the isoline
      AUX=F(1)/AUX/AUX
      Y(1)=Y(1)-F(2)*AUX
      Y(2)=Y(2)-F(3)*AUX
      Y(3)=Y(3)-F(4)*AUX
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE OUTMS(X,Y,DERY,IHLF,NDIM,PRMT)
C
C-----------------------------------------------------------------------
C
      INTEGER IHLF,NDIM
      REAL X,Y(NDIM),DERY(NDIM),PRMT(*)
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
      EXTERNAL ISIDE
      INTEGER ISIDE
C     Auxiliary storage locations:
      REAL XTMP1,XTMP(3),DTMP(3),DUMMY(1),ERR,XINT1,DINT(3)
      REAL X01,X02,X1(3),X2(3)
      INTEGER IDUMMY,IY1(8),IY2(8),ISUR(2)
C.......................................................................
      IF (IHLF.GT.10) THEN
C       MODSRF-17
        CALL ERROR ('MODSRF-17: Wrong IHLF.')
C       This error should not appear. Contact the author.
      ENDIF
      IF (X.GT.PRMT(1)) THEN
        X01=0.
        X02=0.
C       Calling CDE along side 1:
        ERR=PRMT(4)
        XTMP1=X
        XTMP(1)=Y(1)
        XTMP(2)=Y(2)
        XTMP(3)=Y(3)
        DTMP(1)=DERY(1)
        DTMP(2)=DERY(2)
        DTMP(3)=DERY(3)
        IY1(4)=ISBO1
        IY1(5)=ICBO1
        IY1(6)=0
        CALL CDE(ISRFO,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY1,ERR,
     *    X0,X0,YOLD,DYOLD,X,Y,DERY,XTMP1,XTMP,DTMP,XINT1,X1,DINT)
        IF (IY1(6).NE.0) THEN
C         Interface crossed:
          X01=XINT1
        ENDIF
C
C       Calling CDE along side 2:
        XTMP1=X
        XTMP(1)=Y(1)
        XTMP(2)=Y(2)
        XTMP(3)=Y(3)
        DTMP(1)=DERY(1)
        DTMP(2)=DERY(2)
        DTMP(3)=DERY(3)
        IY2(4)=ISBO2
        IY2(5)=ICBO2
        IY2(6)=0
        CALL CDE(ISRFO,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY2,ERR,
     *    X0,X0,YOLD,DYOLD,X,Y,DERY,XTMP1,XTMP,DTMP,XINT1,X2,DINT)
        IF (IY2(6).NE.0) THEN
C         Interface crossed:
          X02=XINT1
        ENDIF
C
        IF ((X01.NE.0.).OR.(X02.NE.0.)) THEN
C         Interface crossed.
C
          IF ((X01.NE.0.).AND.(X02.NE.0.)) THEN
C           Choosing the nearer of the two crossed interfaces:
            IF (X01.LT.X02) THEN
              X02=0.
            ELSEIF (X02.LT.X01) THEN
              X01=0.
            ENDIF
          ENDIF
C
C         Recording indices of blocks computed by CDE:
          IF (X01.NE.0.) THEN
            ISBO1= IY1(4)
            ISRFN1=IY1(6)
            ISBN1= IY1(7)
            ICBN1= IY1(8)
            Y(1)=X1(1)
            Y(2)=X1(2)
            Y(3)=X1(3)
          ENDIF
          IF (X02.NE.0.) THEN
            ISBO2= IY2(4)
            ISRFN2=IY2(6)
            ISBN2= IY2(7)
            ICBN2= IY2(8)
            Y(1)=X2(1)
            Y(2)=X2(2)
            Y(3)=X2(3)
          ENDIF
C         Computing remaining indices of blocks:
          IF (X01.EQ.0.) THEN
            ISUR(1)=ISRFO
            ISUR(2)=ISRFN2
            CALL BLOCKS(X2,2,ISUR,ISBO2,IDUMMY,ISBO1,ICBO1)
            ISRFN1=ISRFN2
            IF (ISIDE(ISRFO,ISBN2).LT.2) THEN
              CALL BLOCKS(X2,2,ISUR,ISBN2,IDUMMY,ISBN1,ICBN1)
            ELSE
              ISBN1=ISBN2
              ICBN1=ICBN2
            ENDIF
          ENDIF
          IF (X02.EQ.0.) THEN
            ISUR(1)=ISRFO
            ISUR(2)=ISRFN1
            CALL BLOCKS(X1,2,ISUR,ISBO1,IDUMMY,ISBO2,ICBO2)
            ISRFN2=ISRFN1
            IF (ISIDE(ISRFO,ISBN1).LT.2) THEN
              CALL BLOCKS(X1,2,ISUR,ISBN1,IDUMMY,ISBN2,ICBN2)
            ELSE
              ISBN2=ISBN1
              ICBN2=ICBN1
            ENDIF
          ENDIF
          PRMT(5)=1.
        ENDIF
      ENDIF
C
C     Storing the old point:
      X0=X
      YOLD(1)=Y(1)
      YOLD(2)=Y(2)
      YOLD(3)=Y(3)
      DYOLD(1)=DERY(1)
      DYOLD(2)=DERY(2)
      DYOLD(3)=DERY(3)
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGLEG(ILEG,IPOIN1,IPOIN2)
C
C-----------------------------------------------------------------------
C
      INTEGER ILEG,IPOIN1,IPOIN2
C Input:
C     ILEG ... Index of gridleg.
C Output:
C     IPOIN1,IPOIN2 ... Indices of gridpoints forming the gridleg.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,ILEG1
C.......................................................................
C
      ILEG1=ILEG-1
      IF (ILEG.LE.NLEG1) THEN
        I31=ILEG1 / N11N2
        I21=(ILEG1 - I31*N11N2) / N11
        I1=ILEG - I31*N11N2 - I21*N11
        IPOIN1=I31*N1N2+I21*N1+I1
        IPOIN2=IPOIN1+1
      ELSEIF (ILEG.LE.NLEG12) THEN
        I31=(ILEG1-NLEG1) / N1N21
        I21=((ILEG1-NLEG1) - I31*N1N21) / N1
        I1=(ILEG-NLEG1) - I31*N1*N21 - I21*N1
        IPOIN1=I31*N1N2+I21*N1+I1
        IPOIN2=IPOIN1+N1
      ELSEIF (ILEG.LE.NLEG) THEN
        I31=(ILEG1-NLEG12) / N1N2
        I21=((ILEG1-NLEG12) - I31*N1N2) / N1
        I1=(ILEG-NLEG12) - I31*N1N2 - I21*N1
        IPOIN1=I31*N1N2+I21*N1+I1
        IPOIN2=IPOIN1+N1*N2
      ELSE
C       MODSRF-18
        CALL ERROR ('MODSRF-18: Wrong ILEG.')
C       This error should not appear. Contact the author.
      ENDIF
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGFAC(IFAC,ILEG1,ILEG2,ILEG3,ILEG4)
C
C-----------------------------------------------------------------------
C
      INTEGER IFAC,ILEG1,ILEG2,ILEG3,ILEG4
C     IFAC ... Index of gridface.
C     ILEG1,ILEG2,ILEG3,ILEG4 ... Indices of gridlegs forming
C              the gridface.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,IFAC1
C.......................................................................
C
      IFAC1=IFAC-1
      IF (IFAC.LE.NFAC1) THEN
        I31= IFAC1 / N1N21
        I21=(IFAC1 - I31*N1N21) / N1
        I1 = IFAC - I31*N1N21 - I21*N1
        ILEG1=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG4=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG2=ILEG4+N1
        ILEG3=ILEG1+N1N21
      ELSEIF (IFAC.LE.NFAC12) THEN
        I31= (IFAC1-NFAC1) / N11N2
        I21=((IFAC1-NFAC1) - I31*N11N2) / N11
        I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11
        ILEG1=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG4=         I31*N11N2+I21*N11+I1
        ILEG2=ILEG4+N11N2
        ILEG3=ILEG1+1
      ELSEIF (IFAC.LE.NFAC) THEN
        I31= (IFAC1-NFAC12) / N11N21
        I21=((IFAC1-NFAC12) - I31*N11N21) / N11
        I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11
        ILEG1=        I31*N11N2+I21*N11+I1
        ILEG4=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG2=ILEG4+1
        ILEG3=ILEG1+N11
      ELSE
C       MODSRF-19
        CALL ERROR ('MODSRF-19: Wrong IFAC.')
C       This error should not appear. Contact the author.
      ENDIF
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGCUB(ICUB,IFAC1,IFAC2,IFAC3,IFAC4,IFAC5,IFAC6)
C
C-----------------------------------------------------------------------
C
      INTEGER ICUB,IFAC1,IFAC2,IFAC3,IFAC4,IFAC5,IFAC6
C     ICUB ... Index of gridcube.
C     IFAC1,...,IFAC6 ... Indices of gridfaces forming the gridcube.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,ICUB1
C.......................................................................
C
      IF ((ICUB.LT.1).OR.(ICUB.GT.NCUB)) THEN
C       MODSRF-20
        CALL ERROR ('MODSRF-20: Wrong ICUB.')
C       This error should not appear. Contact the author.
      ENDIF
      ICUB1=ICUB-1
      I31= ICUB1 / N11N21
      I21=(ICUB1 - I31*N11N21) / N11
      I1=  ICUB  - I31*N11N21 - I21*N11
      IFAC1=       I31*N1N21 +I21*N1 +I1
      IFAC2=NFAC1 +I31*N11N2 +I21*N11+I1
      IFAC3=NFAC12+I31*N11N21+I21*N11+I1
      IFAC4=IFAC1+1
      IFAC5=IFAC2+N11
      IFAC6=IFAC3+N11N21
      RETURN
      END
C
C
C=======================================================================
C
      LOGICAL FUNCTION MSLPIF(XX,IFAC)
C
C-----------------------------------------------------------------------
C     Subroutine for decision, whether the point with coordinates XX
C     lies inside gridface IFAC.
C
      INTEGER IFAC
      REAL XX(3)
C Input:
C     XX   ... Coordinates of the point.
C     IFAC ... Index of the gridface.
C Output:
C     MSLPIF ... .TRUE. when the point lies inside the gridface.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,IFAC1,IP1,IP2
      INTEGER ILEG1,ILEG2
      REAL XA1,XA2,XA3,XB1,XB2,XB3,XC1,XC2,XC3
C.......................................................................
C
      IFAC1=IFAC-1
      IF (IFAC.LE.NFAC1) THEN
        I31= IFAC1 / N1N21
        I21=(IFAC1 - I31*N1N21) / N1
        I1 = IFAC - I31*N1N21 - I21*N1
        ILEG1=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG2=NLEG12 + I31*N1N2+I21*N1+I1 + N1
      ELSEIF (IFAC.LE.NFAC12) THEN
        I31= (IFAC1-NFAC1) / N11N2
        I21=((IFAC1-NFAC1) - I31*N11N2) / N11
        I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11
        ILEG1=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG2=         I31*N11N2+I21*N11+I1 + N11N2
      ELSEIF (IFAC.LE.NFAC) THEN
        I31= (IFAC1-NFAC12) / N11N21
        I21=((IFAC1-NFAC12) - I31*N11N21) / N11
        I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11
        ILEG1=        I31*N11N2+I21*N11+I1
        ILEG2=NLEG1 + I31*N1N21+I21*N1+I1 + 1
      ELSE
C       MODSRF-21
        CALL ERROR ('MODSRF-21: Wrong IFAC.')
C       This error should not appear. Contact the author.
      ENDIF
C
      MSLPIF=.TRUE.
C
      CALL MSGLEG(ILEG1,IP1,IP2)
      CALL MSCP(IP1,XA1,XA2,XA3)
      CALL MSCP(IP2,XB1,XB2,XB3)
      CALL MSGLEG(ILEG2,IP1,IP2)
      CALL MSCP(IP2,XC1,XC2,XC3)
      IF ((XX(1).LT.AMIN1(XA1,XB1,XC1)).OR.
     *    (XX(1).GT.AMAX1(XA1,XB1,XC1))) MSLPIF=.FALSE.
      IF ((XX(2).LT.AMIN1(XA2,XB2,XC2)).OR.
     *    (XX(2).GT.AMAX1(XA2,XB2,XC2))) MSLPIF=.FALSE.
      IF ((XX(3).LT.AMIN1(XA3,XB3,XC3)).OR.
     *    (XX(3).GT.AMAX1(XA3,XB3,XC3))) MSLPIF=.FALSE.
      RETURN
      END
C
C
C=======================================================================
C
      LOGICAL FUNCTION MSLPIL(XX,ILEG)
C
C-----------------------------------------------------------------------
C     Subroutine for decision, whether the point with coordinates XX
C     lies inside gridleg ILEG.
C
      INTEGER ILEG
      REAL XX(3)
C Input:
C     XX   ... Coordinates of the point.
C     ILEG ... Index of the gridleg.
C Output:
C     MSLPIF ... .TRUE. when the point lies inside the gridleg.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER IP1,IP2
      REAL XA1,XA2,XA3,XB1,XB2,XB3
C.......................................................................
C
      CALL MSGLEG(ILEG,IP1,IP2)
      CALL MSCP(IP1,XA1,XA2,XA3)
      CALL MSCP(IP2,XB1,XB2,XB3)
      MSLPIL=.TRUE.
      IF ((XX(1).LT.AMIN1(XA1,XB1)).OR.
     *    (XX(1).GT.AMAX1(XA1,XB1))) MSLPIL=.FALSE.
      IF ((XX(2).LT.AMIN1(XA2,XB2)).OR.
     *    (XX(2).GT.AMAX1(XA2,XB2))) MSLPIL=.FALSE.
      IF ((XX(3).LT.AMIN1(XA3,XB3)).OR.
     *    (XX(3).GT.AMAX1(XA3,XB3))) MSLPIL=.FALSE.
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGCGP(ICUB,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8)
C
C-----------------------------------------------------------------------
C     Subroutine for debugging purposes,  gives indices of points
C     of cube ICUB.
C
C
      INTEGER ICUB,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8
C     ICUB ... Index of gridcube.
C     IP1,...,IP8 ... Indices of gridpoints forming the gridcube.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,ICUB1
C.......................................................................
C
      IF ((ICUB.LT.1).OR.(ICUB.GT.NCUB)) THEN
C       MODSRF-22
        CALL ERROR ('MODSRF-22: Wrong ICUB.')
C       This error should not appear. Contact the author.
      ENDIF
      ICUB1=ICUB-1
      I31= ICUB1 / N11N21
      I21=(ICUB1 - I31*N11N21) / N11
      I1=  ICUB  - I31*N11N21 - I21*N11
      IP1=I31*N1N2+I21*N1+I1
      IP2=IP1+N1
      IP3=IP2+N1N2
      IP4=IP1+N1N2
      IP5=IP1+1
      IP6=IP5+N1
      IP7=IP6+N1N2
      IP8=IP5+N1N2
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSXFAC(IFAC)
C
C-----------------------------------------------------------------------
C     Subroutine for debugging purposes,  gives coordinates of points
C     of face IFAC.
C
      INTEGER IFAC,ILEG1,ILEG2,ILEG3,ILEG4
C     IFAC ... Index of gridface.
C     ILEG1,ILEG2,ILEG3,ILEG4 ... Indices of gridlegs forming
C              the gridface.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,IFAC1,IP1,IP2
      REAL X1,X2,X3
C.......................................................................
C
      IFAC1=IFAC-1
      IF (IFAC.LE.NFAC1) THEN
        I31= IFAC1 / N1N21
        I21=(IFAC1 - I31*N1N21) / N1
        I1 = IFAC - I31*N1N21 - I21*N1
        ILEG1=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG4=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG2=ILEG4+N1
        ILEG3=ILEG1+N1N21
      ELSEIF (IFAC.LE.NFAC12) THEN
        I31= (IFAC1-NFAC1) / N11N2
        I21=((IFAC1-NFAC1) - I31*N11N2) / N11
        I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11
        ILEG1=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG4=         I31*N11N2+I21*N11+I1
        ILEG2=ILEG4+N11N2
        ILEG3=ILEG1+1
      ELSEIF (IFAC.LE.NFAC) THEN
        I31= (IFAC1-NFAC12) / N11N21
        I21=((IFAC1-NFAC12) - I31*N11N21) / N11
        I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11
        ILEG1=        I31*N11N2+I21*N11+I1
        ILEG4=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG2=ILEG4+1
        ILEG3=ILEG1+N11
      ELSE
C       MODSRF-23
        CALL ERROR ('MODSRF-23: Wrong IFAC.')
C       This error should not appear. Contact the author.
      ENDIF
      CALL MSGLEG(ILEG1,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      CALL MSGLEG(ILEG2,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      CALL MSGLEG(ILEG3,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      CALL MSGLEG(ILEG4,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSCP(IPOIN,X1,X2,X3)
C
C-----------------------------------------------------------------------
C
      INTEGER IPOIN
      REAL X1,X2,X3
C     IPOIN ... Index of a point.
C     X1,X2,X3 ... Coordinates of the point.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I11,IPOIN1
C.......................................................................
C
      IPOIN1=IPOIN-1
      I31= IPOIN1 / N1N2
      I21=(IPOIN1 - I31*N1N2) / N1
      I11= IPOIN  - I31*N1N2 - I21*N1 - 1
      X1=O1+FLOAT(I11)*D1
      X2=O2+FLOAT(I21)*D2
      X3=O3+FLOAT(I31)*D3
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSEROR(IERR)
C
C-----------------------------------------------------------------------
C
      INTEGER IERR
C     IERR ... Index of the error.
C.......................................................................
      IF (IERR.EQ.1) THEN
C       MODSRF-24
        CALL ERROR('MODSRF-24: Small array (I)RAM.')
C       Try to enlarge the dimension MRAM in the file
C       ram.inc. If this does not help,
C       contact the author.
      ELSEIF (IERR.EQ.2) THEN
C       MODSRF-25
        CALL ERROR('MODSRF-25: Small array JPT.')
C       Try to enlarge the dimension MJPT at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.3) THEN
C       MODSRF-26
        CALL ERROR('MODSRF-26: Small array IPTE.')
C       Try to enlarge the dimension MIPTE at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.4) THEN
C       MODSRF-27
        CALL ERROR('MODSRF-27: Small array JCN.')
C       Try to enlarge the dimension MJCN at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.5) THEN
C       MODSRF-28
        CALL ERROR('MODSRF-28: Small array JSP.')
C       Try to enlarge the dimension MJSP at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.6) THEN
C       MODSRF-29
        CALL ERROR('MODSRF-29: Small array JPOL.')
C       Try to enlarge the dimension MJPOL at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.7) THEN
C       MODSRF-30
        CALL ERROR('MODSRF-30: Small array JPOL.')
C       Try to enlarge the dimension MLJPOL at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.8) THEN
C       MODSRF-05
        CALL ERROR('MODSRF-05: Interface not found.')
C       This might happen when the whole gridleg coincides with an
C       interface. If possible, slightly change the origin of
C       the grid. If this does not help, contact the author.
      ELSE
C       MODSRF-31
        CALL ERROR('MODSRF-31: Wrong IERR.')
C       This subroutine was invocated with wrong value of IERR.
C       This error should not appear, please contact the author.
      ENDIF
      END
C
C=======================================================================
C
      INCLUDE 'model.for'
C     model.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'rkgs.for'
C     rkgs.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'error.for'
C     error.for
C
C=======================================================================
C