C Program INV1 to evaluate the derivatives of the travel time with
C respect to the model coefficients.
C
C Just a preliminary demo version, illustrating the usage of the
C routines designed to calculate the variations of travel times with
C respect to model parameters (FORTRAN77 source code files 'var.for',
C 'spsp.for', 'ap.for').
C
C Program INV1 assumes all model parameters (coefficients) stored in the
C common block /VALC/ as in the submitted versions of user-defined model
C specification FORTRAN77 source code files 'srfc.for', 'parm.for', and
C 'val.for'.  Thus, unlike the other parts of the complete ray tracing,
C the INV1 program cannot work with user's modifications of the
C subroutines SRFC1, SRFC2, PARM1, and PARM2.
C
C Main input data file read from the interactive device (*):
C (1) 'INV1'
C     'INV1' name of the input file described below.
C     Default: 'INV1'='inv1.out'.
C
C Input data 'INV1':
C     This main data file contains the names of other input and output
C     files.  The data are read in by the list directed input (free
C     format).  In the list of input data below, each numbered paragraph
C     indicates the beginning of a new input operation (new READ
C     statement).  All input variables are of the type CHARACTER.  Only
C     the first 80 characters of the strings are significant.
C (1) 'POINTS','FTT','MODEL','DATA','SOFT','INV1LOG',/
C     'POINTS'... Input data file containing the coordinates of shot and
C             receiver points.  Its structure is described below.
C             Ignored (not opened) if the 'data' filename below is
C             blank.
C     'FTT'... Input data file containing the travel times from the
C             field measurements.  Its structure is described below.
C             Ignored (not opened) if the 'DATA' filename below is
C             blank.
C     'MODEL'... Input data file containing the model parameters.
C             See the file 'model.for' of the package 'MODEL'.
C     'DATA'..Output file containing the computed travel times and their
C             derivatives with respect to the model coefficients.
C             Its structure is described below.
C             If the filename is blank, no file with objective prior
C             information is generated.
C     'SOFT'..Output file containing the subjective prior information.
C             Its structure is described below.
C             If the filename is blank, no file with subjective prior
C             information is generated.
C     'INV1LOG'... Output log file.  Just for your information.
C             Not opened and generated if the 'DATA' filename above is
C             blank.
C     /...    An obligatory slash at the end of line.
C     Default: 'MODEL'='model.dat', 'DATA'='data.out',
C             'SOFT'='soft.out', 'INV1LOG'='inv1log.out'.
C (2) DIST,VPOWER,ERRINT,ERRVEL,/
C     DIST... Maximum distance between the source or receiver point and
C             the initial or end point of a synthetic ray.
C     VPOWER... Velocity power for the computation of the travel-time
C             check sum.  If the VPOWER-th velocity power is expressed,
C             in all blocks of the model, in terms of explicit functions
C             of model coordinates, linearly homogeneous in their
C             coefficients (e.g. B-splines), the travel time minus its
C             check sum (see the log output file) should be zero within
C             rounding errors.  Otherwise, the check sum may have no
C             sense.
C             VPOWER=0: no check sum is evaluated and printed.
C     ERRINT..Weighting factor of the posterior interface B-spline
C             coefficient error.
C     ERRVEL..Weighting factor of the posterior velocity B-spline
C             coefficient error.
C     /...    An obligatory slash at the end of line.
C     Default: VPOWER=0.0, ERRINT=1.0, ERRVEL=1.0.
C (3) Any times the following data:
C     'RAYS','ENDPOINTS','INITIALPOINTS'
C     'RAYS'... File with the quantities stored along rays (see
C             C.R.T.5.5.1).
C     'ENDPOINTS'... File with the quantities at the points of
C             intersection of rays with the specified surface at which
C             the receivers are situated for the case of two-point ray
C             tracing (see C.R.T.5.5.2).
C             If this filename is not blank, just the two-point rays
C             with minimum travel time at each receiver are considered.
C             If this filename is blank, all rays are taken into
C             account.
C             Attention:  All rays taken into account must start in some
C             of the specified sources and terminate in some of the
C             specified receivers, see the input file 'FTT'.
C     'INITIALPOINTS'... File with the quantities at the initial points
C             of rays, corresponding to the above file rays or endpoints
C             (see C.R.T.6.1).
C (4) / (a slash).
C (5) (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.
C     NW1,NW2,NW3... Orders of partial derivatives with respect to
C             X1,X2,X3 coordinates.
C (6) ((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.
C
C Input data file 'POINTS':
C     This data file contains the coordinates of shot and receiver
C     points.  The data are read in by the list directed input
C     (free format).  In the list of input data below, each numbered
C     paragraph indicates the beginning of a new input operation (new
C     READ statement).  The CHARACTER strings are explicitly mentioned
C     in this description.  Otherwise, if the first letter of the
C     symbolic name of the input variable is I-N, the corresponding
C     value in input data must be of the type INTEGER.  Otherwise, the
C     input parameter is of the type REAL.
C (1) Several strings terminated by / (a slash).
C (2) List of the sources and receivers: Any times the following data:
C (2.1) POINT,X1,X2,X3
C     POINT...CHARACTER*11 string containing the name of the source or
C             receiver point.
C     X1,X2,X3... Coordinates of the source or receiver point.
C (3) / (a slash) or the end of file.
C
C Input data file 'FTT':
C     This data file contains the travel times from the field
C     measurements.  The data are read in by the list directed input
C     (free format).  In the list of input data below, each numbered
C     paragraph indicates the beginning of a new input operation (new
C     READ statement).  The CHARACTER strings are explicitly mentioned
C     in this description.  Otherwise, if the first letter of the
C     symbolic name of the input variable is I-N, the corresponding
C     value in input data must be of the type INTEGER.  Otherwise, the
C     input parameter is of the type REAL.
C (1) Several strings terminated by / (a slash).
C (2) List of the travel times: Any times the following data (2.1):
C (2.1) SOURCE,RECEIVER,TFIELD,TERR
C     SOURCE..CHARACTER*11 string containing the name of the source
C             point.
C     RECEIVER... CHARACTER*11 string containing the name of the
C             receiver point.  The source and receiver points may be
C             mutually interchanged.
C     TFIELD..Travel time from a field measurement.
C     TERR... Error of the travel time from a field measurement.
C (3) / (a slash) or the end of file.
C
C Output file 'DATA':
C (1) ND,NM,(INDM(I),I=1,NM)
C     ND...   The number of data (i.e. the number of equations).
C     NM...   Number of unknown model coefficients.
C     INDM... Indices of the model coefficients influencing the travel
C             times.  The indices correspond to the relative location in
C             the memory.
C             B-spline coefficients are listed in the same order as the
C             grid velocities in the file 'MODEL'.
C (2) ND-times the following data (2.1):
C (2.1) KD,RD,ED,WD,(GD(I),I=1,NM)
C     KD...   Index of the field travel time within the file 'ftt'.
C             The field travel times are indexed consecutively 1,2,3,...
C     RD...   Field travel time minus the computed synthetic travel
C             time.  In the case of multiple two-point rays, the first
C             arrival of them is considered.
C     ED...   Prior error of the above travel time difference.
C             Identical to TERR, file 'FTT', part (3).  I.e. the
C             synthetic travel time is assumed to be sufficiently
C             accurate.
C     WD...   Field travel time (stored for the purpose to assess a
C             posterior data misfit covariance weighting matrix).
C     GD...   Derivatives of the synthetic travel time with respect to
C             the model coefficients INDM.
C
C Output file 'SOFT':
C (1) (INDM(I),I=1,NM)
C     INDM... Indices of the model coefficients influencing the travel
C             times.
C (2) (WM(I),I=1,NM)
C     RS...   Parameters (coefficients) of the initial model.
C (3) (WM(I),I=1,NM)
C     WM...   Relative posterior errors of the model parameters.
C             WM=ERRVEL  for material parameters (see (2) of the input
C               data 'INV1'),
C             WM=ERRINT  for structural interfaces (see (2) of the input
C               data 'INV1').
C     Default: ERRVEL=1.0, ERRINT=1.0.
C (4) ((CS(I,J),I=1,J),J=1,NM):
C     CS...   Relative inverse subjective prior information covariance
C             matrix.  Note that subjective data are assumed to be minus
C             the above initial model parameters, and the matrix GS
C             projecting model parameters onto the subjective data is
C             assumed to be identity.
C
C Output log file 'LOG':
C (1) For each considered ray:
C (1.1) SOURCE,RECEIVER,TFIELD,TDIF,SDIST,RDIST,CHECK
C     SOURCE..Name of the source point.
C     RECEIVER... Name of the receiver point.
C     TFIELD..Travel time from a field measurement.
C     TDIF... Field travel time minus the minimum synthetic travel time.
C     SDIST...Distance between the source and the initial point of the
C             synthetic ray.
C     RDIST...Distance between the receiver and the end point of the
C             synthetic ray.
C     CHECK...Synthetic travel time minus the travel time resulting from
C             the derivatives of the theoretical travel time with
C             respect to the model coefficients.  This quantity should
C             not exceed in order the numerical error of the synthetic
C             travel time.
C             In this version defined just for the models described in
C             the terms of velocity.
C
C Date: 1994, January 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Common block /VALC/:
      INCLUDE 'val.inc'
C None of the storage locations of the common block are altered.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4,FILE5,FILE6
C
C     Logical unit numbers:
      INTEGER LU0,LU1,LU2,LU3,LU4,LU5,LU7,IU1,IU2
      PARAMETER (LU0=10)
      PARAMETER (LU1=11)
      PARAMETER (LU2=12)
      PARAMETER (LU3=13)
      PARAMETER (LU4=14)
      PARAMETER (LU5=15)
      PARAMETER (LU7=17)
C
C     Input data:
      INTEGER MP,MT
      PARAMETER (MP=1000)
      PARAMETER (MT=4000)
      CHARACTER*1  TEXT(10)
      CHARACTER*11 POINT(MP),SRC,REC
      INTEGER NP,NT,KS(MT),KR(MT)
      REAL DIST,DIST2,VPOWER,ERRINT,ERRVEL
      REAL COOR(3,MP),TFIELD(MT),TERR(MT)
      INTEGER MW
      PARAMETER (MW=12)
      INTEGER NW,NW1(MW+1),NW2(MW),NW3(MW)
      REAL WCS(MW*(MW+1)/2)
C     POINT...Names of the source and receiver points.
C     NP...   Number of source and receiver points.
C     NT...   Number of field travel times.
C     KS(I)...Index of the source point corresponding to the I-th field
C             travel time.
C     KR(I)...Index of the receiver point corresponding to the I-th
C             field travel time.
C     DIST... Maximum distance between the source or receiver point and
C             the initial or end point of a synthetic ray.
C     DIST2...DIST**2
C     VPOWER... Velocity power for the computation of the travel-time
C             check SUM.
C     COOR... Coordinates of the source or receiver points.
C     TFIELD..Field travel times.
C     TERR... Field travel time errors.
C     NW...   Number of specified partial derivatives.
C     NW1,NW2,NW3... Orders of partial derivatives.
C     WCS(I,J)... Coefficient matrix of the Sobolev scalar product.
C
C     Output data:  variations of the synthetic travel time:
      INTEGER NSUM,NM,INDM(NPAR)
      REAL SUM(NPAR)
C     NM...   Number of the unknown model parameters.
C     INDM... Indices of the unknown model parameters.
C
C     Output data:  subjective prior information:
      INTEGER MM,MCS
C     MM... Maximum number of model parameters.
      PARAMETER (MM=200)
      PARAMETER (MCS=MM*(MM+1)/2)
      REAL CS(MCS)
C
C     Auxiliary storage locations:
      INTEGER IS,IR,IT,ND,IRAYTT,I,J,K
      REAL TTMIN(MT)
      REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6)
      REAL SI,SE,RI,RE,SDIST,RDIST,WEIGHT
C     IS...   Index of the source point.
C     IR...   Index of the receiver point.
C     IT...   Index of the field travel time.
C     ND...   Number of synthetic travel times corresponding to the
C             field travel times.
C     IRAYTT..Index of the last processed ray.
C     I,J,K...Temporary storage locations.
C     TTMIN...Minimum synthetic travel times corresponding to the
C             individual field travel times.
C     TT...   Synthetic travel time.
C     TTAUX...Temporary storage location.
C     TDIF... Field travel time minus the minimum synthetic travel time.
C     XI1,XI2,XI3,XE1,XE2,XE3... Coordinates of the initial and end
C             points of the last processed ray.
C     PI,PE...Slowness vectors at the initial and end points of the
C             last processed ray.
C     SI,SE,RI,RE... Squares of the distances between the source or
C             receiver points and the initial or end points of the ray.
C     SDIST...Distance between the source and the initial point of the
C             synthetic ray.
C     RDIST...Distance between the receiver and the end point of the
C             synthetic ray.
C     WEIGHT..Weighting factor corresponding to the derivatives
C             (average quadratic value of the derivatives on input)
C
C     More than 23*MP+16*MT+4*MM*(MM+1) Bytes of memory required.
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)') ' Enter the name of the main input data file: '
      FILE0='inv1.dat'
      READ(*,*) FILE0
      WRITE(*,'(A)') '+                                            '
C
C     Input data 'INV1':
      OPEN(LU0,FILE=FILE0,STATUS='OLD')
      FILE2='model.dat'
      FILE4='data.out'
      FILE5='soft.out'
      FILE6='inv1log.out'
      READ(LU0,*) FILE0,FILE1,FILE2,FILE4,FILE5,FILE6
      VPOWER=0.
      ERRINT=1.
      ERRVEL=1.
      READ(LU0,*) DIST,VPOWER,ERRINT,ERRVEL
      DIST2=DIST*DIST
      OPEN(LU4,FILE=FILE2,STATUS='OLD')
      CALL MODEL1(LU4)
      CLOSE(LU4)
C
C     Number of unknown model parameters:
      CALL SOFT(2,0,0,0,0,0,0,0.,NM,INDM,CS)
      WRITE(*,'(A,I4,A)') '+',NM,' model parameters'
C
      IF(FILE4.EQ.' ') THEN
C       No objective information file generated.
   10   CONTINUE
          FILE1=' '
          READ(LU0,*,END=70) FILE1,FILE2,FILE3
          IF(FILE1.EQ.' ') THEN
C           GO TO (subjective prior information)
            GO TO 80
          END IF
        GO TO 10
      END IF
C
C.......................................................................
C
C     Reading source and receiver points:
C
      OPEN(LU4,FILE=FILE0,STATUS='OLD')
      NP=-1
      READ(LU4,*,END=2) TEXT
    1 CONTINUE
        NP=NP+1
        IF(NP.GE.MP) THEN
          PAUSE 'Error: Too many source and receiver points'
          STOP
        END IF
        POINT(NP+1)=' '
        READ(LU4,*,END=2) POINT(NP+1),(COOR(I,NP+1),I=1,3)
      IF(POINT(NP+1).NE.' ') GO TO 1
    2 CONTINUE
      CLOSE(LU4)
C
C     Reading field travel times:
C
      OPEN(LU4,FILE=FILE1,STATUS='OLD')
      NT=0
      READ(LU4,*,END=2) TEXT
    3 CONTINUE
        NT=NT+1
        IF(NT.GT.MT) THEN
          PAUSE 'Error: Too many field travel times'
          STOP
        END IF
        SRC=' '
        READ(LU4,*,END=9) SRC,REC,TFIELD(NT),TERR(NT)
        IF(SRC.EQ.' ') THEN
          GO TO 9
        END IF
        DO 4 I=1,NP
          IF(SRC.EQ.POINT(I)) THEN
            KS(NT)=I
            GO TO 5
          END IF
    4   CONTINUE
        PAUSE 'Error: Source name not recognized'
        STOP
    5   CONTINUE
        DO 6 I=1,NP
          IF(REC.EQ.POINT(I)) THEN
            KR(NT)=I
            GO TO 7
          END IF
    6   CONTINUE
        PAUSE 'Error: Receiver name not recognized'
        STOP
    7   CONTINUE
      GO TO 3
    9 CONTINUE
      NT=NT-1
      CLOSE(LU4)
C
C.......................................................................
C
C     Computing quantities describing objective prior information:
C
      OPEN(LU5,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN(LU7,FILE=FILE6)
      WRITE(*,*)
C
      KS(NT+1)=NP+1
      KR(NT+1)=NP+1
      TFIELD(NT+1)=0.
      IRAY=0
      IWAVE=0
      NSUM=IPAR(IPAR(IPAR(2)))
      DO 12 I=1,NT
        TTMIN(I)=999999.
   12 CONTINUE
C
C     Loop for the files with computed rays
   20 CONTINUE
        FILE1=' '
        READ(LU0,*,END=70) FILE1,FILE2,FILE3
        IF(FILE1.EQ.' ') THEN
          GO TO 70
        END IF
        I=INDEX(FILE1,'          ')
        J=INDEX(FILE2,'          ')
        K=INDEX(FILE3,'          ')
        WRITE(*,'(''+Processing: '',3A)') FILE1(1:I),FILE2(1:J),
     *                                                        FILE3(1:K)
        WRITE(*,*)
        OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
        IF(FILE2.EQ.' ') THEN
          IU1=0
          IU2=LU1
        ELSE
          IU1=LU1
          IU2=LU2
          OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
        END IF
        OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD')
        IRAYTT=0
C
C       Loop for the endpoints of rays
   30   CONTINUE
C         Reading the results of the complete ray tracing
          CALL AP00(IU1,IU2,LU3)
          IF(IPT.LE.1.AND.IRAYTT.NE.0)THEN
C           New ray - recording results for the last ray iraytt
C           loop for field travel times - searching for two-point ray
            DO 39 IT=1,NT
              IS=KS(IT)
              IR=KR(IT)
              SI=(COOR(1,IS)-XI1)**2+(COOR(2,IS)-XI2)**2
     *                              +(COOR(3,IS)-XI3)**2
              RI=(COOR(1,IR)-XI1)**2+(COOR(2,IR)-XI2)**2
     *                              +(COOR(3,IR)-XI3)**2
              SE=(COOR(1,IS)-XE1)**2+(COOR(2,IS)-XE2)**2
     *                              +(COOR(3,IS)-XE3)**2
              RE=(COOR(1,IR)-XE1)**2+(COOR(2,IR)-XE2)**2
     *                              +(COOR(3,IR)-XE3)**2
              IF(SE.LE.SI.AND.RI.LE.RE) THEN
C               interchanging source and receiver points
                SI=SE
                RE=RI
                IS=KR(IT)
                IR=KS(IT)
              END IF
              IF((SI.LE.DIST2.AND.RE.LE.DIST2)) THEN
C               Synthetic ray may correspond to the field travel time
C               check for ray directions near the source
C               (allowable angle deviation +-30 deg: cosine**2=0.750)
C               (allowable angle deviation +-15 deg: cosine**2=0.933)
                TTAUX=(COOR(1,IR)-COOR(1,IS))*(XE1-XI1)
     *               +(COOR(2,IR)-COOR(2,IS))*(XE2-XI2)
     *               +(COOR(3,IR)-COOR(3,IS))*(XE3-XI3)
                IF(TTAUX.GT.0..AND.TTAUX*TTAUX.GT.
     *                         0.933*((COOR(1,IR)-COOR(1,IS))**2
     *                               +(COOR(2,IR)-COOR(2,IS))**2
     *                               +(COOR(3,IR)-COOR(3,IS))**2)
     *                  *((XE1-XI1)**2+(XE2-XI2)**2+(XE3-XI3)**2) ) THEN
C                 Synthetic ray corresponds to the field travel time
                  TTAUX=TT-PI(1)*(COOR(1,IS)-XI1)
     *                    -PI(2)*(COOR(2,IS)-XI2)
     *                    -PI(3)*(COOR(3,IS)-XI3)
     *                    +PE(1)*(COOR(1,IR)-XE1)
     *                    +PE(2)*(COOR(2,IR)-XE2)
     *                    +PE(3)*(COOR(3,IR)-XE3)
                  IF(TTAUX.LT.TTMIN(IT)) THEN
                    TTMIN(IT)=TTAUX
C                   Possible minimum synthetic travel time
                    SDIST=SQRT(SI)
                    RDIST=SQRT(RE)
                    WRITE(LU5) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM)
                  END IF
                END IF
              END IF
   39       CONTINUE
            IF(NT.EQ.0) THEN
              WRITE(LU5) NT+1,TT,TT,0.,0.,(SUM(I),I=1,NSUM)
            END IF
            IRAYTT=0
          END IF
          IF(IWAVE.EQ.0) THEN
            GO TO 60
          END IF
C         *** for future extensions (selection of two-point rays):
C         IF(IU1.NE.0) THEN
C           CALL AP30(IREC)
C           IF(IREC.EQ.0) THEN
C             IF(IPT.LE.1.) THEN
C               WRITE(*,'(''+WAVE:'',I3,''    RAY:'',I4,''    POINT:'',
C    *                                              I4)') IWAVE,IRAY,IPT
C             END IF
C             GO TO 30
C           END IF
C         END IF
C         ***
          IF(IPT.EQ.1.OR.MOD(IPT,10).EQ.0) THEN
            WRITE(*,'(''+WAVE:'',I3,''    RAY:'',I4,''    POINT:'',I4)')
     *                                                    IWAVE,IRAY,IPT
          END IF
          IRAYTT=IRAY
          XI1=YI(3)
          XI2=YI(4)
          XI3=YI(5)
          XE1=Y(3)
          XE2=Y(4)
          XE3=Y(5)
          CALL AP01(TT,TTAUX)
          CALL AP02(PI,PE)
          CALL AP29(NSUM,SUM)
        GO TO 30
C       End of the loop for endpoints of rays
C
   60   CONTINUE
        CLOSE(LU1)
        CLOSE(LU2)
      GO TO 20
C
C     All minimum travel times and their derivatives are stored in the
C     scratch file LU5.
C
C.......................................................................
C
C     Writing objective prior information:
C
   70 CONTINUE
      OPEN(LU4,FILE=FILE4)
      I=MAX0(INDEX(FILE4,'          ')-1,11)
      WRITE(*,'(''+Generating the output: '',A)') FILE4(1:I)
      ND=0
      DO 71 I=1,NT
        IF(TTMIN(I).LT.999999.) THEN
          ND=ND+1
        END IF
   71 CONTINUE
C     ND is the number of equations.
C
      REWIND(LU5)
      WRITE(LU4,'(I9,5I13)') ND,NM
      WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM)
      WRITE(LU7,'(2A)') ' SOURCE      RECEIVER       TFIELD    ',
     *                  'TFIELD-TT     SDIST       RDIST    TT-CHECKSUM'
   73 CONTINUE
        READ(LU5,END=79) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM)
        TTMIN(NT+1)=TTAUX
        IF(TTAUX.EQ.TTMIN(IT)) THEN
