C
C Program INVSOFT to evaluate the coefficients of the soft subjective
C a priori information on the perturbations of the model parameters.
C The subjective a priori information is composed of the squares of the
C Sobolev norms of the functions describing the model.
C
C Version: 6.10
C Date: 2007, June 8
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 INVSOFT 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 INVSOFT program cannot work with user's modifications of
C subroutines SRFC1, SRFC2, PARM1, and PARM2.
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 input files:
C     MODEL='string'... String containing the name of the input data
C             file specifying the model.  For description of the data
C             file refer to file 'model.for' of package MODEL.
C             Description of file MODEL
C             Default: MODEL='model.dat'
C     SOBOLEV='string'... String containing the name of the input data
C             file containing the coefficients describing the Sobolev
C             scalar product under consideration.
C             The file is not read if MODSOB=' ', see below.
C             SOBOLEV=' ' means zero Sobolev scalar products.
C             Option SOBOLEV=' ' is allowed only for simultaneous source
C             location, in order to avoid identically zero output matrix
C             MODSOB.
C             Description of file SOBOLEV
C             Default: SOBOLEV=' '
C Weighting factors of surfaces and material parameters:
C     SOBW00=real... Square root of the multiplication factor of the L2
C             and Sobolev scalar products corresponding to the functions
C             describing the surfaces.
C             Default: SOBW00=0.
C     SOBW01=real... Square root of the multiplication factor of the L2
C             and Sobolev scalar products corresponding to the
C             interpolated power of the P wave velocity.
C             Default: SOBW01=0.
C     SOBW02=real... Analogue corresponding to the S wave velocity.
C             Default: SOBW02=0.
C     SOBW03=real... Analogue corresponding to the density.
C             Default: SOBW03=0.
C     SOBW04=real... Analogue corresponding to the P wave loss factor.
C             Default: SOBW04=0.
C     SOBW05=real... Analogue corresponding to the S wave loss factor.
C             Default: SOBW05=0.
C     SOBW06=real to SOBW26=real... Analogues corresponding to the
C             reduced (i.e., divided by the density) anisotropic elastic
C             parameters A11, A12, A22, A13, A23, A33, A14, A24, A34,
C             A44, A15, A25, A35, A45, A55, A16, A26, A36, A46, A56 and
C             A66.
C             Defaults: SOBW06=0. to SOBW26=0.
C     SOBW27=real to SOBW47=real... Analogues corresponding to the
C             reduced (i.e., divided by the density) imaginary parts of
C             anisotropic elastic parameters Q11, Q12, Q22, Q13, Q23,
C             Q33, Q14, Q24, Q34, Q44, Q15, Q25, Q35, Q45, Q55, Q16,
C             Q26, Q36, Q46, Q56 and Q66.
C             Defaults: SOBW27=0. to SOBW47=0.
C Data specifying output files:
C     M1='string'... Name of the output file containing the number NM of
C             model parameters (a single integer).  The same file may be
C             generated by program 'invpts.for', and also by program
C             'invtt.for' with parameter INVSRC=' '.
C             The file is not generated if the value of M1 is blank.
C             Default: M1='m1.out'
C             Note: Default of 'invpts.for' and 'invtt.for' is M1=' '.
C     M1LOC='string'... Name of the input file containing number
C             M1LOC=NM+4*NSRC of model parameters plus the number of
C             parameters (coordinates and time) of NSRC sources.
C             This file shoud be specified if sources are to be located.
C             In such a case, the file should coincide with the output
C             file of program 'invtt.for' given by its parameter M1.
C             Symmetric NM*NM output matrix MODL2, corresponding to
C             model parameters, is then complemented to the M1LOC*M1LOC
C             symmetric matrix by zeros.
C             Symmetric NM*NM output matrix MODSOB, corresponding to
C             model parameters, is then complemented to the M1LOC*M1LOC
C             symmetric matrix by the (M1LOC-NM)*(M1LOC-NM) diagonal
C             matrix corresponding to source parameters,
C             and by zero (M1LOC-NM)*NM and NM*(M1LOC-NM) matrices.
C             Default: M1LOC=' ' means M1LOC=NM.
C     XDLOC=real... Diagonal terms of the (M1LOC-NM)*(M1LOC-NM) diagonal
C             matrix corresponding to source coordinates are set to
C             1/(XDLOC**2) for positive XDLOC.  Otherwise they are set
C             zeros.  The diagonal terms corresponding to hypocenral
C             time are set to zeros.
C     MODIND='string'... Name of the output file containing the indices
C             of model coefficients.  The indices correspond to the
C             relative location in the memory.  B-spline coefficients
C             are listed in the same order as the grid values in input
C             file MODEL.
C             The file is not generated if the value of MODIND is blank.
C             File MODIND is read by program 'modmod.for' when updating
C             the model.
C             The file has the form of an integer vector with NM
C             components.
C             Description of file MODIND
C             Default: MODIND='modind.out'
C     MODPAR='string'... Name of the matrix header file of the output
C             file containing the values of model parameters
C             (coefficients at the model basis functions).
C             The file is not generated if the value of MODPAR is blank.
C             The file has the form of a real-valued vector with NM
C             components.
C             If M1LOC is greater than NM, the vector is complemented
C             by M1LOC-NM zeros.
C             Description of file MODPAR
C             Default: MODPAR=' '
C     MODL2='string'... Name of the matrix header file of the output
C             file containing the symmetric positive-definite matrix
C             of L2 scalar products of the model basis functions.
C             The file is not generated if the value of MODL2 is blank.
C             Description of file MODL2
C             Default: MODL2=' '
C     MODSOB='string'... Name of the matrix header file of the output
C             file containing the symmetric positive-semidefinite matrix
C             of Sobolev scalar products of the model basis functions.
C             The particular kind of the Sobolev scalar product is given
C             by input file SOBOLEV.
C             The file is not generated if the value of MODSOB is blank.
C             Description of file MODSOB
C             Default: MODSOB=' '
C     Recent version of the program cannot deal with sparse matrices.
C     For general description of the files with matrices refer to file
C     forms.htm.
C Form of the output files with matrices:
C     FORMM='string' ... Form of the output files with matrices. Allowed
C             values are FORMM='formatted' and FORMM='unformatted'.
C Optional parameter specifying the form of the output formatted
C matrix data files:
C     NUMLIN=positive integer ... Number of reals or integers in one
C             line of the output file. See the description in file
C             mat.for.
C Optional data specifying the class of functions:
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
C                                                 
C Input data SOBOLEV:
C     This data file contains the coefficients describing the Sobolev
C     scalar product under consideration.
C (1) (NW1(I),NW2(I),NW3(I),I=1,NW),/
C     List of partial derivatives included in the Sobolev scalar product
C     which is assumed to represent subjective prior information about
C     the model, terminated by a slash.
C     NW1,NW2,NW3... Orders of partial derivatives with respect to
C             X1,X2,X3 coordinates.  For (bi-,tri-)cubic splines, the
C             third homogeneous partial derivatives are discontinuous.
C             NWi thus should not exceed 3, allowing for 64 different
C             partial derivatives at the most.
C (2) ((WCS(I,J),I=1,J),J=1,NW)
C     Elements of the constant symmetric weighting matrix of the Sobolev
C     scalar product.
C     WCS(I,J)... Coefficient of the product of
C             (NW1(I),NW2(I),NW3(I))-th and (NW1(J),NW2(J),NW3(J))-th
C             partial derivatives of functions in the Sobolev scalar
C             product.  The product of the derivatives is integrated
C             over the volume (surface, length) of the spline grid and
C             divided by the volume (surface, length) of the grid to
C             yield the average value of the product of the derivatives,
C             The average value is multiplied by WCS(I,J) to form the
C             contribution to the Sobolev scalar product.
C Example of data SOBOLEV
C
C                                                  
C Output file MODIND:
C (1) (INDM(I),I=1,NM),/
C     INDM... Indices of the model parameters considered by this
C             program.  The indices correspond to the relative location
C             in the memory, in array RPAR of common block /VALC/.
C             B-spline coefficients are listed in the same order as the
C             grid velocities in file MODEL.
C             Common block /VALC/
C
C                                                  
C Output file MODPAR:
C (1) (RS(I),I=1,NM)
C     RS...   Parameters (coefficients) of the input model.
C
C                                                   
C Output file MODL2:
C (1) For each column J=1,NM:
C (1.1) (BL2(I,J),I=1,J):
C     BL2...  Symmetric matrix of the L2 scalar products of the basis
C             functions corresponding to the model parameters.
C
C                                                  
C Output file MODSOB:
C (1) For each column J=1,NM:
C (1.1) (BSOB(I,J),I=1,J):
C     BSOB... Symmetric matrix of the Sobolev scalar products of the
C             basis functions corresponding to the model parameters.
C
C-----------------------------------------------------------------------
C
C Common block /VALC/:
      INCLUDE 'val.inc'
