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: 5.40
C Date: 2000, May 29
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.  Otherwise,
C             the filename should be specified.
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 programs 'invpts.for' and 'invtt.for'.
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     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 integer vector of NM components.
C             Description of file MODIND
C             Default: MODIND='modind.out'
C     MODPAR='string'... Name of the output file containing the values
C             of model parameters (coefficients at the model basis
C             functions).
C             The file is not generated if the value of MODPAR is blank.
C             The file has the form real-valued vector of NM components.
C             Description of file MODPAR
C             Default: MODPAR=' '
C     MODL2='string'... Name of the output file containing the symmetric
C             positive-definite matrix of L2 scalar products of the
C             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 output file containing the
C             symmetric positive-semidefinite matrix of Sobolev scalar
C             products of the model basis functions.  The particular
C             kind of the Sobolev scalar product is given by input file
C             SOBOLEV.
C             The file is not generated if the value of MODSOB is blank.
C             Description of file MODSOB
C             Default: MODSOB=' '
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 NW,NM
C     NW...   Number of specified partial derivatives.
C     NM...   Number of the unknown model 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
C
      INTEGER MW,I,J,K
      REAL SOBW(47,2),WEIGHT(47,2),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:
      INDM0=0
      CALL SOFT(2,0,0,0,0,0,0,47,WEIGHT,NM,IRAM(INDM0+1),RAM,1,RAM)
C     (We have just hoped here that array IRAM is sufficiently large.)
C     WRITE(*,'(A,I5,A)') '+INVSOFT:',NM,' model parameters'
      ICS0=INDM0+NM
      IB0=ICS0+NM*(NM+1)/2
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
        CALL WMAT(LU1,FILE1,NM,1,RAM(INDM0+1))
      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 13 I=ICS0+1,IB0
          RAM(I)=0.
   13   CONTINUE
        CALL SOFT(2,0,0,0,0,0,0,47,SOBW,
     *            NM,IRAM(INDM0+1),RAM(ICS0+1),MRAM-IB0,RAM(IB0+1))
        CALL WMAT(LU1,FILE1,NM,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
C         INVSOFT-02
          CALL ERROR('INVSOFT-02: No input file SOBOLEV specified')
C         If parameter MODSOB is not blank, parameter SOBOLEV must be
C         specified and not blank.  There is no default filename.
C         See the input data.
        END IF
        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-03
              CALL ERROR('INVSOFT-03: 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-04
          CALL ERROR('INVSOFT-04: Too many partial derivatives')
C         The number of input triplets of partial derivatives is greater
C         than the maximum number MW defined few lines above.
   23   CONTINUE
        IWCS0=3*NW
        INDM0=IWCS0+NW*(NW+1)/2
        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 input coefficients of
C         the Sobolev scalar product and the output symmetric matrix CS
C         of the Sobolev scalar products of the basis functions
C         corresponding to the model parameters.
        END IF
        READ(LU1,*) (RAM(I),I=IWCS0+1,INDM0)
        CLOSE(LU1)
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
              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
   58     CONTINUE
   59   CONTINUE
C
C       Writing output file MODSOB:
        CALL WMAT(LU1,FILE1,NM,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 '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
C
C=======================================================================
C