C 
C Program INV3 to update the input data for a function describing the
C model.
C
C Version: 5.20
C Date: 1998, October 20
C
C Coded by Ludek Klimes
C
C ......................................................................
C
C Just a preliminary demo version, generating a table of values of a
C given function describing the model.  The given function may be a
C function describing a sooth surface covering a structural interface,
C or a function describing the spatial distribution of a material
C parameter.  The function is evaluated at gridpoints of a given
C rectangular grid of points.  The model parameters (coefficients of
C functions) may be updated by increments resulting from an inversion of
C the system of equations generated by the INV1 program.  Consequently,
C the tables generated by the INV3 program may be used to update the
C input data for the model by means of manual editing.  The inversion
C program (reading results 'DATA' and 'SOFT' of program INV1, and
C generating input 'ANSWER' of program INV3) should be contributed by a
C user.
C
C Program INV3 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 INV3 program cannot work with user's modifications of the
C subroutines SRFC1, SRFC2, PARM1, and PARM2.
C
C.......................................................................
C
C                                                    
C Description of the data files:
C
C Main input data read from the interactive device (*):
C (1) 'MODEL','ANSWER','INV3IN','INV3OUT',/
C     'MODEL','ANSWER','GRID'... Names of the input and output
C             files described below.
C     'MODEL'... Input data file containing the model parameters.
C     'ANSWER'... Updates to the model.  If blank, model is not updated.
C     'INV3IN'... Input file specifying the grid in which a selected
C             function is discretized.
C     'INV3OUT'... Copy of the specification of the grid in which the
C             selected function is discretized from INV3IN,
C             supplemented with the function values at gridpoints.
C     /...    Obligatory slash to enable future compatible extensions.
C     Default: 'MODEL'='model.dat', 'ANSWER'=' ', 'INV3IN'='inv3.dat',
C             'INV3OUT'='inv3.out'.
C
C Input file MODEL:
C Input data file containing the model parameters.  See file 'model.for'
C of package MODEL.
C Description of file MODEL
C
C Input file ANSWER:
C (1) NM
C     NM...   Number of model parameters.
C             NM=0: initial model is discretized into the given grid.
C             Default: NM=0.
C (2) NM-times the following data:
C     INDM,RM
C     INDM... Index of a model parameter.
C     RM...   Increment of the model parameter.
C (3) (CM(I,J),I=1,J),J=1,NM)
C     CM...   Model parameter covariance matrix including the subjective
C             prior information.
C
C Input file INV3IN:
C (1) TEXTG,IGROUP
C     Identification of the group.
C     TEXTG...A string.  If its first character is 'I' or 'S', the
C             computed function describes a surface.  Otherwise, it
C             describes a material parameter.
C     IGROUP..Index of a surface, or of a complex block.
C (2) TEXTF
C     This input is not performed for a surface, i.e. if the first
C     character of TEXTG (see above) is 'I' or 'S'.
C     TEXTF...String identifying a material parameter.
C (3) K1,K2,K3
C     K1,K2,K3... Indices of coordinates.
C (4) N1,N2,N3
C     N1,N2,N3... Numbers of grid lines.
C (5) X1(1),...,X1(N1)
C     The grid coordinates corresponding to the first independent
C     variable.
C (6) X2(1),...,X2(N2)
C     The grid coordinates corresponding to the second independent
C     variable.
C (7) X3(1),...,X3(N3)
C     The grid coordinates corresponding to the third independent
C     variable.
C
C Output file INV3OUT:
C (1)-(7): Copy of (1)-(7) from input file 'INV3IN'.
C (8) (((W(I1,I2,I3),I1=1,N1),I2=1,N2),I3=1,N3)
C     The values of function W at grid points. Function value
C     W(I1,I2,I3) corresponds to point (X1(I1),X2(I2),X3(I3)).
C (9) 'STANDARD DEVIATIONS ....................' (a string)
C (10) (((E(I1,I2,I3),I1=1,N1),I2=1,N2),I3=1,N3)
C     Standard deviations of function values.
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 Subroutines and external functions required:
      EXTERNAL LENGTH
      INTEGER  LENGTH
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working arrays:
      INTEGER MX,ME
      PARAMETER (MX=100,ME=8000,MCM=MRAM-NPAR-4*MX-ME)
      INTEGER INDM(NPAR)
      REAL X(MX,3),W(MX),E(ME)
      REAL CM(MCM)
      EQUIVALENCE (INDM ,RAM                )
      EQUIVALENCE (X    ,RAM(NPAR        +1))
      EQUIVALENCE (W    ,RAM(NPAR+3*MX   +1))
      EQUIVALENCE (E    ,RAM(NPAR+4*MX   +1))
      EQUIVALENCE (CM   ,RAM(NPAR+4*MX+ME+1))
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILE1,FILE2,FILE3,FILE4
C
C     Logical unit numbers:
      INTEGER LU1,LU2,LU3
      PARAMETER (LU1=11)
      PARAMETER (LU2=12)
      PARAMETER (LU3=13)
