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: 5.40
C Date: 2000, May 19
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.......................................................................
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 component of the normal to the interface
C               (or to the 2-D slice) at the vertex.
C             'NORM2' ... Second component of the normal.
C             'NORM3' ... Third 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=(D1+D2+D3)/3000. (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
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,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 ... 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
      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),MJSP,NJSP,
     * MJPOL,NJPOL,MLJPOL,NLJPOL,NIND
      PARAMETER (MJCN=50,MJSP=100,MJPOL=10,MLJPOL=10)
      INTEGER JCN(3,MJCN),JSP(3,MJSP),JPOL(0:MLJPOL,MJPOL),IND(MJCN)
C     JCN(1 to 3,i) ... Connection i of the face: address of the first
C       point, address of the second point, index of the interface.
C     NNJCN(1:7) ... Number of connections for first gridface, second
C     gridface, ... , sixth gridface, between edges.
C     JSP(1 to 3,i) ... Array of starting points as they go along
C       a polygon: previous point, point, following point.
C     JPOL(0,i) ... Number of points in polygon i.
C     JPOL(1 to JPOL(0,i),i) ... Points in the polygon.
      INTEGER MIPTE
      PARAMETER (MIPTE=5)
      INTEGER IPTE(MIPTE),NIPTE
C     IPTE(i) ... address of the point to be used as starting point
C       when searching for the edge. IPTE(i).LT.0 means that the point
C       is to be used in reverse order than it is stored in array IRAM.
      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-32
        CALL ERROR('MODSRF-32: 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-01
        CALL ERROR('MODSRF-01: Wrong specification of the grid.')
C       This specification of the grid may cause problems. Please,
C       specify D1,D2,D3 nonzero and N1,N2,N3 greater than 0. Di may
C       equal zero in case that corresponding Ni equals 1.
      ENDIF
C
C     Reading the data for the model:
      CALL RSEP3T('MODEL',FMOD,' ')
      IF (FMOD.EQ.' ') THEN
C       MODSRF-33
        CALL ERROR('MODSRF-33: 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-02
            CALL ERROR('MODSRF-02: Wrong value of COLUMN.')
C           See allowed values of COLUMNii in the
C           description of file SEP.
          ENDIF
          I1=I1+1
          IF (I1.GT.69) THEN
            CALL RSEP3T('COLUMN70',TEXT,' ')
            IF (TEXT.NE.' ') THEN
C             MODSRF-03
              CALL ERROR('MODSRF-03: More than 69 COLUMNs.')
C             Currently up to 69 values of COLUMNii may be specified.
C             See allowed values of COLUMNii in the
C             description of file SEP.
            ENDIF
          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:
      CALL RSEP3R('ERRSRF',ERRSRF,(D1+D2+D3)/3000.)
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-04
              CALL ERROR ('MODSRF-04: More than 100 intersections.')
C             More than 100 points of intersection of the
C             gridline element with the structural interfaces.
C             Check the input data and then contact the author.
            ENDIF
C
C           Determining the interface between points XOLD, XNEW:
            IY(4)=ISBOLD
            IY(5)=ICBOLD
            IY(6)=0
            XTMP(I4)=XNEW(I4)
            XOLD1=XOLD(I4)
            XNEW1=XNEW(I4)
            XTMP1=XNEW1
            CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERRSRF,
     *               XOLD1,XOLD1,XOLD,DOLD,XNEW1,XNEW,DNEW,
     *                     XTMP1,XTMP,DTMP,XINT1,XINT,DINT)
            IF ((ICBOLD.NE.ICBNEW).AND.(IY(6).EQ.0)) THEN
C             The interface was not found between the points with
C             different ICB. This may happen when the NEW point is
C             situated directly at the interface.
C             Repeating the search in reverse direction:
              DOLD(I4)=-DOLD(I4)
              DNEW(I4)=-DNEW(I4)
              DTMP(I4)=-DTMP(I4)
              XTMP1=XOLD1
              XTMP(1)=XOLD(1)
              XTMP(2)=XOLD(2)
              XTMP(3)=XOLD(3)
              IY(4)=ISBNEW
              IY(5)=ICBNEW
              CALL CDE(0,0,IDUMMY,0,IDUMMY,DUMMY,1,3,3,IY,ERRSRF,
     *               (XTMP1-ERRSRF),XNEW1,XNEW,DNEW,XOLD1,XOLD,DOLD,
     *                     XTMP1,XTMP,DTMP,XINT1,XINT,DINT)
              IF (IY(6).EQ.0) THEN
C               The interface was not found even in reverse direction.
C               This may happen when the whole gridleg coincides with an
C               interface:
                CALL BLOCK(XOLD,0,ISBOLD,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBOLD) CALL MSEROR(8)
                CALL BLOCK(XOLD,0,ISBNEW,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBNEW) CALL MSEROR(8)
                CALL BLOCK(XNEW,0,ISBNEW,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBNEW) CALL MSEROR(8)
                CALL BLOCK(XNEW,0,ISBOLD,IDUMMY,IB,IDUMMY)
                IF (IB.NE.ISBOLD) CALL MSEROR(8)
C               Whole gridleg coincides with an interface:
                CALL SEPAR(ISBNEW,ISBOLD,IDUMMY,IY(6))
                IF (IDUMMY(1).LT.1) CALL MSEROR(8)
                IY(7)=ISBOLD
                IY(8)=ICBOLD
              ENDIF
              IF (IY(8).NE.ICBOLD) THEN
C               There should be only one interface between the points.
C               
C               MODSRF-06
                CALL ERROR('MODSRF-06: More than one interface.')
C               This error should not appear. Contact the author.
              ENDIF
C             Reverting the points to the original order:
              DOLD(I4)=-DOLD(I4)
              DNEW(I4)=-DNEW(I4)
              DTMP(I4)=-DTMP(I4)
              XTMP1=XNEW1
              XTMP(1)=XNEW(1)
              XTMP(2)=XNEW(2)
              XTMP(3)=XNEW(3)
              ISBOLD=IY(7)
              ICBOLD=IY(8)
              ISBNEW=IY(4)
              ICBNEW=IY(5)
              IY(4)=ISBOLD
              IY(5)=ICBOLD
              IY(6)=-IY(6)
              IY(7)=ISBNEW
              IY(8)=ICBNEW
            ENDIF
            IF (IY(6).NE.0) THEN
C             Structural interface, recording the point of intersection:
C             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-07
                CALL ERROR('MODSRF-07: Point not on the gridleg.')
C               This error should not appear. Contact the author.
              ENDIF
              IF (NPOI+NQPOI.GT.MRAM) CALL MSEROR(1)
              RAM(NPOI+1)=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)=PRMT(2)/30.
      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 (IABS(JPT(4,I2)).EQ.IABS(JPT(4,I3))) THEN
