C 
C Program INVPTS to calculate the derivatives of functions, describing
C interfaces or material parameters, with respect to the model B-spline
C coefficients
C
C Version: 6.00
C Date: 2006, March 3
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C Program INVPTS assumes all model parameters (coefficients) stored in
C the common block /VALC/ as in the submitted versions of user-defined
C model specification FORTRAN77 source code files 'srfc.for', 'parm.for'
C and 'val.for'.  Thus, unlike the other parts of the complete ray
C tracing, the INVTT program cannot work with user's modifications of
C subroutines SRFC1, SRFC2, PARM1, and PARM2.
C
C The material parameters at the given points are converted into such
C their powers which are interpolated by B-splines.  Then the inversion
C of the functions of coordinates is linear and one iteration yields
C the exact solution within the rounding errors.  On the other hand,
C inversion of a material parameter depending on another material
C parameters, e.g., RHO=W(VP(X1,X2,X3)), require the second iteration.
C More than one iteration may also be required to reduce the rounding
C errors.
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:
C     MODEL='string'... Name of the input data file describing the
C             model.  For the description of the data refer to
C             subroutine file 'model.for'.
C             Default: MODEL='model.dat'
C Input files with the values to be used as the data for the inversion:
C     PTS='string'... String with the name of the input data file
C             containing the individual given points with the values to
C             be fitted by the model.  Points situated outside the model
C             box are not considered.
C             This program may be executed several times to take into
C             account several files with individual points, or the same
C             file with several values at each point (e.g., both P and
C             S wave velocities at each point).
C             If the filename is blank, no individual points are given.
C             Description of file PTS
C             Default: PTS=' '
C     LIN='string'... Alternative to PTS, if the given points are
C             arranged into lines.  Differs from PTS just by the
C             form of the file.
C             If the filename is blank, no lines of points are given.
C             Default: LIN=' '
C     GRD='string'... String with the name of the input data file
C             containing the grided values to be fitted by the model.
C             The grid may also contain undefined values.
C             If the whole grid cell centred at a point is situated
C             outside the model box, the point is not considered.
C             This program may be executed several times to take into
C             account several grids of values.
C             If the filename is blank, no gridded values are assumed.
C             Default: GRD=' '
C     GRDERR='string'... String with the name of the input data file
C             containing the grided standard deviations of the given
C             values.
C             Considered only if GRD is specified.
C             If the filename is blank, unit standard deviations are
C             assumed.
C             The standard deviations are multiplied by the value of
C             ERRMUL, described later on.
C             Default: GRDERR=' '
C     GRDICB='string'... String with the name of the input data file
C             containing the grided indices of the complex blocks
C             corresponding to file GRD, if file GRD contains the
C             material parameters of different complex blocks
C             (considered only if GRD is specified).
C             If GRDICB=' ', the index of the surface or complex block
C             corresponding to file GRD is common for all gridpoints
C             and is given by parameter INDFUN.
C             Default: GRDICB=' '
C     For general description of the files with gridded data refer to
C     file forms.htm.
C Data specifying the function corresponding to the given values:
C     ICLASS=integer... Class of model parameters to be inverted:
C             ICLASS=0: All model parameters are inverted.
C             ICLASS=1: Only model parameters describing interfaces are
C               inverted.
C             ICLASS=2: Only model parameters describing material
C               parameters are inverted.
C             Default: ICLASS=0
C     MPAR=integer... Material parameter corresponding to the given
C             values.
C             MPAR=0:  The given values correspond to the functions
C               describing the surfaces covering structural interfaces.
C             MPAR=positive integer:  The given values describe a
C               material parameter:
C               MPAR=1: P wave velocity,
C               MPAR=2: S wave velocity,
C               MPAR=3: density,
C               MPAR=4: P wave loss factor,
C               MPAR=5: S wave loss factor,
C               MPAR=6 to 26: Reduced (i.e., divided by the density)
C                 anisotropic elastic parameters A11, A12, A22, A13,
C                 A23, A33, A14, A24, A34, A44, A15, A25, A35, A45, A55,
C                 A16, A26, A36, A46, A56 or A66.
C               MPAR=27-47: Reduced (i.e., divided by the density)
C                 imaginary parts of anisotropic elastic parameters Q11,
C                 Q12, Q22, Q13, Q23, Q33, Q14, Q24, Q34, Q44, Q15, Q25,
C                 Q35, Q45, Q55, Q16, Q26, Q36, Q46, Q56 or Q66.
C             Default: MPAR=0
C     POWERM=real... The given values correspond to the POWERMth power
C             of the material parameter determined by MPAR.
C             Default: POWERM=1.
C     KOLFUN=integer... Specifies the column in input file 'PTS' or
C             'LIN' containing the index of the surface (for MPAR=0) or
C             the complex block.
C             If KOLFUN=0, the index is common for all points and is
C             given by parameter INDFUN.
C             Default: KOLFUN=0
C     INDFUN=integer... Index of the surface (for MPAR=0) or of the
C             complex block to which the given values correspond.
C             If KOLFUN=0 for individual points or GRDICB=' ' for
C             gridded values, INDFUN must be specified and positive.
C             Otherwise, INDFUN is not used.
C             Default: INDFUN=0
C Data specifying the values to be processed and their accuracy:
C     KOLUMN=integer... Specifies the column in input file 'PTS' or
C             'LIN' containing the given values.  Note that the first
C             3 columns contain coordinates of the points.
C             If KOLUMN=0, zero values are assumed (option often used
C             for the points at interfaces, where the function
C             describing the corresponding surface should be zero).
C             Default: KOLUMN=0
C     KOLERR=integer... Specifies the column in input file 'PTS' or
C             'LIN' containing the standard deviations of the given
C             values.
C             If KOLERR=0, unit standard deviations are assumed.
C             Note that the standard deviations are multiplied by
C             ERRMUL, described below.
C             Default: KOLERR=0
C     ERRMUL=real... Multiplication factor for the standard deviations
C             of the given values.  Often used with GRDERR=' ' to
C             specify the standard deviation constant over the grid or
C             with KOLERR=0 to specify the standard deviation common for
C             all points.
C             Default: ERRMUL=1.
C     Data specifying the grid parameters of the gridded values
C     (considered only if GRD is specified):
C     O1=real, O2=real, O3=real ... Coordinates of the origin of the
C             grid.
C             Defaults: O1=0., O2=0., O3=0.
C     D1=real, D2=real, D3=real ... Grid intervals.
C             Default: D1=1., D2=1., D3=1.
C     N1=positive integer, N2=positive integer, N3=positive integer ...
C             Numbers of gridpoints.
C             Default: N1=1, N2=1, N3=1
C Data specifying input/output files:
C     M1='string'... Name of the output file containing number M1 of
C             model parameters (a single integer).  The same file may be
C             generated by programs 'invsoft.for' and 'invtt.for'.
C             The file is not generated if the value of M1 is blank.
C             Default: M1=' '
C             Note: Default of 'invsoft.for' is M1='m1.out',
C                   default of 'invtt.for' is M1=' '.
C     M2='string'... Name of the output file containing number M2 of
C             accumulated values (a single integer).  M2 is M2IN
C             increased by the number of given values.
C             The file is not generated if the value of M2 is blank.
C             Default: M2='m2.out'
C     M2IN='string'... Name of the input file containing number M2IN of
C             values already processed by programs 'invpts.for' or
C             'invtt.for', i.e., the name of output file M2 from the
C             previous execution of 'invpts.for' or 'invtt.for'.
C             The corresponding output values of this program will be
C             appended after M1*M2IN input values of file GM1 and after
C             M2IN input values of files GM2, GM3 and DM1, respectively.
C             Default: M2IN=' ' means M2IN=0 (no previous values).
C Data specifying input/output files with matrix and vector components:
C     If the values in input files PTS, LIN or GRD correspond to
C     a material parameter (specified by MPAR=1 to MPAR=47),
C     they may be converted:  Input files PTS, LIN or GRD contain
C     the values of the material parameter powered to POWERM.
C     On the other hand, the material parameter powered to POWER is
C     interpolated by B-splines.  Exponent POWER is specified in the
C     input data for the model.  The data
C     for inversion of the material parameter are thus the values
C     given in input files PTS, LIN or GRD powered to POWER/POWERM.
C     This conversion of input data makes the inversion of material
C     parameters linear in most cases.  Then one iteration yields the
C     exact solution within the rounding errors.  On the other hand,
C     inversion of a material parameter depending on another material
C     parameter, e.g., RHO=W(VP(X1,X2,X3)), requires the second
C     iteration.
C     GM1='string'... String with the name of the input/output data file
C             to accumulate the matrix of the derivatives of
C             (POWER/POWERM)th powers of M2 given values with respect to
C             M1 model coefficients.
C             The matrix has M1*M2IN components on input and M1*M2
C             components on output.  The formatted file is composed of
C             real-valued matrix components to be read at once by the
C             list-directed input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: GM1=' '
C     GM2='string'... String with the name of the input/output data file
C             containing the vector composed of differences between the
C             given values powered to POWER/POWERM and the corresponding
C             values calculated in the model.
C             The vector has M2IN components on input and M2 components
C             on output.  The formatted file is composed of real-valued
C             vector components to be read at once by list-directed
C             input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: GM2=' '
C     GM3='string'... String with the name of the input/output data file
C             containing the vector composed of the values calculated in
C             the model.  The vector has M2IN components on input and M2
C             components on output.  The formatted file is composed of
C             real-valued vector components to be read at once by
C             list-directed input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: GM3=' '
C     DM1='string'... String with the name of the input/output data file
C             containing the diagonal matrix composed of the variances
C             (squares of standard deviations) of the given values
C             powered to POWER/POWERM.
C             The diagonal has M2IN components on input and M2
C             components on output.  The formatted file is composed
C             of real-valued diagonal components to be read at once by
C             list-directed input.
C             If the filename is blank, no file is read nor written.
C             The file is not read, just written if M2IN=0.
C             Default: DM1=' '
C     For general description of the files with matrices refer to file
C     forms.htm.
C Form of the files with matrices:
C     FORMM='string' ... Form of the files with matrices. Allowed values
C             are FORMM='formatted' and FORMM='unformatted'. If the form
C             differs for input and for output files, FORMMR and FORMMW
C             should be used instead of FORMM.
C             Default: FORMM='formatted'
C     FORMMR='string'... Form of the files with matrices to be read.
C             Default: FORMMR=FORMM
C     FORMMW='string'... Form of the files with matrices to be written.
C             Default: FORMMW=FORMM
C
C                                                     
C Input file PTS with the points:
C (1) None to several strings terminated by / (a slash)
C (2) For each point data:
C (2.1) 'NAME',X1,X2,X3,V1,...,VN,/
C     'NAME'... Name of the point.  Not considered.  May be blank.
C     X1,X2,X3... Coordinates of the point
C     V1,...,VN...Optional given values and their standard deviations,
C             see input parameters KOLUMN, KOLERR and KOLFUN.
C     /...    Values must be terminated by a slash.
C (3) / or end of file.
C
C                                                     
C Input file LIN with the lines:
C (1) None to several strings terminated by / (a slash)
C (2) For each line data (2.1), (2.2) and (2.3):
C (2.1) 'NAME',X1,X2,X3,/
C     'NAME'... Name of the line.  Not considered.  May be blank but
C             must be different from '$'.
C     X1,X2,X3... Optional coordinates of the reference point of the
C             line.  Not considered.
C     /...    List of values must be terminated by a slash.
C (2.2) For each point of the line data (2.2.1):
C (2.2.1) X1,X2,X3,V1,...,VN,/
C     X1,X2,X3... Coordinates of the point of the line.
C     V1,...,VN...Optional given values and their standard deviations,
C             see input parameters KOLUMN, KOLERR and KOLFUN.
C     /...    List of values must be terminated by a slash.
C (2.3) /
C (3) / or end of file.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C Common blocks /MODELC/ and /VALC/:
      INCLUDE 'model.inc'