C
C     A line of the copied file:
      CHARACTER*160 LINE
C
C     Input data:
      INTEGER NM
      INTEGER IGROUP,K1,K2,K3,N1,N2,N3
C
C     Auxiliary storage locations:
      CHARACTER*3 TEXT
      INTEGER MFUN
      PARAMETER (MFUN=64)
      INTEGER I1,I2,I3,IVAL,IFUN(MFUN),NFUN,IE,II,JJ,IJJ
      REAL COOR(3),UP(10),US(10),AUX0,AUX1,AUX2
      REAL FUN(MFUN),B0I,B1I,B2I,B3I
C
      INTEGER MM
C     MM... Maximum number of model parameters.
      MM=INT(SQRT(FLOAT(2*MCM)+0.25)-0.5)
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') ' Enter names of input and output files: '
      FILE1='model.dat'
      FILE2=' '
      FILE3='inv3.dat'
      FILE4='inv3.out'
      READ(*,*) FILE1,FILE2,FILE3
      WRITE(*,'(A)') '+                                       '
C
C     Reading input data for model:
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
C
C     Updating the model:
      NM=0
      IF(FILE2.NE.' ') THEN
        OPEN(LU2,FILE=FILE2,STATUS='OLD')
        READ(LU2,*) NM
        IF(NM.GT.MM) THEN
C         INV3-01
          CALL ERROR('INV3-01: Too many model parameters')
        END IF
        DO 1 I1=1,NM
          READ(LU2,*) INDM(I1),AUX0
          IF(INDM(I1).GT.0) THEN
            RPAR(INDM(I1))=RPAR(INDM(I1))+AUX0
          END IF
    1   CONTINUE
        IF(NM.GT.0) THEN
          READ(LU2,*) (CM(I1),I1=1,NM*(NM+1)/2)
        END IF
        CLOSE(LU2)
      END IF
C
C     Copying input file 'INV3IN' to output file 'INV3OUT':
      OPEN(LU2,FILE=FILE3,STATUS='OLD')
      OPEN(LU3,FILE=FILE4)
    9 CONTINUE
        READ(LU2,'(A)',END=10) LINE
        WRITE(LU3,'(A)') LINE(1:LENGTH(LINE))
      GO TO 9
   10 CONTINUE
      CLOSE(LU2)
      REWIND(LU3)
C     CLOSE(LU3)
C
C     Reading function and grid specifications:
C     OPEN(LU3,FILE=FILE4,STATUS='OLD')
      READ(LU3,*) TEXT,IGROUP
      IF(TEXT(1:1).NE.'I'.AND.TEXT(1:1).NE.'S') THEN
        READ(LU3,*) TEXT
      END IF
      READ(LU3,*) K1,K2,K3
      READ(LU3,*) N1,N2,N3
      IF(MAX0(N1,N2,N3).GT.MX) THEN
C       INV3-02
        CALL ERROR('INV3-02: Too many grid lines')
      END IF
      IF(N1*N2*N3.GT.ME) THEN