C         Minimum synthetic travel time
C
C         System of linear equations:
          TDIF=TFIELD(IT)-TTAUX
          WRITE(LU4,'(I9,4X,6G13.6)') IT,TDIF,TERR(IT),TFIELD(IT)
          WRITE(LU4,'(6G13.6)') (SUM(INDM(I)),I=1,NM)
C
C         Check sums and log output:
          IF(VPOWER.NE.0.) THEN
            TTAUX=0.
            DO 74 I=1,NM
              J=INDM(I)
              IF(IPAR(IPAR(IPAR(1))).LT.J) THEN
                IF(SUM(J).NE.0.) THEN
                  TTAUX=TTAUX+RPAR(J)*SUM(J)
                END IF
              END IF
   74       CONTINUE
            IS=KS(IT)
            IR=KR(IT)
            TTAUX=TT+VPOWER*TTAUX
            WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR),
     *           TFIELD(IT),TDIF,SDIST,RDIST,TTAUX
          ELSE
            WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR),
     *           TFIELD(IT),TDIF,SDIST,RDIST
          END IF
        END IF
      GO TO 73
C
   79 CONTINUE
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU7)
C
C.......................................................................
C
C     Subjective prior information:
C
   80 CONTINUE
      IF(FILE5.NE.' ') THEN