C     val.inc
C None of the storage locations of the common block are altered.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILE1,FILE2
C
C     Logical unit number:
      INTEGER LU1
      PARAMETER (LU1=1)
C
      INTEGER ICLASS,NW,NM,M1LOC
C     NW...   Number of specified partial derivatives.
C     NM...   Number of the unknown model parameters.
C     M1LOC...Number of the unknown model and source parameters.
C
C     Addresses in array RAM:
      INTEGER IWCS0,INDM0,ICS0,IB0
C     IRAM(1:3),IRAM(4:6),...,IRAM(3*NW-2:3*NW)... Orders of partial
C             derivatives.
C     IWCS0=3*NW... Origin of array WCS(I,J) of the weights describing
C             the Sobolev scalar product.
C     INDM0=IWCS0+NW*(NW+1)/2... Origin of array INDM of the indices of
C             model parameters.
C     ICS0=INDM0+NM... Origin of symmetric matrix CS of the Sobolev
C             scalar products of the basis functions corresponding to
C             the model parameters.
C     IB0=ICS0+NM*(NM+1)/2... Origin of the working array.
C
      INTEGER MW,I,J,K
      REAL SOBW(47,2),WEIGHT(47,2),XDLOC,W
C
C.......................................................................
C
C     Opening data files and reading the input data:
C
C     Main input data file read from the interactive device (*):
      WRITE(*,'(A)') '+INVSOFT: Enter input filename: '
      FILE1=' '
      READ(*,*) FILE1
      IF(FILE1.EQ.' ') THEN