C               Same interface - connection found:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                IF (JPT(4,I2).LT.0) THEN
                  IRAM(NCON+1)=IABS(JPT(1,I2))
                  IRAM(NCON+2)=IABS(JPT(1,I3))
                  IRAM(NCON+3)=JPT(4,I2)
                  NCON=NCON+3
                ELSE
                  IRAM(NCON+1)=IABS(JPT(1,I3))
                  IRAM(NCON+2)=IABS(JPT(1,I2))
                  IRAM(NCON+3)=JPT(4,I3)
                  NCON=NCON+3
                ENDIF
                GOTO 45
              ENDIF
            ENDIF
  36      CONTINUE
          DO 37, I3=IRAM(OFAC+I1-1),NCON-3,3
            IF ((IRAM(I3+1).EQ.IABS(JPT(1,I2))).OR.
     *          (IRAM(I3+2).EQ.IABS(JPT(1,I2))))
C             The point I2 is already connected:
     *        GOTO 45
  37      CONTINUE
C         Point I2 is not connected, looking for point of intersection
C         of structural edges with gridface. (Following the intersection
C         of the interface with the gridface to the edge.)
          IFACE=I1
          IPTE(1)=JPT(1,I2)
          NIPTE=1
C         Loop for identification of all edges connected with point I2:
  39      CONTINUE
C           Initiating the search for the edge:
            I3=IABS(IPTE(NIPTE))-NQPOI
            X1=    RAM(I3+1)
            X2=    RAM(I3+2)
            X3=    RAM(I3+3)
            IF (IPTE(NIPTE).LT.0) THEN
              ICBO1=IRAM(I3+8)
              ISBO1=IRAM(I3+7)
              ISRFO=-IRAM(I3+6)
              ISBO2=IRAM(I3+5)
              ICBO2=IRAM(I3+4)
            ELSE
              ICBO1=IRAM(I3+4)
              ISBO1=IRAM(I3+5)
              ISRFO=IRAM(I3+6)
              ISBO2=IRAM(I3+7)
              ICBO2=IRAM(I3+8)
            ENDIF
            NIPTE=NIPTE-1
C           Computing the point of intersection of the edge
C           with the gridface:
            CALL RKGS(PRMT,XX,DTMP,3,I7,FCTMS,OUTMS,AUX)
            IF (PRMT(5).NE.1.) THEN
C             MODSRF-08
              CALL ERROR('MODSRF-08: PRMT(5) not equal to 1.')
C             This error should not appear. Contact the author.
            ENDIF
            IF (.NOT.MSLPIF(XX,IFACE)) THEN
C             MODSRF-09
              CALL ERROR('MODSRF-09: Point outside the gridface.')
C             This error should not appear. Contact the author.
            ENDIF
C           Storing the first point on the edge:
            IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
            RAM(NPOI+1)=X1
            RAM(NPOI+2)=X2
            RAM(NPOI+3)=X3
            IRAM(NPOI+4)=ICBO1
            IRAM(NPOI+5)=ISBO1
            IRAM(NPOI+6)=ISRFO
            IRAM(NPOI+7)=ISBO2
            IRAM(NPOI+8)=ICBO2
            NPOI=NPOI+NQPOI
            J1=NPOI
C           Recording the connection:
            IF (NCON+3.GT.MRAM) CALL MSEROR(1)
            IF (ISRFO.LT.0) THEN
              IRAM(NCON+1)=IABS(IPTE(NIPTE+1))
              IRAM(NCON+2)=NPOI
              IRAM(NCON+3)=ISRFO
            ELSE
              IRAM(NCON+1)=NPOI
              IRAM(NCON+2)=IABS(IPTE(NIPTE+1))
              IRAM(NCON+3)=-ISRFO
            ENDIF
            NCON=NCON+3
C
            IF (ICBO1.NE.ICBN1) THEN
C             Storing the second point on the edge:
              IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
              RAM(NPOI+1)=X1
              RAM(NPOI+2)=X2
              RAM(NPOI+3)=X3
              IRAM(NPOI+4)=ICBN1
              IRAM(NPOI+5)=ISBN1
              IRAM(NPOI+6)=-ISRFN1
              IRAM(NPOI+7)=ISBO1
              IRAM(NPOI+8)=ICBO1
              NPOI=NPOI+NQPOI
              IF (NCUB.EQ.0) THEN