C     model.inc
      INCLUDE 'val.inc'
C     val.inc
C None of the storage locations of the common block are altered.
C
C External procedures directly referred:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,RARRAI,RARRAY,WARRAY,
     *  UARRAY,LENGTH,MODEL1,SRFC2,PARM2,PARM3,SOFT,VAR6,RMAT,WMAT
      INTEGER  LENGTH
      REAL UARRAY
C
C.......................................................................
C
C     Constants:
      INTEGER LU1
      PARAMETER (LU1=11)
      REAL UNDEF
C
C     Filenames:
      CHARACTER*80 FILE1
C
C     Data:
      INTEGER N1,N2,N3,M1,M2,M2IN
      INTEGER ICLASS,MPAR,KOLUMN,KOLERR,KOLFUN,INDFUN
      REAL D1,D2,D3,O1,O2,O3,POWERM,ERRMUL
C
C     Other variables:
      CHARACTER*1 TEXT
      INTEGER N1N2N3,KOLMAX,NFREE,NUNDEF,NOFF,NPTS6,NFUN
      INTEGER I,I1,I2,I3,J1,J2,J3,K1,K2,K3
      REAL F(10,47),W(47,2),POWER(47),X1,X2,X3,B0I,AUX,CS(1)
      EQUIVALENCE (F(1,1),W(1,1))