C
        IF(NM*(NM+1)/2.GT.MCS) THEN
          PAUSE 'Error: Too many model parameters'
          STOP
        END IF
C
C       Opening output file:
        OPEN(LU4,FILE=FILE5)
        I=MAX0(INDEX(FILE5,'          ')-1,11)
        WRITE(*,'('' Generating the output: '',A)') FILE5(1:I)
C
C       Reading prior subjective information coefficients:
        DO 81 I=1,MW+1
          NW1(I)=-1
   81   CONTINUE
        READ(LU0,*) (NW1(I),NW2(I),NW3(I),I=1,MW),NW1(MW+1)
        DO 82 I=1,MW
          IF(NW1(I).EQ.-1) THEN
            NW=I-1
            GO TO 83
          END IF
   82   CONTINUE
        PAUSE 'Error: Too large energy coefficient matrix'
        STOP
   83   CONTINUE
        READ(LU0,*) (WCS(I),I=1,NW*(NW+1)/2)
C
C       Generating prior subjective information covariance matrix:
        DO 84 I=1,MCS
          CS(I)=0.
   84   CONTINUE
        DO 89 J=1,NW
          DO 88 I=1,J
            WEIGHT=WCS(J*(J-1)/2+I)
            IF(WEIGHT.NE.0.) THEN
              CALL SOFT(2,NW1(I),NW2(I),NW3(I),NW1(J),NW2(J),NW3(J),
     *                                                WEIGHT,NM,INDM,CS)
            END IF
   88     CONTINUE
   89   CONTINUE