C               'Connection' between first and second point on the edge:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                J2=NPOI-NQPOI
                IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
                J3=NPOI
                IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
                IRAM(NCON+1)=J2
                IRAM(NCON+2)=J3
                IRAM(NCON+3)=0
                NCON=NCON+3
              ENDIF
              DO 40, I3=1,NJPT
                IF ((ICBN1.EQ.JPT(2,I3)).AND.
     *              (ISBN1.EQ.JPT(3,I3)).AND.
     *              (ISBO1.EQ.JPT(5,I3)).AND.
     *              (ICBO1.EQ.JPT(6,I3)).AND.
     *              (IABS(ISRFN1).EQ.IABS(JPT(4,I3)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (JPT(4,I3).LT.0) THEN
                    IRAM(NCON+1)=IABS(JPT(1,I3))
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=JPT(4,I3)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=IABS(JPT(1,I3))
                    IRAM(NCON+3)=-JPT(4,I3)
                  ENDIF
                  NCON=NCON+3
                  GOTO 405
                ENDIF
  40          CONTINUE
              DO 401, I3=1,NIPTE
                I4=IABS(IPTE(I3))-NQPOI
                IF ((ICBN1.EQ.IRAM(I4+8)).AND.
     *              (ISBN1.EQ.IRAM(I4+7)).AND.
     *              (ISBO1.EQ.IRAM(I4+5)).AND.
     *              (ICBO1.EQ.IRAM(I4+4)).AND.
     *              (IABS(ISRFN1).EQ.IABS(IRAM(I4+6)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (IRAM(I4+6).LT.0) THEN
                    IRAM(NCON+1)=I4+NQPOI
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=IRAM(I4+6)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=I4+NQPOI
                    IRAM(NCON+3)=-IRAM(I4+6)
                  ENDIF
                  NCON=NCON+3
                  GOTO 405
                ENDIF
  401         CONTINUE
C             More than one edge, current point will become a starting
C             point of the search for new edge:
              NIPTE=NIPTE+1
              IF (NIPTE.GT.MIPTE) CALL MSEROR(3)
              IPTE(NIPTE)=-NPOI
  405         CONTINUE
            ENDIF
C
            IF (ICBN2.NE.ICBN1) THEN
C             Storing the third point on the edge:
              IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
              RAM(NPOI+1)=X1
              RAM(NPOI+2)=X2
              RAM(NPOI+3)=X3
              IRAM(NPOI+4)=ICBN2
              IRAM(NPOI+5)=ISBN2
              IRAM(NPOI+6)=-ISRFO
              IRAM(NPOI+7)=ISBN1
              IRAM(NPOI+8)=ICBN1
              NPOI=NPOI+NQPOI
              IF (NCUB.EQ.0) THEN
C               'Connection' between second and third point on the edge:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                J2=NPOI-NQPOI
                IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
                J3=NPOI
                IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
                IRAM(NCON+1)=J2
                IRAM(NCON+2)=J3
                IRAM(NCON+3)=0
                NCON=NCON+3
              ENDIF
              DO 42, I3=1,NJPT
                IF ((ICBN2.EQ.JPT(2,I3)).AND.
     *              (ISBN2.EQ.JPT(3,I3)).AND.
     *              (ISBN1.EQ.JPT(5,I3)).AND.
     *              (ICBN1.EQ.JPT(6,I3)).AND.
     *              (IABS(ISRFO).EQ.IABS(JPT(4,I3)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (JPT(4,I3).LT.0) THEN
                    IRAM(NCON+1)=IABS(JPT(1,I3))
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=JPT(4,I3)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=IABS(JPT(1,I3))
                    IRAM(NCON+3)=-JPT(4,I3)
                  ENDIF
                  NCON=NCON+3
                  GOTO 425
                ENDIF
  42          CONTINUE
              DO 421, I3=1,NIPTE
                I4=IABS(IPTE(I3))-NQPOI
                IF ((ICBN2.EQ.IRAM(I4+8)).AND.
     *              (ISBN2.EQ.IRAM(I4+7)).AND.
     *              (ISBN1.EQ.IRAM(I4+5)).AND.
     *              (ICBN1.EQ.IRAM(I4+4)).AND.
     *              (IABS(ISRFO).EQ.IABS(IRAM(I4+6)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (IRAM(I4+6).LT.0) THEN
                    IRAM(NCON+1)=I4+NQPOI
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=IRAM(I4+6)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=I4+NQPOI
                    IRAM(NCON+3)=-IRAM(I4+6)
                  ENDIF
                  NCON=NCON+3
                  GOTO 425
                ENDIF
  421         CONTINUE
C             More than one edge, current point will become a starting
C             point of the search for new edge:
              NIPTE=NIPTE+1
              IF (NIPTE.GT.MIPTE) CALL MSEROR(3)
              IPTE(NIPTE)=-NPOI
  425         CONTINUE
            ENDIF
C
            IF (ICBO2.NE.ICBN2) THEN
C             Storing the fourth point on the edge:
              IF (NPOI+NQPOI.GT.OCON) CALL MSEROR(1)
              RAM(NPOI+1)=X1
              RAM(NPOI+2)=X2
              RAM(NPOI+3)=X3
              IRAM(NPOI+4)=ICBO2
              IRAM(NPOI+5)=ISBO2
              IRAM(NPOI+6)=ISRFN2
              IRAM(NPOI+7)=ISBN2
              IRAM(NPOI+8)=ICBN2
              NPOI=NPOI+NQPOI
              IF (NCUB.EQ.0) THEN
C               'Connection' between third and fourth point on the edge:
                IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                J2=NPOI-NQPOI
                IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
                J3=NPOI
                IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
                IRAM(NCON+1)=J2
                IRAM(NCON+2)=J3
                IRAM(NCON+3)=0
                NCON=NCON+3
              ENDIF
              DO 41, I3=1,NJPT
                IF ((ICBO2.EQ.JPT(2,I3)).AND.
     *              (ISBO2.EQ.JPT(3,I3)).AND.
     *              (ISBN2.EQ.JPT(5,I3)).AND.
     *              (ICBN2.EQ.JPT(6,I3)).AND.
     *              (IABS(ISRFN2).EQ.IABS(JPT(4,I3)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (JPT(4,I3).LT.0) THEN
                    IRAM(NCON+1)=IABS(JPT(1,I3))
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=JPT(4,I3)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=IABS(JPT(1,I3))
                    IRAM(NCON+3)=-JPT(4,I3)
                  ENDIF
                  NCON=NCON+3
                  GOTO 415
                ENDIF
  41          CONTINUE
              DO 411, I3=1,NIPTE
                I4=IABS(IPTE(I3))-NQPOI
                IF ((ICBO2.EQ.IRAM(I4+8)).AND.
     *              (ISBO2.EQ.IRAM(I4+7)).AND.
     *              (ISBN2.EQ.IRAM(I4+5)).AND.
     *              (ICBN2.EQ.IRAM(I4+4)).AND.
     *              (IABS(ISRFN2).EQ.IABS(IRAM(I4+6)))) THEN
C                 Recording the connection:
                  IF (NCON+3.GT.MRAM) CALL MSEROR(1)
                  IF (IRAM(I4+6).LT.0) THEN
                    IRAM(NCON+1)=I4+NQPOI
                    IRAM(NCON+2)=NPOI
                    IRAM(NCON+3)=IRAM(I4+6)
                  ELSE
                    IRAM(NCON+1)=NPOI
                    IRAM(NCON+2)=I4+NQPOI
                    IRAM(NCON+3)=-IRAM(I4+6)
                  ENDIF
                  NCON=NCON+3
                  GOTO 415
                ENDIF
  411         CONTINUE
C             More than one edge, current point will become a starting
C             point of the search for new edge:
              NIPTE=NIPTE+1
              IF (NIPTE.GT.MIPTE) CALL MSEROR(3)
              IPTE(NIPTE)=-NPOI
  415         CONTINUE
            ENDIF
C
            IF (NCUB.EQ.0) THEN
C             'Connection' between fourth and first point on the edge:
              IF (NCON+3.GT.MRAM) CALL MSEROR(1)
              J2=NPOI
              IF (IRAM(J2-NQPOI+6).LT.0) J2=-J2
              J3=J1
              IF (IRAM(J3-NQPOI+6).GT.0) J3=-J3
              IRAM(NCON+1)=J2
              IRAM(NCON+2)=J3
              IRAM(NCON+3)=0
              NCON=NCON+3
            ENDIF
          IF (NIPTE.GE.1)
C           Continuing with the search for next edge:
     *      GOTO 39
C
C         No more edges, the point is connected.
  45      CONTINUE
C         Continuing with next point of the gridface:
  46    CONTINUE
        IRAM(OFAC+I1)=NCON
  49  CONTINUE
      NPOIE=NPOI
C
C
      IF (NCUB.NE.0) THEN
C       Loop along all gridcubes,
C       merging the connections between the points on structural
C       interfaces into polygons approximating the interfaces.
        WRITE(*,'(A)')
     *'+MODSRF: Making polygons approximating the interfaces.'
        NPOL=MRAM+1
        DO 100, I1=1,NCUB
C         Indices of gridfaces of the gridcube:
          CALL MSGCUB(I1,IF1,IF2,IF3,IF4,IF5,IF6)
C
C         Forming array with connections for the gridcube:
          NJCN=0
          NNJCN(0)=NJCN
          DO 51, I2=IRAM(OFAC+IF1-1),IRAM(OFAC+IF1)-3,3
C           Connections of the gridface IF1:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+1)
            JCN(2,NJCN)=IRAM(I2+2)
            JCN(3,NJCN)=IRAM(I2+3)
  51      CONTINUE
          NNJCN(1)=NJCN
          DO 52, I2=IRAM(OFAC+IF2-1),IRAM(OFAC+IF2)-3,3
C           Connections of the gridface IF2:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+1)
            JCN(2,NJCN)=IRAM(I2+2)
            JCN(3,NJCN)=IRAM(I2+3)
  52      CONTINUE
          NNJCN(2)=NJCN
          DO 53, I2=IRAM(OFAC+IF3-1),IRAM(OFAC+IF3)-3,3
C           Connections of the gridface IF3:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+1)
            JCN(2,NJCN)=IRAM(I2+2)
            JCN(3,NJCN)=IRAM(I2+3)
  53      CONTINUE
          NNJCN(3)=NJCN
          DO 54, I2=IRAM(OFAC+IF4-1),IRAM(OFAC+IF4)-3,3
C           Connections of the gridface IF4:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+2)
            JCN(2,NJCN)=IRAM(I2+1)
            JCN(3,NJCN)=IRAM(I2+3)
  54      CONTINUE
          NNJCN(4)=NJCN
          DO 55, I2=IRAM(OFAC+IF5-1),IRAM(OFAC+IF5)-3,3
C           Connections of the gridface IF5:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+2)
            JCN(2,NJCN)=IRAM(I2+1)
            JCN(3,NJCN)=IRAM(I2+3)
  55      CONTINUE
          NNJCN(5)=NJCN
          DO 56, I2=IRAM(OFAC+IF6-1),IRAM(OFAC+IF6)-3,3
C           Connections of the gridface IF6:
            NJCN=NJCN+1
            IF (NJCN.GT.MJCN) CALL MSEROR(4)
            JCN(1,NJCN)=IRAM(I2+2)
            JCN(2,NJCN)=IRAM(I2+1)
            JCN(3,NJCN)=IRAM(I2+3)
  56      CONTINUE
          NNJCN(6)=NJCN
          DO 58, I2=1,NNJCN(6)
C           Connections between points on edges:
            IF (JCN(2,I2).GT.NPOIL) THEN
              ICBO1=IRAM(JCN(2,I2)-NQPOI+4)
              ICBO2=IRAM(JCN(2,I2)-NQPOI+8)
              LCN=.FALSE.
              DO 57, I3=1,NNJCN(6)
                IF (I3.NE.I2) THEN
                  IF ((JCN(1,I3).GT.NPOIL).AND.
     *                (IABS(JCN(3,I2)).EQ.IABS(JCN(3,I3)))) THEN
                    ICBN1=IRAM(JCN(1,I3)-NQPOI+4)
                    ICBN2=IRAM(JCN(1,I3)-NQPOI+8)
                    IF (((ICBO1.EQ.ICBN1).AND.(ICBO2.EQ.ICBN2)).OR.
     *                  ((ICBO2.EQ.ICBN1).AND.(ICBO1.EQ.ICBN2))) THEN
                      NJCN=NJCN+1
                      IF (NJCN.GT.MJCN) CALL MSEROR(4)
                      JCN(1,NJCN)=JCN(2,I2)
                      JCN(2,NJCN)=JCN(1,I3)
                      JCN(3,NJCN)=JCN(3,I3)
                      LCN=.TRUE.
                    ENDIF
                  ENDIF
                ENDIF
  57          CONTINUE
              IF (.NOT.LCN) THEN
C               Connection for edge not found, this connection is to be
C               canceled by means of creation of linear 'triangle':
                NJCN=NJCN+1
                IF (NJCN.GT.MJCN) CALL MSEROR(4)
                JCN(1,NJCN)=JCN(2,I2)
                JCN(2,NJCN)=JCN(1,I2)
                JCN(3,NJCN)=JCN(3,I2)
              ENDIF
            ENDIF
  58      CONTINUE
          NNJCN(7)=NJCN
C
C         Forming array with starting points
C         (from point, point, to point)
          NJSP=0
          DO 60, I2=1,NJCN
            DO 59, I3=1,NJCN
              IF ((JCN(2,I2).EQ.JCN(1,I3)).AND.
     *            (JCN(2,I3).NE.JCN(1,I2))) THEN
                NJSP=NJSP+1
                IF (NJSP.GT.MJSP) CALL MSEROR(5)
                JSP(1,NJSP)=JCN(1,I2)
                JSP(2,NJSP)=JCN(2,I2)
                JSP(3,NJSP)=JCN(2,I3)
              ENDIF
  59        CONTINUE
  60      CONTINUE
C
C         Looking for polygons on the faces of the cube.
          NJPOL=0
C         Loop for faces of the cube:
          DO 70, I2=1,6
C           Initiating indices of connections:
            DO 61, I3=1,MJCN
              IND(I3)=0
  61        CONTINUE
            NIND=NNJCN(I2)-NNJCN(I2-1)
C           Loop along connections:
            I3=1
  62        CONTINUE
              IF (NJPOL+1.GT.MJPOL) CALL MSEROR(6)
C             Search for first unused connection:
              IF (I3.LE.NIND-2) THEN
                IF (IND(I3).GE.0) THEN
C                 Unused connection, recording first two points
C                 of the polygon:
                  NLJPOL=2
                  JPOL(1,NJPOL+1)=JCN(1,NNJCN(I2-1)+I3)
                  JPOL(2,NJPOL+1)=JCN(2,NNJCN(I2-1)+I3)
                  IND(I3)=-1
                ELSE
                  I3=I3+1
                  GOTO 62
                ENDIF
C               Search for next point of polygon:
  63            CONTINUE
                DO 69, I4=I3+1,NIND
                  IF ( JCN(1,NNJCN(I2-1)+I4)       .EQ.
     *                 JPOL(NLJPOL,NJPOL+1) ) THEN
C                   Recording next point of polygon:
                    NLJPOL=NLJPOL+1
                    IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7)
                    JPOL(NLJPOL,NJPOL+1)=JCN(2,NNJCN(I2-1)+I4)
                    IND(I4)=-1
                    DO 68, I5=NLJPOL-3,1,-1
C                     Examining whether the polygon is closed:
                      IF (JPOL(I5,NJPOL+1).EQ.JPOL(NLJPOL,NJPOL+1)) THEN
C                       Polygon is closed.
C                       Recording the polygon:
                        NJPOL=NJPOL+1
                        NLJPOL=NLJPOL-I5
                        JPOL(0,NJPOL)=NLJPOL
                        DO 65, I6=1,NLJPOL
                          JPOL(I6,NJPOL)=JPOL(I5+I6,NJPOL)
  65                    CONTINUE
C                       Removing used starting points:
                        DO 67, I6=1,NLJPOL
                          IF (I6.EQ.1) THEN
                            J1=JPOL(NLJPOL,NJPOL)
                          ELSE
                            J1=JPOL(I6-1,NJPOL)
                          ENDIF
                          J2=JPOL(I6,NJPOL)
                          IF (I6.EQ.NJPOL) THEN
                            J3=JPOL(1,NJPOL)
                          ELSE
                            J3=JPOL(I6+1,NJPOL)
                          ENDIF
                          DO 66, I7=1,NJSP
                            IF ((J1.EQ.JSP(1,I7)).AND.(J2.EQ.JSP(2,I7))
     *                          .AND.(J3.EQ.JSP(3,I7))) THEN
                              JSP(1,I7)=0
                              JSP(2,I7)=0
                              JSP(3,I7)=0
                            ENDIF
  66                      CONTINUE
  67                    CONTINUE
C                       Continuing with next polygon:
                        GOTO 62
                      ENDIF
  68                CONTINUE
                    GOTO 63
                  ENDIF
  69            CONTINUE
C               Polygon is opened, searching for other polygon:
                GOTO 62
              ENDIF
C             No more unused connections,
C             no more polygons on this face of the cube.
C           End of the loop along connections.
  70      CONTINUE
C
C         If there are two or more polygons on the faces of the cube,
C         the polygons are to be omitted:
          IF (NJPOL.GE.2) THEN
C           MODSRF-10
            CALL ERROR('MODSRF-10: Two or more polygons on gridface.')
C           This situation cannot be handled by current version of
C           MODSRF. Try to reduce the gridstep (D1,D2,D3).
          ENDIF
C
C         If there is one polygon on the faces of the cube,
C         the polygon is to be recorded as a polygon approximating
C         a structural interface:
          IF (NJPOL.EQ.1) THEN
C           Searching, whether the polygon is to be recorded:
C           Loop for already recorded polygons:
            I2=MRAM+1
  73        CONTINUE
              IF (I2.GT.NPOL) THEN
                J1=IRAM(I2-1)
                I2=I2-J1-1
                IF (J1.EQ.JPOL(0,1)) THEN
                  DO 74, I3=1,J1
                    IF (IRAM(I2+I3).NE.JPOL(I3,1)) GOTO 73
  74              CONTINUE
C                 The polygon is already recorded:
                  GOTO 77
                ELSE
                  GOTO 73
                ENDIF
              ENDIF
C           End of the loop for already recorded polygons.
C           Recording the polygon:
            IF (NPOL-JPOL(0,1)-1.LE.NCON) CALL MSEROR(1)
            NPOL=NPOL-1
            IRAM(NPOL)=JPOL(0,1)
            DO 75, I2=1,JPOL(0,1)
              NPOL=NPOL-1
              IRAM(NPOL)=JPOL(I2,1)
  75        CONTINUE
  77        CONTINUE
          ENDIF
C
C         Recording polygons approximating structural interfaces:
          NJPOL=0
  80      CONTINUE
C         Loop for adding new polygons
C           Initializing a polygon:
            DO 82, I2=1,NJSP
C             Prefering points on edges:
              IF (JSP(1,I2).GT.NPOIL) THEN
                NJPOL=NJPOL+1
                IF (NJPOL.GT.MJPOL) CALL MSEROR(6)
                NLJPOL=3
                JPOL(0,NJPOL)=0
                JPOL(1,NJPOL)=JSP(1,I2)
                JPOL(2,NJPOL)=JSP(2,I2)
                JPOL(3,NJPOL)=JSP(3,I2)
                JSP(1,I2)=0
                JSP(2,I2)=0
                JSP(3,I2)=0
                GOTO  85
              ENDIF
  82        CONTINUE
            DO 84, I2=1,NJSP
C             No points on edges, taking remaining points:
              IF (JSP(1,I2).NE.0) THEN
                NJPOL=NJPOL+1
                IF (NJPOL.GT.MJPOL) CALL MSEROR(6)
                NLJPOL=3
                JPOL(0,NJPOL)=0
                JPOL(1,NJPOL)=JSP(1,I2)
                JPOL(2,NJPOL)=JSP(2,I2)
                JPOL(3,NJPOL)=JSP(3,I2)
                JSP(1,I2)=0
                JSP(2,I2)=0
                JSP(3,I2)=0
                GOTO  85
              ENDIF
  84        CONTINUE
C           No other point to initialize a polygon:
            GOTO 95
C
  85        CONTINUE
C           Loop for adding a points to the polygon:
              DO 93, I2=1,NJSP
                IF ((JSP(1,I2).EQ.JPOL(NLJPOL-1,NJPOL)).AND.
     *              (JSP(2,I2).EQ.JPOL(NLJPOL  ,NJPOL))) THEN
                  DO 87, I3=2,NLJPOL-2
                    IF (JSP(3,I2).EQ.JPOL(I3,NJPOL)) THEN
C                     This point would cause overlooping of the polygon:
                      JSP(1,I2)=0
                      JSP(2,I2)=0
                      JSP(3,I2)=0
                      DO 86, I4=1,NJSP
                        IF ((JSP(1,I4).EQ.JPOL(NLJPOL,NJPOL)).AND.
     *                      (JSP(2,I4).EQ.JPOL(I3    ,NJPOL)).AND.
     *                      (JSP(3,I4).EQ.JPOL(I3+1  ,NJPOL))) THEN
                          JSP(1,I4)=0
                          JSP(2,I4)=0
                          JSP(3,I4)=0
                          GOTO 85
                        ENDIF
  86                  CONTINUE
C                     
C                     MODSRF-11
                      CALL ERROR('MODSRF-11: Disorder in array JSP.')
C                     This error should not appear. Contact the author.
                    ENDIF
  87              CONTINUE
                  IF (JSP(3,I2).EQ.JPOL(1,NJPOL)) THEN
C                   This point would close the polygon:
                    DO 91, I3=I2+1,NJSP
C                     Looking for other possible continuation, which
C                     would make polygon longer:
                      IF ((JSP(1,I3).EQ.JPOL(NLJPOL-1,NJPOL)).AND.
     *                    (JSP(2,I3).EQ.JPOL(NLJPOL  ,NJPOL)).AND.
     *                    (JSP(3,I3).NE.JPOL(1,NJPOL))) THEN
C                       Other possible continuation of the polygon:
                        DO 90, I4=2,NLJPOL-2
                          IF (JSP(3,I3).EQ.JPOL(I4,NJPOL)) THEN
C                           This point would cause overlooping of the
C                           polygon:
                            JSP(1,I3)=0
                            JSP(2,I3)=0
                            JSP(3,I3)=0
                            DO 89, I5=1,NJSP
                              IF (   (JSP(1,I5).EQ.JPOL(NLJPOL,NJPOL))
     *                          .AND.(JSP(2,I5).EQ.JPOL(I4    ,NJPOL))
     *                          .AND.(JSP(3,I5).EQ.JPOL(I4+1  ,NJPOL)))
     *                          THEN
                                JSP(1,I5)=0
                                JSP(2,I5)=0
                                JSP(3,I5)=0
                                GOTO 85
                              ENDIF
  89                        CONTINUE
C                           
C                           MODSRF-12
                            CALL ERROR
     *                      ('MODSRF-12: Disorder in array JSP.')
C                           This error should not appear.
C                           Contact the author.
                          ENDIF
  90                    CONTINUE
C                       Continuation I2 is to be omitted because
C                       continuation I3 makes the polygon longer:
                        JSP(1,I2)=0
                        JSP(2,I2)=0
                        JSP(3,I2)=0
                        IF (JPOL(0,NJPOL).EQ.0) THEN
C                         Removing the second point which closes the
C                         polygon from JSP:
                          DO 88, I4=1,NJSP
                            IF ((JSP(1,I4).EQ.JPOL(NLJPOL,NJPOL)).AND.
     *                          (JSP(2,I4).EQ.JPOL(1,NJPOL)).AND.
     *                          (JSP(3,I4).EQ.JPOL(2,NJPOL))) THEN
                              JSP(1,I4)=0
                              JSP(2,I4)=0
                              JSP(3,I4)=0
                              JPOL(0,NJPOL)=0
                              GOTO 85
                            ENDIF
  88                      CONTINUE
C                         
C                         MODSRF-13
                          CALL ERROR
     *                    ('MODSRF-13: Disorder in array JSP.')
C                         This error should not appear.
C                         Contact the author.
                        ELSE
                          JPOL(0,NJPOL)=0
                          GOTO 85
                        ENDIF
                      ENDIF
  91                CONTINUE
C                   No other continuation, polygon is to be closed:
                    JSP(1,I2)=0
                    JSP(2,I2)=0
                    JSP(3,I2)=0
                    IF (JPOL(0,NJPOL).EQ.0) THEN
C                     First point which closes the polygon.
                      JPOL(0,NJPOL)=NLJPOL
                      GOTO 85
                    ELSE
C                     Second point which closes the polygon.
                      GOTO 80
                    ENDIF
                  ELSE
C                   Adding a point to the polygon:
                    NLJPOL=NLJPOL+1
                    IF (NLJPOL.GT.MLJPOL) CALL MSEROR(7)
                    JPOL(NLJPOL,NJPOL)=JSP(3,I2)
                    JSP(1,I2)=0
                    JSP(2,I2)=0
                    JSP(3,I2)=0
                    GOTO 85
                  ENDIF
                ELSEIF ((JSP(1,I2).EQ.JPOL(NLJPOL,NJPOL)).AND.
     *                  (JSP(2,I2).EQ.JPOL(1,NJPOL)).AND.
     *                  (JSP(3,I2).EQ.JPOL(2,NJPOL))) THEN
C                 The polygon is closed.
                  JSP(1,I2)=0
                  JSP(2,I2)=0
                  JSP(3,I2)=0
                  IF (JPOL(0,NJPOL).EQ.0) THEN
C                   First point which closes the polygon.
                    JPOL(0,NJPOL)=NLJPOL
                    GOTO 85
                  ELSE
C                   Second point which closes the polygon.
                    GOTO 80
                  ENDIF
                ENDIF
  93          CONTINUE
C             MODSRF-14
              CALL ERROR ('MODSRF-14: Polygon opened.')
C             This error should not appear. Contact the author.
C           End of the loop for adding new points to the polygon.
C         End of the loop for adding new polygons.
C
  95      CONTINUE
C         Writing polygons to IRAM:
          DO 97, I2=1,NJPOL
            IF (NPOL-JPOL(0,I2)-1.LE.NCON) CALL MSEROR(1)
            NPOL=NPOL-1
            IRAM(NPOL)=JPOL(0,I2)
            DO 96, I3=1,JPOL(0,I2)
              NPOL=NPOL-1
              IRAM(NPOL)=JPOL(I3,I2)
  96        CONTINUE
  97      CONTINUE
          NJPOL=0
C
C       End of the loop over gridcubes.
 100    CONTINUE
C
C
      ELSE
C       Loop over 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-15
              CALL ERROR ('MODSRF-15: Polygon not closed.')
C             This error should not appear. Contact the author.
            ENDIF
 130      CONTINUE
C       End of the loop along gridfaces of the 2-D slice.
 200    CONTINUE
C
C       Preparing the normal to the 2-D slice:
        F(2)=0.
        F(3)=0.
        F(4)=0.
        IF (N1.EQ.1) F(2)=1.
        IF (N2.EQ.1) F(3)=1.
        IF (N3.EQ.1) F(4)=1.
C
      ENDIF
C
C
C     Storing points on structural interfaces and polygons:
      WRITE(*,'(A)')
     *'+MODSRF: Writing.                                     '
C
C     Storing points:
      OPEN(LU,FILE=VRTX,FORM='FORMATTED')
      WRITE(LU,'(A)') '/'
      NPOINT=(NPOI-OPOI)/NQPOI
C     Length of the names of points:
      LEN1=LENGTH(TEXTP)
      LEN2=0
 202  CONTINUE
        IF (NPOINT.GE.10**LEN2) THEN
          LEN2=LEN2+1
          GOTO 202
        ENDIF
      LENG=LEN1+LEN2
C     Loop over points:
      DO 218, I1=1,NPOINT
C       Address in RAM just before the current point:
        J1=OPOI+(I1-1)*NQPOI
C       Name of the point:
        DO 204, I2=0,LEN2-1
          TEXTP(LENG-I2:LENG-I2)=
     *    CHAR(ICHAR('0')+MOD(I1,10**(I2+1))/10**I2)
 204    CONTINUE
C       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-16
              CALL ERROR ('MODSRF-16: Disorder in polygons.')
C             This error should not appear. Contact the author.
            ENDIF
            DO 224, I2=I1-1,I1-J1,-1
C             Index of the point (from 1 to NPOINT):
              J2=(IRAM(I2)-OPOI)/NQPOI
C             Name of the point:
              I4=I1-I2
              DO 222, I3=0,LEN2-1
                TEXTS(I4)(1:LEN1)=TEXTP(1:LEN1)
                TEXTS(I4)(LENG-I3:LENG-I3)=
     *          CHAR(ICHAR('0')+MOD(J2,10**(I3+1))/10**I3)
 222          CONTINUE
 224        CONTINUE
            WRITE(LU2,FORMA2)
     *           ('''',TEXTS(I2)(1:LENG),'''',I2=1,J1),'/'
          ENDIF
        ENDIF
        I1=I1-(J1+1)
        GOTO 220
      ENDIF
C     End of the loop over polygons.
      IF (PLGN.NE.' ') THEN
        WRITE(LU1,'(A)') '/'
        CLOSE(LU1)
      ENDIF
      IF (PLGNS.NE.' ') THEN
        WRITE(LU2,'(A)') '/'
        CLOSE(LU2)
      ENDIF
C
      WRITE(*,'(A)')
     *'+MODSRF: 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 Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
      EXTERNAL ISIDE
      INTEGER ISIDE
C     Auxiliary storage locations:
      REAL XTMP1,XTMP(3),DTMP(3),DUMMY(1),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-17
        CALL ERROR ('MODSRF-17: Wrong IHLF.')
C       This error should not appear. Contact the author.
      ENDIF
      IF (X.GT.PRMT(1)) THEN
        X01=0.
        X02=0.
C       Calling CDE along side 1:
        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
          PRMT(5)=1.
        ENDIF
      ENDIF
C
C     Storing the old point:
      X0=X
      YOLD(1)=Y(1)
      YOLD(2)=Y(2)
      YOLD(3)=Y(3)
      DYOLD(1)=DERY(1)
      DYOLD(2)=DERY(2)
      DYOLD(3)=DERY(3)
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGLEG(ILEG,IPOIN1,IPOIN2)
C
C-----------------------------------------------------------------------
C
      INTEGER ILEG,IPOIN1,IPOIN2
C Input:
C     ILEG ... Index of gridleg.
C Output:
C     IPOIN1,IPOIN2 ... Indices of gridpoints forming the gridleg.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,ILEG1
C.......................................................................
C
      ILEG1=ILEG-1
      IF (ILEG.LE.NLEG1) THEN
        I31=ILEG1 / N11N2
        I21=(ILEG1 - I31*N11N2) / N11
        I1=ILEG - I31*N11N2 - I21*N11
        IPOIN1=I31*N1N2+I21*N1+I1
        IPOIN2=IPOIN1+1
      ELSEIF (ILEG.LE.NLEG12) THEN
        I31=(ILEG1-NLEG1) / N1N21
        I21=((ILEG1-NLEG1) - I31*N1N21) / N1
        I1=(ILEG-NLEG1) - I31*N1*N21 - I21*N1
        IPOIN1=I31*N1N2+I21*N1+I1
        IPOIN2=IPOIN1+N1
      ELSEIF (ILEG.LE.NLEG) THEN
        I31=(ILEG1-NLEG12) / N1N2
        I21=((ILEG1-NLEG12) - I31*N1N2) / N1
        I1=(ILEG-NLEG12) - I31*N1N2 - I21*N1
        IPOIN1=I31*N1N2+I21*N1+I1
        IPOIN2=IPOIN1+N1*N2
      ELSE
C       MODSRF-18
        CALL ERROR ('MODSRF-18: Wrong ILEG.')
C       This error should not appear. Contact the author.
      ENDIF
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGFAC(IFAC,ILEG1,ILEG2,ILEG3,ILEG4)
C
C-----------------------------------------------------------------------
C
      INTEGER IFAC,ILEG1,ILEG2,ILEG3,ILEG4
C     IFAC ... Index of gridface.
C     ILEG1,ILEG2,ILEG3,ILEG4 ... Indices of gridlegs forming
C              the gridface.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,IFAC1
C.......................................................................
C
      IFAC1=IFAC-1
      IF (IFAC.LE.NFAC1) THEN
        I31= IFAC1 / N1N21
        I21=(IFAC1 - I31*N1N21) / N1
        I1 = IFAC - I31*N1N21 - I21*N1
        ILEG1=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG4=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG2=ILEG4+N1
        ILEG3=ILEG1+N1N21
      ELSEIF (IFAC.LE.NFAC12) THEN
        I31= (IFAC1-NFAC1) / N11N2
        I21=((IFAC1-NFAC1) - I31*N11N2) / N11
        I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11
        ILEG1=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG4=         I31*N11N2+I21*N11+I1
        ILEG2=ILEG4+N11N2
        ILEG3=ILEG1+1
      ELSEIF (IFAC.LE.NFAC) THEN
        I31= (IFAC1-NFAC12) / N11N21
        I21=((IFAC1-NFAC12) - I31*N11N21) / N11
        I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11
        ILEG1=        I31*N11N2+I21*N11+I1
        ILEG4=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG2=ILEG4+1
        ILEG3=ILEG1+N11
      ELSE
C       MODSRF-19
        CALL ERROR ('MODSRF-19: Wrong IFAC.')
C       This error should not appear. Contact the author.
      ENDIF
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGCUB(ICUB,IFAC1,IFAC2,IFAC3,IFAC4,IFAC5,IFAC6)
C
C-----------------------------------------------------------------------
C
      INTEGER ICUB,IFAC1,IFAC2,IFAC3,IFAC4,IFAC5,IFAC6
C     ICUB ... Index of gridcube.
C     IFAC1,...,IFAC6 ... Indices of gridfaces forming the gridcube.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,ICUB1
C.......................................................................
C
      IF ((ICUB.LT.1).OR.(ICUB.GT.NCUB)) THEN
C       MODSRF-20
        CALL ERROR ('MODSRF-20: Wrong ICUB.')
C       This error should not appear. Contact the author.
      ENDIF
      ICUB1=ICUB-1
      I31= ICUB1 / N11N21
      I21=(ICUB1 - I31*N11N21) / N11
      I1=  ICUB  - I31*N11N21 - I21*N11
      IFAC1=       I31*N1N21 +I21*N1 +I1
      IFAC2=NFAC1 +I31*N11N2 +I21*N11+I1
      IFAC3=NFAC12+I31*N11N21+I21*N11+I1
      IFAC4=IFAC1+1
      IFAC5=IFAC2+N11
      IFAC6=IFAC3+N11N21
      RETURN
      END
C
C
C=======================================================================
C
      LOGICAL FUNCTION MSLPIF(XX,IFAC)
C
C-----------------------------------------------------------------------
C     Subroutine for decision, whether the point with coordinates XX
C     lies inside gridface IFAC.
C
      INTEGER IFAC
      REAL XX(3)
C Input:
C     XX   ... Coordinates of the point.
C     IFAC ... Index of the gridface.
C Output:
C     MSLPIF ... .TRUE. when the point lies inside the gridface.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,IFAC1,IP1,IP2
      INTEGER ILEG1,ILEG2
      REAL XA1,XA2,XA3,XB1,XB2,XB3,XC1,XC2,XC3
C.......................................................................
C
      IFAC1=IFAC-1
      IF (IFAC.LE.NFAC1) THEN
        I31= IFAC1 / N1N21
        I21=(IFAC1 - I31*N1N21) / N1
        I1 = IFAC - I31*N1N21 - I21*N1
        ILEG1=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG2=NLEG12 + I31*N1N2+I21*N1+I1 + N1
      ELSEIF (IFAC.LE.NFAC12) THEN
        I31= (IFAC1-NFAC1) / N11N2
        I21=((IFAC1-NFAC1) - I31*N11N2) / N11
        I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11
        ILEG1=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG2=         I31*N11N2+I21*N11+I1 + N11N2
      ELSEIF (IFAC.LE.NFAC) THEN
        I31= (IFAC1-NFAC12) / N11N21
        I21=((IFAC1-NFAC12) - I31*N11N21) / N11
        I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11
        ILEG1=        I31*N11N2+I21*N11+I1
        ILEG2=NLEG1 + I31*N1N21+I21*N1+I1 + 1
      ELSE
C       MODSRF-21
        CALL ERROR ('MODSRF-21: Wrong IFAC.')
C       This error should not appear. Contact the author.
      ENDIF
C
      MSLPIF=.TRUE.
C
      CALL MSGLEG(ILEG1,IP1,IP2)
      CALL MSCP(IP1,XA1,XA2,XA3)
      CALL MSCP(IP2,XB1,XB2,XB3)
      CALL MSGLEG(ILEG2,IP1,IP2)
      CALL MSCP(IP2,XC1,XC2,XC3)
      IF ((XX(1).LT.AMIN1(XA1,XB1,XC1)).OR.
     *    (XX(1).GT.AMAX1(XA1,XB1,XC1))) MSLPIF=.FALSE.
      IF ((XX(2).LT.AMIN1(XA2,XB2,XC2)).OR.
     *    (XX(2).GT.AMAX1(XA2,XB2,XC2))) MSLPIF=.FALSE.
      IF ((XX(3).LT.AMIN1(XA3,XB3,XC3)).OR.
     *    (XX(3).GT.AMAX1(XA3,XB3,XC3))) MSLPIF=.FALSE.
      RETURN
      END
C
C
C=======================================================================
C
      LOGICAL FUNCTION MSLPIL(XX,ILEG)
C
C-----------------------------------------------------------------------
C     Subroutine for decision, whether the point with coordinates XX
C     lies inside gridleg ILEG.
C
      INTEGER ILEG
      REAL XX(3)
C Input:
C     XX   ... Coordinates of the point.
C     ILEG ... Index of the gridleg.
C Output:
C     MSLPIF ... .TRUE. when the point lies inside the gridleg.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER IP1,IP2
      REAL XA1,XA2,XA3,XB1,XB2,XB3
C.......................................................................
C
      CALL MSGLEG(ILEG,IP1,IP2)
      CALL MSCP(IP1,XA1,XA2,XA3)
      CALL MSCP(IP2,XB1,XB2,XB3)
      MSLPIL=.TRUE.
      IF ((XX(1).LT.AMIN1(XA1,XB1)).OR.
     *    (XX(1).GT.AMAX1(XA1,XB1))) MSLPIL=.FALSE.
      IF ((XX(2).LT.AMIN1(XA2,XB2)).OR.
     *    (XX(2).GT.AMAX1(XA2,XB2))) MSLPIL=.FALSE.
      IF ((XX(3).LT.AMIN1(XA3,XB3)).OR.
     *    (XX(3).GT.AMAX1(XA3,XB3))) MSLPIL=.FALSE.
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSGCGP(ICUB,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8)
C
C-----------------------------------------------------------------------
C     Subroutine for debugging purposes,  gives indices of points
C     of cube ICUB.
C
C
      INTEGER ICUB,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8
C     ICUB ... Index of gridcube.
C     IP1,...,IP8 ... Indices of gridpoints forming the gridcube.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,ICUB1
C.......................................................................
C
      IF ((ICUB.LT.1).OR.(ICUB.GT.NCUB)) THEN
C       MODSRF-22
        CALL ERROR ('MODSRF-22: Wrong ICUB.')
C       This error should not appear. Contact the author.
      ENDIF
      ICUB1=ICUB-1
      I31= ICUB1 / N11N21
      I21=(ICUB1 - I31*N11N21) / N11
      I1=  ICUB  - I31*N11N21 - I21*N11
      IP1=I31*N1N2+I21*N1+I1
      IP2=IP1+N1
      IP3=IP2+N1N2
      IP4=IP1+N1N2
      IP5=IP1+1
      IP6=IP5+N1
      IP7=IP6+N1N2
      IP8=IP5+N1N2
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSXFAC(IFAC)
C
C-----------------------------------------------------------------------
C     Subroutine for debugging purposes,  gives coordinates of points
C     of face IFAC.
C
      INTEGER IFAC,ILEG1,ILEG2,ILEG3,ILEG4
C     IFAC ... Index of gridface.
C     ILEG1,ILEG2,ILEG3,ILEG4 ... Indices of gridlegs forming
C              the gridface.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I1,IFAC1,IP1,IP2
      REAL X1,X2,X3
C.......................................................................
C
      IFAC1=IFAC-1
      IF (IFAC.LE.NFAC1) THEN
        I31= IFAC1 / N1N21
        I21=(IFAC1 - I31*N1N21) / N1
        I1 = IFAC - I31*N1N21 - I21*N1
        ILEG1=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG4=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG2=ILEG4+N1
        ILEG3=ILEG1+N1N21
      ELSEIF (IFAC.LE.NFAC12) THEN
        I31= (IFAC1-NFAC1) / N11N2
        I21=((IFAC1-NFAC1) - I31*N11N2) / N11
        I1 = (IFAC-NFAC1) - I31*N11N2 - I21*N11
        ILEG1=NLEG12 + I31*N1N2+I21*N1+I1
        ILEG4=         I31*N11N2+I21*N11+I1
        ILEG2=ILEG4+N11N2
        ILEG3=ILEG1+1
      ELSEIF (IFAC.LE.NFAC) THEN
        I31= (IFAC1-NFAC12) / N11N21
        I21=((IFAC1-NFAC12) - I31*N11N21) / N11
        I1 = (IFAC-NFAC12) - I31*N11N21 - I21*N11
        ILEG1=        I31*N11N2+I21*N11+I1
        ILEG4=NLEG1 + I31*N1N21+I21*N1+I1
        ILEG2=ILEG4+1
        ILEG3=ILEG1+N11
      ELSE
C       MODSRF-23
        CALL ERROR ('MODSRF-23: Wrong IFAC.')
C       This error should not appear. Contact the author.
      ENDIF
      CALL MSGLEG(ILEG1,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      CALL MSGLEG(ILEG2,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      CALL MSGLEG(ILEG3,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      CALL MSGLEG(ILEG4,IP1,IP2)
      CALL MSCP(IP1,X1,X2,X3)
      CALL MSCP(IP2,X1,X2,X3)
C
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSCP(IPOIN,X1,X2,X3)
C
C-----------------------------------------------------------------------
C
      INTEGER IPOIN
      REAL X1,X2,X3
C     IPOIN ... Index of a point.
C     X1,X2,X3 ... Coordinates of the point.
C.......................................................................
C Common block /RAMC/ to store the information about points on
C structural interfaces and common block /MSC/ to store auxiliary
C quantities:
      INCLUDE 'modsrf.inc'
C modsrf.inc.
C ...........................
C     Auxiliary storage locations:
      INTEGER I31,I21,I11,IPOIN1
C.......................................................................
C
      IPOIN1=IPOIN-1
      I31= IPOIN1 / N1N2
      I21=(IPOIN1 - I31*N1N2) / N1
      I11= IPOIN  - I31*N1N2 - I21*N1 - 1
      X1=O1+FLOAT(I11)*D1
      X2=O2+FLOAT(I21)*D2
      X3=O3+FLOAT(I31)*D3
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE MSEROR(IERR)
C
C-----------------------------------------------------------------------
C
      INTEGER IERR
C     IERR ... Index of the error.
C.......................................................................
      IF (IERR.EQ.1) THEN
C       MODSRF-24
        CALL ERROR('MODSRF-24: Small array (I)RAM.')
C       Try to enlarge the dimension MRAM in the file
C       ram.inc. If this does not help,
C       contact the author.
      ELSEIF (IERR.EQ.2) THEN
C       MODSRF-25
        CALL ERROR('MODSRF-25: Small array JPT.')
C       Try to enlarge the dimension MJPT at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.3) THEN
C       MODSRF-26
        CALL ERROR('MODSRF-26: Small array IPTE.')
C       Try to enlarge the dimension MIPTE at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.4) THEN
C       MODSRF-27
        CALL ERROR('MODSRF-27: Small array JCN.')
C       Try to enlarge the dimension MJCN at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.5) THEN
C       MODSRF-28
        CALL ERROR('MODSRF-28: Small array JSP.')
C       Try to enlarge the dimension MJSP at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.6) THEN
C       MODSRF-29
        CALL ERROR('MODSRF-29: Small array JPOL.')
C       Try to enlarge the dimension MJPOL at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.7) THEN
C       MODSRF-30
        CALL ERROR('MODSRF-30: Small array JPOL.')
C       Try to enlarge the dimension MLJPOL at the beginning of this file.
C       If this does not help, contact the author.
      ELSEIF (IERR.EQ.8) THEN
C       MODSRF-05
        CALL ERROR('MODSRF-05: Interface not found.')
C       This might happen when the whole gridleg coincides with an
C       interface. If possible, slightly change the origin of
C       the grid. If this does not help, contact the author.
      ELSE
C       MODSRF-31
        CALL ERROR('MODSRF-31: Wrong IERR.')
C       This subroutine was invocated with wrong value of IERR.
C       This error should not appear, please contact the author.
      ENDIF
      END
C
C=======================================================================
C
      INCLUDE '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