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