C       INVSOFT-01
        CALL ERROR('INVSOFT-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)') '+INVSOFT: Working...            '
C
C     Reading main input data file:
      CALL RSEP1(LU1,FILE1)
C
C     Reading input data MODEL for the model:
      CALL RSEP3T('MODEL',FILE1,'model.dat')
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
      DO 11 I=1,47
        SOBW(I,1)=0.
        WEIGHT(I,1)=0.
        WEIGHT(I,2)=0.
   11 CONTINUE
      CALL RSEP3R('SOBW00',SOBW(01,1),0.)
      CALL RSEP3R('SOBW01',SOBW(01,2),0.)
      CALL RSEP3R('SOBW02',SOBW(02,2),0.)
      CALL RSEP3R('SOBW03',SOBW(03,2),0.)
      CALL RSEP3R('SOBW04',SOBW(04,2),0.)
      CALL RSEP3R('SOBW05',SOBW(05,2),0.)
      CALL RSEP3R('SOBW06',SOBW(06,2),0.)
      CALL RSEP3R('SOBW07',SOBW(07,2),0.)
      CALL RSEP3R('SOBW08',SOBW(08,2),0.)
      CALL RSEP3R('SOBW09',SOBW(09,2),0.)
      CALL RSEP3R('SOBW10',SOBW(10,2),0.)
      CALL RSEP3R('SOBW11',SOBW(11,2),0.)
      CALL RSEP3R('SOBW12',SOBW(12,2),0.)
      CALL RSEP3R('SOBW13',SOBW(13,2),0.)
      CALL RSEP3R('SOBW14',SOBW(14,2),0.)
      CALL RSEP3R('SOBW15',SOBW(15,2),0.)
      CALL RSEP3R('SOBW16',SOBW(16,2),0.)
      CALL RSEP3R('SOBW17',SOBW(17,2),0.)
      CALL RSEP3R('SOBW18',SOBW(18,2),0.)
      CALL RSEP3R('SOBW19',SOBW(19,2),0.)
      CALL RSEP3R('SOBW20',SOBW(20,2),0.)
      CALL RSEP3R('SOBW21',SOBW(21,2),0.)
      CALL RSEP3R('SOBW22',SOBW(22,2),0.)
      CALL RSEP3R('SOBW23',SOBW(23,2),0.)
      CALL RSEP3R('SOBW24',SOBW(24,2),0.)
      CALL RSEP3R('SOBW25',SOBW(25,2),0.)
      CALL RSEP3R('SOBW26',SOBW(26,2),0.)
      CALL RSEP3R('SOBW27',SOBW(27,2),0.)
      CALL RSEP3R('SOBW28',SOBW(28,2),0.)
      CALL RSEP3R('SOBW29',SOBW(29,2),0.)
      CALL RSEP3R('SOBW30',SOBW(30,2),0.)
      CALL RSEP3R('SOBW31',SOBW(31,2),0.)
      CALL RSEP3R('SOBW32',SOBW(32,2),0.)
      CALL RSEP3R('SOBW33',SOBW(33,2),0.)
      CALL RSEP3R('SOBW34',SOBW(34,2),0.)
      CALL RSEP3R('SOBW35',SOBW(35,2),0.)
      CALL RSEP3R('SOBW36',SOBW(36,2),0.)
      CALL RSEP3R('SOBW37',SOBW(37,2),0.)
      CALL RSEP3R('SOBW38',SOBW(38,2),0.)
      CALL RSEP3R('SOBW39',SOBW(39,2),0.)
      CALL RSEP3R('SOBW40',SOBW(40,2),0.)
      CALL RSEP3R('SOBW41',SOBW(41,2),0.)
      CALL RSEP3R('SOBW42',SOBW(42,2),0.)
      CALL RSEP3R('SOBW43',SOBW(43,2),0.)
      CALL RSEP3R('SOBW44',SOBW(44,2),0.)
      CALL RSEP3R('SOBW45',SOBW(45,2),0.)
      CALL RSEP3R('SOBW46',SOBW(46,2),0.)
      CALL RSEP3R('SOBW47',SOBW(47,2),0.)
      SOBW(1,1)=SOBW(1,1)*SOBW(1,1)
      DO 12 I=1,47
        SOBW(I,2)=SOBW(I,2)*SOBW(I,2)
   12 CONTINUE