C
C       (1,2) Writing the initial model parameters:
        WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM)
        WRITE(LU4,'(6G13.6)') (RPAR(INDM(I)),I=1,NM)
C
C       (3) Writing relative posterior errors of the model parameters:
        DO 93 I=1,NM
          J=INDM(I)
          IF(IPAR(IPAR(IPAR(1))).LT.J) THEN
            SUM(I)=ERRVEL
          ELSE
            SUM(I)=ERRINT
          END IF
   93   CONTINUE
        WRITE(LU4,'(24F3.0)') (SUM(I),I=1,NM)
C
C       (4) Writing the prior subjective information covariance matrix:
        DO 94 J=1,NM
          WRITE(LU4,'(6G13.6)') (CS(I),I=J*(J-1)/2+1,J*(J+1)/2)
   94   CONTINUE
        CLOSE(LU4)
C
      END IF
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE SOFT(NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,WEIGHT,NM,INDM,CS)
C
      INTEGER NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,NM,INDM(*)
      REAL    WEIGHT,CS(*)
C
C This subroutine accumulates the prior subjective information
C covariance matrix describing the smoothness of the functions
C interpolated by means of subroutines of the file 'val.for'.
C
C Input:
C     NCLASS..Number of classes required.
C     NA1,NA2,NA3... Orders of X1,X2,X3-partial derivatives of the first
C             basis function in the product being integrated.
C     NB1,NB2,NB3... Orders of X1,X2,X3-partial derivatives of the
C             second basis function in the product being integrated.
C     WEIGHT..Weighting factor corresponding to the derivatives
C             (average quadratic value of the derivatives on input)
C     CS...   Symmetric matrix to be regularized:
C             Contribution corresponding to the derivatives
C             is added to the matrix CS
C None of the input parameters are altered.
C
C Output:
C     NM...   Number of the unknown model parameters.
C     INDM... Indices of the unknown model parameters.
C     CS...   Resulting symmetric subjective prior information
C             covariance matrix.
C
C Common block:
      INCLUDE 'val.inc'
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL SPSP
C
C Error messages:
C     359a..  Incorrect number of classes:
C             The number of classes required is zero, negative or
C             greater than the number of classes defined.
C     359b..  Insufficient memory:
C             The dimension of the buffer B is not sufficient.
C
C Date: 1992, September 7
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      INTEGER MM,MB
C     MM... Maximum number of model parameters.
      PARAMETER (MM=200)
      PARAMETER (MB=MM*(MM+1)/2+3*(20*21/2))
      REAL B(MB)
      INTEGER ICLASS,JGROUP,LFUNCT,MFUNCT,JFUNCT,LADR,MADR,IADR,NVAR
      INTEGER NX(3),NX1,NX2,NX3
      EQUIVALENCE (NX(1),NX1),(NX(2),NX2),(NX(3),NX3)
      INTEGER JADR(7),JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7
      EQUIVALENCE (JADR(1),JADR1),(JADR(2),JADR2),(JADR(3),JADR3)
      EQUIVALENCE (JADR(4),JADR4),(JADR(5),JADR5),(JADR(6),JADR6)
      EQUIVALENCE (JADR(7),JADR7)
      INTEGER I1,I2,I3,I,J,N
      REAL WDENS,SIGMA