C       INV3-03
        CALL ERROR('INV3-03: Too large grid')
      END IF
      READ(LU3,*) (X(I1,K1),I1=1,N1)
      READ(LU3,*) (X(I2,K2),I2=1,N2)
      READ(LU3,*) (X(I3,K3),I3=1,N3)
C
C     Evaluating grid values of the given function:
      IE=0
      DO 23 I3=1,N3
        COOR(K3)=X(I3,K3)
        DO 22 I2=1,N2
          COOR(K2)=X(I2,K2)
          DO 21 I1=1,N1
C
C           Evaluating the functional value:
            COOR(K1)=X(I1,K1)
            IF(TEXT(1:1).EQ.'I'.OR.TEXT(1:1).EQ.'S') THEN
              CALL SRFC2(IGROUP,COOR,UP)
              IVAL=1
              W(I1)=UP(1)
            ELSE
              CALL PARM2(IGROUP,COOR,UP,US,AUX0,AUX1,AUX2)
              IF(TEXT.EQ.'VP ') THEN
                IVAL=1
                W(I1)=UP(1)
              ELSE IF(TEXT.EQ.'VS ') THEN
                IVAL=2
                W(I1)=US(1)
              ELSE IF(TEXT.EQ.'DEN') THEN
                IVAL=3
                W(I1)=AUX0
              ELSE IF(TEXT.EQ.'QP ') THEN
                IVAL=4
                W(I1)=AUX1
              ELSE IF(TEXT.EQ.'QS ') THEN
                IVAL=5
                W(I1)=AUX2
              ELSE
C               INV3-04
                CALL ERROR
     *              ('INV3-04: Name of medium parameter not recognized')
              END IF
            END IF
C
C           Evaluating the standard deviation:
            IF(NM.GT.0) THEN
              II=0
   11         CONTINUE
                II=II+1
                CALL VAR6(IVAL,II,NFUN,IFUN(II),B0I,B1I,B2I,B3I)
                IF(II.LE.NFUN) THEN
                  IF(NFUN.GT.MFUN) THEN
C                   
C                   INV3-05
                    CALL ERROR('INV3-05: Array index out of range')
                  END IF
                  FUN(II)=B0I
                END IF
              IF(II.LT.NFUN) GO TO 11
              DO 14 JJ=1,NFUN
                DO 12 II=1,NM
                  IF(INDM(II).EQ.IFUN(JJ)) THEN
                    IFUN(JJ)=II
                    GO TO 13
                  END IF
   12           CONTINUE
C               INV3-06
                CALL ERROR
     *                 ('INV3-06: Model parameter index not recognized')
   13           CONTINUE
   14         CONTINUE
              AUX0=0.
              DO 17 JJ=1,NFUN
                IJJ=IFUN(JJ)*(IFUN(JJ)-1)/2
                AUX2=0.
                DO 16 II=1,JJ-1
                  AUX2=AUX2+     CM(IJJ+IFUN(II))*FUN(II)
   16           CONTINUE
                AUX0=AUX0+FUN(JJ)*(CM(IJJ+IFUN(JJ))*FUN(JJ)+2.*AUX2)
   17         CONTINUE
              IE=IE+1
              E(IE)=SQRT(AUX0)
            END IF
C
   21     CONTINUE
C         WRITE(LU3,'(8F10.6)') (W(I1),I1=1,N1)
          CALL WARRAY(LU3,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,N1,W)
   22   CONTINUE
        IF(N1.NE.1.AND.N2.NE.1) WRITE(LU3,*)
   23 CONTINUE
C
      IF(NM.GT.0) THEN
        WRITE(LU3,'(A)') '''STANDARD DEVIATIONS ....................'''
        IE=0
        DO 33 I3=1,N3
          DO 32 I2=1,N2
            WRITE(LU3,'(8F10.6)') (E(I1),I1=IE+1,IE+N1)
            IE=IE+N1
   32     CONTINUE
          IF(N1.NE.1.AND.N2.NE.1) WRITE(LU3,*)
   33   CONTINUE
      END IF
C
      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 '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
C
C=======================================================================
C