C
C     Number and indices of unknown model parameters:
      CALL RSEP3I('ICLASS',ICLASS,0)
      IF(ICLASS.LT.0.OR.2.LT.ICLASS) THEN
C       INVSOFT-02
        CALL ERROR('INVSOFT-02: Incorrect class index ICLASS')
C       The value of ICLASS must be 0, 1 or 2.
C       Check the input data.
      END IF
      INDM0=0
      NM=0
      IF(ICLASS.LE.1) THEN
        CALL SOFT(1,0,0,0,0,0,0,47,WEIGHT,NM,IRAM(INDM0+1),RAM,1,RAM)
      END IF
      IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN
        CALL SOFT(2,0,0,0,0,0,0,47,WEIGHT,NM,IRAM(INDM0+1),RAM,1,RAM)
      END IF
C     WRITE(*,'(A,I5,A)') '+INVSOFT:',NM,' model parameters'
C     (Last two actual arguments RAM are not used since WEIGHT(*,*)=0.)
C     (We have just hoped here that array IRAM is sufficiently large.)
C
C     Reading input file M1LOC:
      M1LOC=NM
      CALL RSEP3T('M1LOC',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        I=MAX0(INDEX(FILE1,'          ')-1,11)
        WRITE(*,'(2A)') '+INVSOFT: Reading file ',FILE1(1:I)
        OPEN(LU1,FILE=FILE1)
        READ(LU1,*) M1LOC
        CLOSE(LU1)
        IF(MOD(M1LOC-NM,4).NE.0) THEN
C         INVSOFT-03
          CALL ERROR('INVSOFT-03: Incorrect number M1LOC of parameters')
C         Number M1LOC-NM of source parameters must be multiple of 4,
C         because each source has 4 parameters: 3 coordinates and the
C         hypocentral time.
        END IF
      END IF
C
C     Writing output file M1:
      CALL RSEP3T('M1',FILE1,'m1.out')
      IF(FILE1.NE.' ') THEN
        I=MAX0(INDEX(FILE1,'          ')-1,11)
        WRITE(*,'(2A)') '+INVSOFT: Writing file ',FILE1(1:I)
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(I10)') NM
        CLOSE(LU1)
      END IF
C
C     Writing output file MODIND:
      CALL RSEP3T('MODIND',FILE1,'modind.out')
      IF(FILE1.NE.' ') THEN
        I=MAX0(INDEX(FILE1,'          ')-1,11)
        WRITE(*,'(2A)') '+INVSOFT: Writing file ',FILE1(1:I)
        OPEN(LU1,FILE=FILE1)
        WRITE(LU1,'(10(I7,1X))') (IRAM(I),I=INDM0+1,INDM0+NM)
        CLOSE(LU1)
      END IF
C
C     Writing output file MODPAR:
      CALL RSEP3T('MODPAR',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        IF(INDM0+NM+M1LOC.GE.MRAM) THEN
C         INVSOFT-04
          CALL ERROR('INVSOFT-04: Too small array RAM')
C         Dimension MRAM of array RAM in include file
C         ram.inc
C         should be increased to accommodate the input coefficients of
C         the Sobolev scalar product, the indices of model parameters,
C         and if MODPAR is not blank, also the values of model
C         parameters.
C         Consider simultaneously the memory requirements described at
C         errors INVSOFT-05,
C         INVSOFT-06,
C         INVSOFT-10,
C         and INVSOFT-11.
        END IF
        DO 13 I=INDM0+1,INDM0+NM
          RAM(I+NM)=RPAR(IRAM(I))
   13   CONTINUE
        DO 14 I=INDM0+NM+1,INDM0+M1LOC
          RAM(I+NM)=0.
   14   CONTINUE
        CALL WMAT(LU1,FILE1,M1LOC,1,RAM(INDM0+NM+1))
      END IF
C
C     Allocating memory for L2 scalar products:
      ICS0=INDM0+NM
      IB0=ICS0+NM*(NM+1)/2
      IF(IB0.GE.MRAM) THEN
C       INVSOFT-05
        CALL ERROR('INVSOFT-05: Too small array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc
C       should be increased to accommodate the indices of model
C       parameters, the output symmetric matrix CS of subroutine SOFT
C       containing the L2 scalar products of the basis functions
C       corresponding to the model parameters, and working array B of
C       subroutine SOFT.
C       Consider simultaneously the memory requirement described at
C       errors INVSOFT-06,
C       INVSOFT-10,
C       and INVSOFT-11.
      END IF
      IF(ICS0+M1LOC*(M1LOC+1)/2.GE.MRAM) THEN
C       INVSOFT-06
        CALL ERROR('INVSOFT-06: Too small array RAM')
C       Dimension MRAM of array RAM in include file
C       ram.inc
C       should be increased to accommodate NM indices of model
C       parameters and the output M1LOC*M1LOC symmetric matrix of the
C       Sobolev scalar products.
      END IF
C
C     Calculating and writing the L2 scalar products:
      CALL RSEP3T('MODL2',FILE1,' ')
      IF(FILE1.NE.' ') THEN
        WRITE(*,'(A)') '+INVSOFT: Calculating L2 scalar products'
        DO 15 I=ICS0+1,IB0
          RAM(I)=0.
   15   CONTINUE
        NM=0
        IF(ICLASS.LE.1) THEN
          CALL SOFT(1,0,0,0,0,0,0,47,SOBW,
     *              NM,IRAM(INDM0+1),RAM(ICS0+1),MRAM-IB0,RAM(IB0+1))
        END IF
        IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN
          CALL SOFT(2,0,0,0,0,0,0,47,SOBW,
     *              NM,IRAM(INDM0+1),RAM(ICS0+1),MRAM-IB0,RAM(IB0+1))
        END IF
        DO 16 I=IB0+1,ICS0+M1LOC*(M1LOC+1)/2
          RAM(I)=0.
   16   CONTINUE
        CALL WMAT(LU1,FILE1,M1LOC,0,RAM(ICS0+1))
      END IF
C
C.......................................................................
C
      CALL RSEP3T('MODSOB',FILE1,' ')
      IF(FILE1.NE.' ') THEN
C
C       Input data SOBOLEV:
        CALL RSEP3T('SOBOLEV',FILE2,' ')
        IF(FILE2.EQ.' ') THEN
          IF(M1LOC.EQ.NM) THEN
C           INVSOFT-07
            CALL ERROR('INVSOFT-07: No input file SOBOLEV specified')
C           If parameter MODSOB is not blank and no simultaneous
C           source location is performed, parameter SOBOLEV must be
C           specified and not blank.
C           See the input data.
          END IF
          NW=0
          INDM0=0
        ELSE
          I=MAX0(INDEX(FILE2,'          ')-1,11)
          WRITE(*,'(2A)') '+INVSOFT: Reading file ',FILE2(1:I)
          OPEN(LU1,FILE=FILE2,STATUS='OLD')
C         Reading prior subjective information coefficients:
C         Maximum number MW of different partial derivatives
          MW=MIN0(64,(MRAM-1)/3)
          DO 21 I=1,3*MW+1
            IRAM(I)=-1
   21     CONTINUE
          READ(LU1,*) (IRAM(I),I=1,3*MW+1)
          DO 22 I=1,3*MW+1
            IF(IRAM(I).LT.0) THEN
              NW=(I-1)/3
              IF(3*NW.NE.I-1) THEN
C             INVSOFT-08
                CALL ERROR('INVSOFT-08: Wrong partial derivatives')
C               The input partial derivatives do not form triplets,
C               or some of the derivatives is of a negative order.
              END IF
              GO TO 23
            END IF
   22     CONTINUE
C           INVSOFT-09
            CALL ERROR('INVSOFT-09: Too many partial derivatives')
C           The number of input triplets of partial derivatives is
C           greater than the maximum number MW defined few lines above.
   23     CONTINUE
          IWCS0=3*NW
          INDM0=IWCS0+NW*(NW+1)/2
          IF(IB0.GE.MRAM) THEN
C           INVSOFT-10
            CALL ERROR('INVSOFT-10: Too small array RAM')
C           Dimension MRAM of array RAM in include file
C           ram.inc
C           should be increased to accommodate the input coefficients of
C           the Sobolev scalar product, the indices of model parameters,
C           the output symmetric matrix CS of subroutine SOFT containing
C           the Sobolev scalar products of the basis functions
C           corresponding to the model parameters, and working array B
C           of subroutine SOFT.
          END IF
          READ(LU1,*) (RAM(I),I=IWCS0+1,INDM0)
          CLOSE(LU1)
        END IF
        ICS0=INDM0+NM
        IB0=ICS0+NM*(NM+1)/2
        IF(ICS0+M1LOC*(M1LOC+1)/2.GE.MRAM) THEN
C         INVSOFT-11
          CALL ERROR('INVSOFT-11: Too small array RAM')
C         Dimension MRAM of array RAM in include file
C         ram.inc
C         should be increased to accommodate the input coefficients of
C         the Sobolev scalar product, the indices of model parameters,
C         and the output M1LOC*M1LOC symmetric matrix of the Sobolev
C         scalar products.
        END IF
C
C       Generating the Sobolev scalar products:
        WRITE(*,'(A)') '+INVSOFT: Calculating Sobolev scalar products'
        DO 51 I=ICS0+1,IB0
          RAM(I)=0.
   51   CONTINUE
        DO 59 J=1,NW
          DO 58 I=1,J
            W=RAM(IWCS0+J*(J-1)/2+I)
            IF(W.NE.0.) THEN
              WEIGHT(1,1)=SOBW(1,1)*W
              DO 52 K=1,47
                WEIGHT(K,2)=SOBW(K,2)*W
   52         CONTINUE
              NM=0
              IF(ICLASS.LE.1) THEN
                CALL SOFT(1,IRAM(3*I-2),IRAM(3*I-1),IRAM(3*I),
     *                      IRAM(3*J-2),IRAM(3*J-1),IRAM(3*J),47,WEIGHT,
     *                 NM,IRAM(INDM0+1),RAM(ICS0+1),MRAM-IB0,RAM(IB0+1))
              END IF
              IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN
                CALL SOFT(2,IRAM(3*I-2),IRAM(3*I-1),IRAM(3*I),
     *                      IRAM(3*J-2),IRAM(3*J-1),IRAM(3*J),47,WEIGHT,
     *                 NM,IRAM(INDM0+1),RAM(ICS0+1),MRAM-IB0,RAM(IB0+1))
              END IF
            END IF
   58     CONTINUE
   59   CONTINUE
C
C       Writing output file MODSOB:
        DO 61 I=IB0+1,ICS0+M1LOC*(M1LOC+1)/2
          RAM(I)=0.
   61   CONTINUE
        CALL RSEP3R('XDLOC',XDLOC,0.)
        IF(XDLOC.GT.0.) THEN
          W=1./XDLOC**2
          DO 63 J=NM+1,M1LOC,4
            DO 62 I=J,J+2
              RAM(ICS0+I*(I+1)/2)=W
   62       CONTINUE
   63     CONTINUE
        END IF
        CALL WMAT(LU1,FILE1,M1LOC,0,RAM(ICS0+1))
      END IF
C
      WRITE(*,'(A)') '+INVSOFT: Done.                              '
      STOP
      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 'mat.for'
C     mat.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 'spsp.for'
C     spsp.for
      INCLUDE 'soft.for'
C     soft.for
      INCLUDE 'indexi.for'
C     indexi.for
C
C=======================================================================
C