C
C     B...    Auxiliary storage locations for partial covariance
C             matrices.
C     ICLASS..Index of the class.
C     JGROUP..Address of a current group.
C     LFUNCT,MFUNCT,JFUNCT... Addresses of the first, last and arbitrary
C             functions of the group.
C     LADR,MADR,IADR... Addresses of the first, last and arbitrary
C             parameters of the current function.
C     NVAR... Number of the independent variables A1, A2, A3 of the
C             interpolated function W.
C     NX=(NX1,NX2,NX3)... Numbers of grid lines.
C     JADR=(JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7)... Addresses of
C             parameters describing the interpolated function (grid
C             coordinates, B-spline coefficients, B-spline basis
C             functions).
C     I1,I2,I3,I,J,N... Local auxiliary variables.
C     WDENS...Weighting factor divided by volume (area,length).
C     SIGMA...Tension factor.
C
C.......................................................................
C
      NM=0
      IF(NCLASS.LT.1.OR.IPAR(0).LT.NCLASS) THEN
        PAUSE 'Error 359a in SOFT: Incorrect number of classes'
        STOP
      END IF
C     Loop for classes:
      DO 92 ICLASS=1,NCLASS
C       Loop for groups of the current class:
        DO 91 JGROUP=IPAR(ICLASS-1)+1,IPAR(ICLASS)
          LFUNCT=IPAR(JGROUP-1)+1
          MFUNCT=IPAR(JGROUP)
          MADR  =IPAR(LFUNCT-1)