C
C     Usage of RAM:
C     RAM(1:6*NPTS):
C             RAM(1,*)=X1
C             RAM(2,*)=X2
C             RAM(3,*)=X3
C             RAM(4,*)=given value, later given minus reference value
C             IRAM(4,*)=MPAR, later RAM(4,*)=reference value
C             RAM(6,*)=standard deviation, later the variance
C     RAM(6*NPTS+1:6*NPTS+M1*M2): Derivatives
C     RAM(6*NPTS+M1*M2+1:6*NPTS+M1*M2+M1): Indices of model coefficients
C
      UNDEF=UARRAY()
C
C-----------------------------------------------------------------------
C
C     Opening data files and reading the input data:
C
C     Reading main input data:
      WRITE(*,'(A)') '+INVPTS: Enter input filename: '
      FILE1=' '
      READ (*,*) FILE1
      IF(FILE1.EQ.' ') THEN
C       INVPTS-01
        CALL ERROR('INVPTS-01: No input file specified')
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.
      END IF
      WRITE(*,'(A)') '+INVPTS: Working...            '
      CALL RSEP1(LU1,FILE1)
C
C     Reading input data for the model:
      CALL RSEP3T('MODEL' ,FILE1 ,'model.dat')
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
C     Reading input parameters describing the input values:
      CALL RSEP3I('KOLUMN',KOLUMN,0)
      CALL RSEP3I('KOLERR',KOLERR,0)
      CALL RSEP3I('KOLFUN',KOLFUN,0)
      CALL RSEP3I('INDFUN',INDFUN,0)
      KOLMAX=MAX0(3,KOLUMN,KOLERR,KOLFUN)