C
C         Loop for functions of the current group:
          DO 90 JFUNCT=LFUNCT,MFUNCT
C           Starting and end addresses of the parameters describing the
C           function
            LADR=MADR+1
            MADR=IPAR(JFUNCT)
            IF(LADR.LE.MADR) THEN
              WDENS=WEIGHT
C             Tension factor
              SIGMA=RPAR(LADR+5)
C
C             The number, types and values of the independent variables
C             of the interpolated function:
C             Initial address
              IADR=LADR+6
C             Initial number of the independent variables
              NVAR=0
              NX1=1
              NX2=1
              NX3=1
              JADR1=0
              JADR2=0
              JADR3=0
              JADR4=0
C             Loop for the possible independent variables:
              DO 20 I=LADR+2,LADR+4
C               Type of the possible independent variable:
                J=IPAR(I)
                IF(J.GT.0) THEN
                  N=IABS(IPAR(IADR))
                  IF(N.GE.2) THEN
                    NVAR=NVAR+1
                    NX(NVAR)=N
                  ELSE IF(N.EQ.1) THEN
                    JADR(NVAR+1)=JADR(NVAR+1)+1
                  END IF
                  IADR=IADR+1
                END IF
   20         CONTINUE
              JADR5=0
              JADR6=0
              JADR7=0
C
C             Interpolated function W:
              JADR1=IADR+JADR1
              IF(NVAR.LE.0) THEN
C               No independent variable.
                JADR4=JADR1
              ELSE
                JADR2=JADR1+NX1+JADR2
                WDENS=WDENS/(RPAR(JADR2-1)-RPAR(JADR1))
                IF(NVAR.EQ.1) THEN
C                 One independent variable:
                  JADR4=JADR2
                  JADR5=JADR2+NX1
                ELSE
                  JADR3=JADR2+NX2+JADR3
                  WDENS=WDENS/(RPAR(JADR3-1)-RPAR(JADR2))
                  IF(NVAR.EQ.2) THEN
C                   Two independent variables:
                    JADR4=JADR3
                    JADR5=JADR3+NX1*NX2
                    JADR6=JADR4+5*NX1
                  ELSE
C                   Three independent variables:
                    JADR4=JADR3+NX3+JADR4
                    JADR5=JADR4+NX1*NX2*NX3
                    JADR6=JADR5+5*NX1
                    JADR7=JADR6+5*NX2
                    WDENS=WDENS/(RPAR(JADR4-1)-RPAR(JADR3))
                  END IF
                END IF
              END IF
C             WDENS is weighting factor divided by volume (area,length).
C
              JADR4=JADR4-1
              N=NX1*NX2*NX3
              IF(WEIGHT.EQ.0.) THEN
                DO 72 I2=1,N
                  INDM(NM+I2)=JADR4+I2
   72           CONTINUE
                NM=NM+N
              ELSE
                IF(NA1.EQ.NB1.AND.NA2.EQ.NB2.AND.NA3.EQ.NB3) THEN
                  I=N*(N+1)/2
                  J=0
                ELSE
                  I=N*N
                  J=1
                END IF
                IF(NA1.EQ.NB1) THEN
                  I1=NX1*(NX1+1)/2
                ELSE
                  I1=NX1*NX1
                END IF
                IF(NA2.EQ.NB2) THEN
                  I2=NX2*(NX2+1)/2
                ELSE
                  I2=NX2*NX2
                END IF
                IF(NA3.EQ.NB3) THEN
                  I3=NX3*(NX3+1)/2
                ELSE
                  I3=NX3*NX3
                END IF
                I1=I1+I
                I2=I2+I1
                I3=I3+I2
                IF(I3.GT.MB) THEN
                  PAUSE 'Error 359b in SOFT: Insufficient memory'
                  STOP
                END IF
                CALL SPSP(NA1,NA2,NA3,NB1,NB2,NB3,NX1,NX2,NX3,
     *                    RPAR(JADR1),RPAR(JADR2),RPAR(JADR3),
     *                    RPAR(JADR5),RPAR(JADR6),RPAR(JADR7),
     *                          SIGMA,WDENS,B(1),B(I+1),B(I1+1),B(I2+1))
                I=NM*(NM+1)/2
                DO 82 I2=1,N
                  INDM(NM+I2)=JADR4+I2
                  I=I+NM
                  DO 81 I1=1,I2
                    I=I+1
                    IF(J.EQ.0) THEN
                      CS(I)=CS(I)+B(I2*(I2-1)/2+I1)
                    ELSE
                      CS(I)=CS(I)+B(N*(I2-1)+I1)+B(N*(I1-1)+I2)
                    END IF
   81             CONTINUE
   82           CONTINUE
                NM=NM+N
C               Contribution corresponding to the derivatives
C               is added to the matrix CS.
              END IF
C
            END IF
   90     CONTINUE
   91   CONTINUE
   92 CONTINUE
C     End of loops for functions.
C
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'modelv.for'
      INCLUDE 'metric.for'
      INCLUDE 'srfc.for'
      INCLUDE 'parmv.for'
      INCLUDE 'valv.for'
      INCLUDE 'fitv.for'
      INCLUDE 'var.for'
      INCLUDE 'spsp.for'
      INCLUDE 'means.for'
      INCLUDE 'ap.for'
      INCLUDE 'apvar.for'
C
C=======================================================================
C