C
C.......................................................................
C
C     Reading the points:
C
C     Number of valid points
      NPTS6=0
C     Number of points with zero indeces (e.g., situated in free space)
      NFREE=0
C     Number of points with positive indices but undefined values
      NUNDEF=0
C     Number of defined points with positive indices outside the model
      NOFF=0
C
C     Reading gridded values:
      CALL RSEP3T('GRD',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        CALL RSEP3I('N1',N1,1)
        CALL RSEP3I('N2',N2,1)
        CALL RSEP3I('N3',N3,1)
        CALL RSEP3R('D1',D1,1.)
        CALL RSEP3R('D2',D2,1.)
        CALL RSEP3R('D3',D3,1.)
        CALL RSEP3R('O1',O1,0.)
        CALL RSEP3R('O2',O2,0.)
        CALL RSEP3R('O3',O3,0.)
        N1N2N3=N1*N2*N3
        IF(6*N1N2N3.GT.MRAM) THEN
          GO TO 99
        END IF
        CALL RARRAY(LU1,FILE1,'FORMATTED',.TRUE.,UNDEF,N1N2N3,RAM(1))
C       Reading indices of surfaces or complex blocks
        CALL RSEP3T('GRDICB',FILE1 ,' ')
        IF(FILE1.EQ.' '.AND.INDFUN.LE.0) THEN
C         INVPTS-03
          CALL ERROR('INVPTS-03: Neither GRDICB nor INDFUN specified')
C         See the description of the input data.
        END IF
        IF(FILE1.EQ.' ') THEN
          DO 11 I=N1N2N3+1,2*N1N2N3
            IRAM(I)=INDFUN
   11     CONTINUE
        ELSE
          CALL RARRAI(LU1,FILE1,'FORMATTED',.TRUE.,0,
     *                                            N1N2N3,IRAM(N1N2N3+1))
        END IF
C       Reading errors
        CALL RSEP3T('GRDERR',FILE1 ,' ')
        IF(FILE1.EQ.' ') THEN
          DO 12 I=2*N1N2N3+1,3*N1N2N3
            RAM(I)=1.
   12     CONTINUE
        ELSE
          CALL RARRAY(LU1,FILE1,'FORMATTED',.TRUE.,1.,
     *                                           N1N2N3,RAM(2*N1N2N3+1))
        END IF
C       Rearranging the values in the memory
        DO 13 I=N1N2N3+1,2*N1N2N3
          RAM (3*I-2)=RAM (I-N1N2N3)
          IRAM(3*I-1)=IRAM(I       )
          RAM (3*I  )=RAM (I+N1N2N3)
   13   CONTINUE
C       Removing undefined values and values outside the model
        J1=NINT((BOUNDM(1)-O1)/D1)
        K1=NINT((BOUNDM(2)-O1)/D1)
        J2=NINT((BOUNDM(3)-O2)/D2)
        K2=NINT((BOUNDM(4)-O2)/D2)
        J3=NINT((BOUNDM(5)-O3)/D3)
        K3=NINT((BOUNDM(6)-O3)/D3)
        IF(D1.LT.0.) THEN
          I1=J1
          J1=K1
          K1=I1
        END IF
        IF(D2.LT.0.) THEN
          I2=J2
          J2=K2
          K2=I2
        END IF
        IF(D3.LT.0.) THEN
          I3=J3
          J3=K3
          K3=I3
        END IF
        I=3*N1N2N3
        DO 23 I3=0,N3-1
          X3=O3+D3*FLOAT(I3)
          DO 22 I2=0,N2-1
            X2=O2+D2*FLOAT(I2)
            DO 21 I1=0,N1-1
              X1=O1+D1*FLOAT(I1)
              I=I+3
              IF(IRAM(I-1).LE.0) THEN
                NFREE=NFREE+1
              ELSE IF(RAM(I-2).EQ.UNDEF) THEN
                NUNDEF=NUNDEF+1
              ELSE
                IF(J3.LE.I3.AND.I3.LE.K3) THEN
                  IF(J2.LE.I2.AND.I2.LE.K2) THEN
                    IF(J1.LE.I1.AND.I1.LE.K1) THEN
                      RAM (NPTS6+1)=X1
                      RAM (NPTS6+2)=X2
                      RAM (NPTS6+3)=X3
                      RAM (NPTS6+4)=RAM (I-2)
                      IRAM(NPTS6+5)=IRAM(I-1)
                      RAM (NPTS6+6)=RAM (I  )
                      NPTS6=NPTS6+6
                    END IF
                  END IF
                END IF
              END IF
   21       CONTINUE
   22     CONTINUE
   23   CONTINUE
        NOFF=N1N2N3-NFREE-NUNDEF-NPTS6/6
      END IF
C
C     Reading individual points:
      CALL RSEP3T('PTS',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(KOLFUN.LE.0.AND.INDFUN.LE.0) THEN
C         INVPTS-04
          CALL ERROR('INVPTS-04: Neither KOLFUN nor INDFUN specified')
C         See the description of the input data.
        END IF
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        READ(LU1,*,END=39) (TEXT,I=1,20)
C       Loop over points
   31   CONTINUE
          IF(NPTS6+3+KOLMAX.GT.MRAM) THEN
            GO TO 99
          END IF
          TEXT='$'
          X1=0.
          X2=0.
          X3=0.
          READ(LU1,*,END=39) TEXT,X1,X2,X3,
     *                                 (RAM(I),I=NPTS6+7,NPTS6+3+KOLMAX)
          IF(TEXT.EQ.'$') THEN
            GO TO 39
          END IF
          NOFF=NOFF+1
          IF(BOUNDM(1).LE.X1.AND.X1.LE.BOUNDM(2)) THEN
            IF(BOUNDM(3).LE.X2.AND.X2.LE.BOUNDM(4)) THEN
              IF(BOUNDM(5).LE.X3.AND.X3.LE.BOUNDM(6)) THEN
                RAM(NPTS6+1)=X1
                RAM(NPTS6+2)=X2
                RAM(NPTS6+3)=X3
                IF(KOLUMN.GT.0) THEN
                  RAM(NPTS6+4)=RAM(NPTS6+3+KOLUMN)
                ELSE
                  RAM(NPTS6+4)=0.
                END IF
                IF(KOLFUN.GT.0) THEN
                  IRAM(NPTS6+5)=NINT(RAM(NPTS6+3+KOLFUN))
                ELSE
                  IRAM(NPTS6+5)=INDFUN
                END IF
                IF(KOLERR.GT.0) THEN
                  RAM(NPTS6+6)=RAM(NPTS6+3+KOLERR)
                ELSE
                  RAM(NPTS6+6)=1.
                END IF
                NPTS6=NPTS6+6
                NOFF=NOFF-1
              END IF
            END IF
          END IF
        GO TO 31
C       End of the loop over points
   39   CONTINUE
        CLOSE(LU1)
      END IF
C
C     Reading the lines of points:
      CALL RSEP3T('LIN',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(KOLFUN.LE.0.AND.INDFUN.LE.0) THEN
C         INVPTS-05
          CALL ERROR('INVPTS-05: Neither KOLFUN nor INDFUN specified')
C         See the description of the input data.
        END IF
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        READ(LU1,*,END=49) (TEXT,I=1,20)
C       Loop over lines
   41   CONTINUE
          TEXT='$'
          READ(LU1,*,END=90) TEXT,AUX,AUX,AUX
          IF(TEXT.EQ.'$') THEN
            GO TO 49
          END IF
C         Loop over points
   42     CONTINUE
            IF(NPTS6+3+KOLMAX.GT.MRAM) THEN
              GO TO 99
            END IF
            X1=UNDEF
            X2=0.
            X3=0.
            READ(LU1,*,END=49) TEXT,X1,X2,X3,
     *                              (RAM(I),I=NPTS6+7,NPTS6+3+KOLMAX)
            IF(X1.EQ.UNDEF) THEN
C             End of the line
              GO TO 48
            END IF
            NOFF=NOFF+1
            IF(BOUNDM(1).LE.X1.AND.X1.LE.BOUNDM(2)) THEN
              IF(BOUNDM(3).LE.X2.AND.X2.LE.BOUNDM(4)) THEN
                IF(BOUNDM(5).LE.X3.AND.X3.LE.BOUNDM(6)) THEN
                  RAM(NPTS6+1)=X1
                  RAM(NPTS6+2)=X2
                  RAM(NPTS6+3)=X3
                  IF(KOLUMN.GT.0) THEN
                    RAM(NPTS6+4)=RAM(NPTS6+3+KOLUMN)
                  ELSE
                    RAM(NPTS6+4)=0.
                  END IF
                  IF(KOLFUN.GT.0) THEN
                    IRAM(NPTS6+5)=NINT(RAM(NPTS6+3+KOLFUN))
                  ELSE
                    IRAM(NPTS6+5)=INDFUN
                  END IF
                  IF(KOLERR.GT.0) THEN
                    RAM(NPTS6+6)=RAM(NPTS6+3+KOLERR)
                  ELSE
                    RAM(NPTS6+6)=1.
                  END IF
                  NPTS6=NPTS6+6
                  NOFF=NOFF-1
                END IF
              END IF
            END IF
          GO TO 42
C         End of the loop over points
   48     CONTINUE
        GO TO 41
C       End of the loop over lines
   49   CONTINUE
        CLOSE(LU1)
      END IF
C
C.......................................................................
C
C     Matrix dimensions:
C
C     Reading number of points processed previously:
      M2IN=0
      CALL RSEP3T('M2IN',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        READ(LU1,*) M2IN
        CLOSE(LU1)
      END IF
C
C     Writing the total number of points:
      M2=M2IN+NPTS6/6
      CALL RSEP3T('M2',FILE1,'m2.out')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(I10)') M2
        CLOSE(LU1)
      END IF
C
C     Number of unknown model parameters:
      CALL RSEP3I('ICLASS',ICLASS,0)
      IF(ICLASS.LT.0.OR.2.LT.ICLASS) THEN
C       INVPTS-06
        CALL ERROR('INVPTS-06: Incorrect class index ICLASS')
C       The value of ICLASS must be 0, 1 or 2.
C       Check the input data.
      END IF
      DO 51 I=1,47
        W(I,1)=0.
        W(I,2)=0.
   51 CONTINUE
      M1=0
      IF(ICLASS.LE.1) THEN
        CALL SOFT(1,0,0,0,0,0,0,47,W,M1,IRAM(NPTS6+1),CS(1),1,CS(1))
      END IF
      IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN
        CALL SOFT(2,0,0,0,0,0,0,47,W,M1,IRAM(NPTS6+1),CS(1),1,CS(1))
      END IF
      IF(NPTS6+M1*M2+M1.GT.MRAM) THEN
        GO TO 99
      END IF
      I1=0
      IF(ICLASS.LE.1) THEN
        CALL SOFT(1,0,0,0,0,0,0,47,W,I1,IRAM(NPTS6+M1*M2+1),
     *                                                CS(1),1,CS(1))
      END IF
      IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN
        CALL SOFT(2,0,0,0,0,0,0,47,W,I1,IRAM(NPTS6+M1*M2+1),
     *                                                CS(1),1,CS(1))
      END IF
      CALL RSEP3T('M1',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(I10)') M1
        CLOSE(LU1)
      END IF
C
C     Information written to the screen:
      IF(NFREE.GT.0) THEN
        WRITE(*,'(A,I8,A)')'+INVPTS:',NFREE,' gridpoints in free space,'
        WRITE(*,'(A)') ' INVPTS: Working...'
      END IF
      IF(NUNDEF.GT.0) THEN
        WRITE(*,'(A,I8,A)') '+INVPTS:',NUNDEF,' undefined grid values,'
        WRITE(*,'(A)') ' INVPTS: Working...'
      END IF
      IF(NOFF.GT.0) THEN
        WRITE(*,'(A,I8,A)') '+INVPTS:',NOFF,' points outside model,'
        WRITE(*,'(A)') ' INVPTS: Working...'
      END IF
C
C.......................................................................
C
C     Calculating the functional values and the derivatives with respect
C     to the model coefficients at the stored points:
C
      CALL RSEP3I('MPAR'  ,MPAR  ,0)
      IF(MPAR.LT.0.OR.47.LT.MPAR) THEN
C       INVPTS-07
        CALL ERROR('INVPTS-07: Wrong material parameter')
C       Material parameters are indexed by integers from 1 to 47,
C       but a negative index or an index greater than 47 has been
C       encountered.  Check the input data.
      END IF
      CALL RSEP3R('POWERM',POWERM,1.)
      CALL RSEP3R('ERRMUL',ERRMUL,1.)
C
C     Loop over the points:
      DO 69 I3=0,NPTS6-6,6
        IF((ICLASS.LE.1.AND.MPAR.LE.0).OR.
     *     ((ICLASS.EQ.0.OR.ICLASS.EQ.2).AND.MPAR.GE.1)) THEN
C
C         Evaluation of the function:
          IF(MPAR.LE.0) THEN
C           Functions describing surfaces:
C           CALL SRFC2(IRAM(I3+5),RAM(I3+1),F)
            CALL VAL2(1,IRAM(I3+5),1,RAM(I3+1),F,POWER)
            RAM(I3+4)=RAM(I3+4)-F(1,1)
            RAM(I3+5)=F(1,1)
          ELSE
C           Material parameters:
            CALL VAL2(2,IRAM(I3+5),47,RAM(I3+1),F,POWER)
            IF(POWER(MPAR).NE.POWERM) THEN
              AUX=POWER(MPAR)/POWERM
              IF(RAM(I3+4).LE.0.) THEN
C               INVPTS-09
                CALL ERROR('INVPTS-09: Negative material parameter')
C               Different power of material parameter than that
C               interpolated by B-splines is given, and the given value
C               is not positive.  Check input data file PTS, LIN or GRD
C               for non-positive given values of the POWERMth power of
C               material parameter MPAR.
              END IF
C             Standard deviation of the power of the given value
              RAM(I3+6)=RAM(I3+6)*AUX*RAM(I3+4)**(AUX-1.)
C             Power of the given value
              RAM(I3+4)=RAM(I3+4)**AUX
            END IF
            RAM(I3+4)=RAM(I3+4)-F(1,MPAR)
            RAM(I3+5)=F(1,MPAR)
          END IF
C
C         Variance (square of the standard deviation):
          RAM(I3+6)=(RAM(I3+6)*ERRMUL)**2
C
C         Derivatives of the function with respect to model
C         coefficients:
          DO 61 I1=NPTS6+M1*(M2IN+I3/6)+1,NPTS6+M1*(M2IN+I3/6)+M1
            RAM(I1)=0.
   61     CONTINUE
          I2=0
   62     CONTINUE
            I2=I2+1
            CALL VAR6(MAX0(1,MPAR),I2,NFUN,I,B0I,AUX,AUX,AUX)
            IF(I2.LE.NFUN) THEN
              DO 63 I1=NPTS6+M1*M2+1,NPTS6+M1*M2+M1
                IF(I.EQ.IRAM(I1)) THEN
                  RAM(M1*(I3-NPTS6)/6+I1)=B0I
                  GO TO 64
                END IF
   63         CONTINUE
C             INVPTS-08
              CALL ERROR
     *                 ('INVPTS-08: Incorrect index of model parameter')
C               This error should not apperar.  Contact the author.
   64         CONTINUE
            END IF
          IF(I2.LT.NFUN) GO TO 62
C
        END IF
   69 CONTINUE
C
C.......................................................................
C
C     Writing output files:
C
C     Writing the derivatives:
      CALL RSEP3T('GM1',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF (M2IN.GT.0) THEN
          CALL RMAT(LU1,FILE1,M1,M2IN,RAM(NPTS6+1))
        END IF
        CALL WMAT(LU1,FILE1,M1,M2,RAM(NPTS6+1))
      END IF
C
C     Writing the given values minus the reference values:
      CALL RSEP3T('GM2',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RMAT(LU1,FILE1,M2IN,1,RAM(NPTS6+1))
        END IF
        DO 81 I=1,NPTS6/6
          RAM(NPTS6+M2IN+I)=RAM(6*I-2)
   81   CONTINUE
        CALL WMAT(LU1,FILE1,M2,1,RAM(NPTS6+1))
      END IF
C
C     Writing the reference values:
      CALL RSEP3T('GM3',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RMAT(LU1,FILE1,M2IN,1,RAM(NPTS6+1))
        END IF
        DO 82 I=1,NPTS6/6
          RAM(NPTS6+M2IN+I)=RAM(6*I-1)
   82   CONTINUE
        CALL WMAT(LU1,FILE1,M2,1,RAM(NPTS6+1))
      END IF
C
C     Writing the variances:
      CALL RSEP3T('DM1',FILE1 ,' ')
      IF(FILE1.NE.' ') THEN
        IF(M2IN.GT.0) THEN
          CALL RMAT(LU1,FILE1,M2IN,1,RAM(NPTS6+1))
        END IF
        DO 83 I=1,NPTS6/6
          RAM(NPTS6+M2IN+I)=RAM(6*I)
   83   CONTINUE
        CALL WMAT(LU1,FILE1,M2,1,RAM(NPTS6+1))
      END IF
C
   90 CONTINUE
      WRITE(*,'(A,I8,A)') '+INVPTS:',NPTS6/6,' points processed.'
      STOP
C
   99 CONTINUE
C     INVPTS-02
      CALL ERROR('INVPTS-02: Too small dimension MRAM of array RAM')
C     Dimension MRAM of array RAM in include file
C     ram.inc should be increased.
C     MRAM should be MAX(6*N1*N2*N3,6*NPTS+M1*M2+M1) at the least.
C     Here NPTS=M2-M2IN is the number of valid points, i.e., of points
C     situated inside the model box, with defined given values.
      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 'modelv.for'
C     modelv.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parmv.for'
C     parmv.for
      INCLUDE 'valv.for'
C     valv.for
      INCLUDE 'fitv.for'
C     fitv.for
      INCLUDE 'var.for'
C     var.for
      INCLUDE 'soft.for'
C     soft.for
      INCLUDE 'spsp.for'
C     spsp.for
C
C=======================================================================
C