C 
C Program MODMOD to modify the model (update or change parametrization)
C
C Version: 5.40
C Date: 2000, February 7
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 MODMOD 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 MODMOD 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 the original model:
C     MODEL='string'... String containing the name of the input data
C             file specifying the original model.
C             Description of file MODEL
C             Default: MODEL='model.dat'
C Data specifying the modification of the model:
C     M1='string'... Name of the input file containing the number NM of
C             model parameters (a single integer).
C             The file is generated by program 'invsoft.for' and is
C             used just if MODNEW is specified and is not blank or if
C             OLDMOD differs from 1.
C             Default: M1='m1.out'
C     MODIND='string'... Name of the input file containing the indices
C             of model parameters.
C             The file is generated by program 'invsoft.for' and is
C             used just if MODNEW is specified and is not blank or if
C             OLDMOD differs from 1.0.
C             File MODIND
C             Default: MODIND='modind.out'
C     MODNEW='string'... Name of the input file containing the updates
C             of the values of model parameters (coefficients at the
C             model basis functions).
C             If blank, original model MODEL is not updated.
C             File MODNEW
C             Default: MODNEW=' '
C     OLDMOD=real... Percentage of the original model MODEL kept in the
C             model.  For OLDMOD=1.0, original model is updated by the
C             values stored in file MODPAR.  For OLDMOD=0.0, original
C             model is discarded and replaced by the values stored in
C             file MODPAR (in such a case, MODPAR should contain
C             complete parameter values instead of their updates).
C             Default: OLDMOD=1.0
C Data specifying the form and name of the output model file:
C     MODIN='string'... Name of the file describing the form of the
C             parametrization of the output model.  If no changes in the
C             parametrization of the model are required, the default
C             value (value of parameter MODEL) is appropriate.
C             The functional values describing surfaces and material
C             parameters in file MODIN do not influence the resulting
C             model and may thus be arbitrary.
C             Description of file MODIN
C             Default: MODIN=MODEL
C     MODOUT='string'... Name of the output file describing the new
C             model.  File MODOUT is a copy of file MODIN, with the
C             functional values replaced by new ones.
C             Description of file MODOUT
C             Default: MODOUT='model.out'
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
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.......................................................................
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,MODEL1,NEWMOD,NEWVAL
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FILE1,FILE2,FILE3
      INTEGER LU1,LU2
      PARAMETER (LU1=1,LU2=2)
C
C     Auxiliary storage locations:
      INTEGER NSRF,NCB,M1,I
      REAL OLDMOD
      CHARACTER*3 TSRF(1),TCB(47)
      DATA TSRF/'   '/
      DATA TCB/'VP ','VS ','DEN','QP ','QS ',
     *'A11','A12','A22','A13','A23','A33','A14','A24','A34','A44',
     *'A15','A25','A35','A45','A55','A16','A26','A36','A46','A56','A66',
     *'Q11','Q12','Q22','Q13','Q23','Q33','Q14','Q24','Q34','Q44',
     *'Q15','Q25','Q35','Q45','Q55','Q16','Q26','Q36','Q46','Q56','Q66'/
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+MODMOD: Enter input filename: '
      FILE1=' '
      READ(*,*) FILE1
      IF(FILE1.EQ.' ') THEN
C       MODMOD-01
        CALL ERROR('MODMOD-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)') '+MODMOD: Working...            '
      CALL RSEP1(LU1,FILE1)
C
C     Checking input data MODIN:
      CALL RSEP3T('MODEL',FILE2,'model.dat')
      CALL RSEP3T('MODIN',FILE1,FILE2)
      IF(FILE1.NE.FILE2) THEN
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        CALL MODEL1(LU1)
        CLOSE(LU1)
        IPAR(0)=0
      END IF
C
C     Reading input data MODEL for the model to be updated:
      OPEN(LU2,FILE=FILE2,STATUS='OLD')
      CALL MODEL1(LU2)
      CLOSE(LU2)
C
C     Updating the model corresponding to data MODEL:
      CALL RSEP3T('MODNEW',FILE2,' ')
C     Reading percentage of old model parameters
      CALL RSEP3R('OLDMOD',OLDMOD,1.00)
      IF(FILE2.NE.' '.OR.OLDMOD.NE.1.0) THEN
        CALL RSEP3T('M1',FILE3,'m1.out')
        OPEN(LU2,FILE=FILE3,STATUS='OLD')
        READ(LU2,*) M1
        CLOSE(LU2)
        IF(2*M1.GT.MRAM) THEN
C         MODMOD-02
          CALL ERROR('MODMOD-02: Too many model parameters')
        END IF
C       Reading indices of model parameters
        CALL RSEP3T('MODIND',FILE3,'modind.out')
        OPEN(LU2,FILE=FILE3,STATUS='OLD')
        READ(LU2,*) (IRAM(I),I=1,M1)
        CLOSE(LU2)
        IF(FILE2.NE.' ') THEN
C         Reading increments of model parameters
          OPEN(LU2,FILE=FILE2,STATUS='OLD')
          READ(LU2,*) (RAM(I),I=M1+1,2*M1)
          CLOSE(LU2)
C         Updating the model
          DO 11 I=1,M1
            RPAR(IRAM(I))=RPAR(IRAM(I))*OLDMOD+RAM(M1+I)
   11     CONTINUE
        ELSE
C         Updating the model
          DO 12 I=1,M1
            RPAR(IRAM(I))=RPAR(IRAM(I))*OLDMOD
   12     CONTINUE
        END IF
      END IF
C
C     Converting input data MODIN into output data MODOUT:
      CALL RSEP3T('MODOUT',FILE2,'model.out' )
      OPEN(LU1,FILE=FILE1,STATUS='OLD')
      OPEN(LU2,FILE=FILE2)
      CALL NEWMOD(LU1,LU2,NSRF,NCB)
      CALL NEWVAL(LU1,LU2,1,NSRF,1,TSRF)
      CALL NEWVAL(LU1,LU2,2,NCB,47,TCB)
      CLOSE(LU1)
      CLOSE(LU2)
      WRITE(*,'(A)') '+MODMOD: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE NEWMOD(LU1,LU2,NSRF,NCB)
      INTEGER LU1,LU2,NSRF,NCB
C
C Subroutines and external functions required:
      EXTERNAL NEWLIN
C
C-----------------------------------------------------------------------
C
      CHARACTER*1 TEXTM
      INTEGER I,J,K,N,NEXPV,NEXPQ,NSB
      REAL AUX
C
C.......................................................................
C
C     Model description:
      N=0
   11 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=11) TEXTM
C
C     Model indices:
      N=0
   12 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
        NEXPV=1
        NEXPQ=1
      READ(LU2,*,END=12) I,NEXPV,NEXPQ,I
C
C     Model boundaries:
      N=0
   13 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=13) (AUX,I=1,6)
C
C     Number of surfaces:
      N=0
   14 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=14) NSRF
C
C     Number of simple blocks:
      N=0
   20 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=20) NSB
C
C     Indices of surfaces bounding simple blocks:
      DO 22 J=1,NSB
        N=0
   21   CONTINUE
          CALL NEWLIN(LU1,LU2,N)
        READ(LU2,*,END=21) (K,I=1,99)
   22 CONTINUE
C
C     Number of complex blocks:
      N=0
   30 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
      READ(LU2,*,END=30) NCB
C
C     Indices of simple blocks forming complex blocks:
      DO 32 J=1,NCB
        N=0
   31   CONTINUE
          CALL NEWLIN(LU1,LU2,N)
        READ(LU2,*,END=31) (K,I=1,99)
   32 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE NEWVAL(LU1,LU2,ICLASS,NGROUP,NFUNCT,TFUNCT)
      INTEGER LU1,LU2,ICLASS,NGROUP,NFUNCT
      CHARACTER*(*) TFUNCT(NFUNCT)
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 ERROR,WARRAY,VAL2,NEWLIN
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
      CHARACTER*3 TEXT
      LOGICAL WHAT
      INTEGER NVAR,IVAR(3),NX(3),MX,IGROUP,IFUNCT,JFUNCT,IADR
      INTEGER I1,I2,I3,I,N
      REAL GROUP,POWERW,COOR(3),F(10,47),POWER(47),AUX
C
C.......................................................................
C
C     Flag if the physical meaning of the functions is included in the
C     input data:
      WHAT=.FALSE.
      DO 10 I=1,NFUNCT
        IF(TFUNCT(I).NE.'   ') WHAT=.TRUE.
   10 CONTINUE
C
C     Loop for groups of functions:
      N=0
   11 CONTINUE
        CALL NEWLIN(LU1,LU2,N)
        GROUP=1.
      READ(LU2,*,END=11) TEXT,GROUP
      DO 90 IGROUP=1,NGROUP
C
C       Loop for functions of the current group:
        DO 80 IFUNCT=1,NFUNCT
C
C         Physical meaning of the function:
          IF(WHAT) THEN
            N=0
   21       CONTINUE
              CALL NEWLIN(LU1,LU2,N)
              GROUP=1.
            READ(LU2,*,END=21) TEXT,GROUP
            DO 22 I=1,NFUNCT
              IF(TFUNCT(I).EQ.TEXT) THEN
                JFUNCT=I
                GO TO 23
              END IF
   22       CONTINUE
            GO TO 89
   23       CONTINUE
          ELSE
            JFUNCT=IFUNCT
          END IF
C
C         Initial address of the function parameters:
          I2=IPAR(ICLASS-1)+IGROUP
          DO 25 I1=IPAR(I2-1)+1,IPAR(I2-1)+NFUNCT
            IADR=IPAR(I1-1)
            IF(IPAR(IADR+1).EQ.JFUNCT) THEN
              GO TO 26
            END IF
   25     CONTINUE
C           MODMOD-04
            CALL ERROR('MODMOD-04: Function not found')
C           Function specified in data MODIN has not been specified in
C           data MODEL.
   26     CONTINUE
C
C         Reading spline grid:
          DO 31 I=1,3
            IVAR(I)=0
            NX(I)=1
   31     CONTINUE
          N=0
   32     CONTINUE
            CALL NEWLIN(LU1,LU2,N)
            IVAR(1)=0
            IVAR(2)=0
            IVAR(3)=0
            POWERW=1.
          READ(LU2,*,END=32) (IVAR(I),I=1,3),AUX,POWERW
          NVAR=3
          I2=0
   41     CONTINUE
            I2=I2+1
            IF(IVAR(I2).LE.0) THEN
              NVAR=NVAR-1
              DO 42 I1=I2,NVAR
                IVAR(I1)=IVAR(I1+1)
   42         CONTINUE
              I2=I2-1
            END IF
          IF(I2.LT.NVAR) GO TO 41
          IF(NVAR.GT.0) THEN
            N=0
   44       CONTINUE
              CALL NEWLIN(LU1,LU2,N)
            READ(LU2,*,END=44) (NX(I),I=1,NVAR)
          END IF
          MX=MAX0(NX(1),NX(2),NX(3))
          RAM(     1)=0.
          RAM(  MX+1)=0.
          RAM(2*MX+1)=0.
          IF(4*MX.GT.MRAM) THEN
C           MODMOD-03
            CALL ERROR('MODMOD-03: Small array RAM')
          END IF
          DO 46 I2=1,NVAR
            IF(NX(I2).GT.0) THEN
              N=0
   45         CONTINUE
                CALL NEWLIN(LU1,LU2,N)
              READ(LU2,*,END=45)
     *                         (RAM(I1),I1=(I2-1)*MX+1,(I2-1)*MX+NX(I2))
            ELSE
              NX(I2)=1
            END IF
   46     CONTINUE
          READ(LU1,*) (AUX,I=1,NX(1)*NX(2)*NX(3))
C
C         Changing coordinate indices to 1,2,3:
          DO 53 I2=3,5
            IF(IPAR(IADR+I2).LE.0) THEN
              IPAR(IADR+I2)=0
            ELSE
              DO 51 I1=1,NVAR
                IF(IPAR(IADR+I2).EQ.IVAR(I1)) THEN
                  IPAR(IADR+I2)=I1
                  GO TO 52
                END IF
   51         CONTINUE
C               MODMOD-05
                CALL ERROR('MODMOD-05: Wrong independent variable')
C               Function in data MODEL depends on different variables
C               than the corresponding function in data MODIN.
   52         CONTINUE
            END IF
   53     CONTINUE
C
C         Calculating and writing grid values of the given function:
          DO 63 I3=1,NX(3)
            IF(NX(1).NE.1.AND.NX(2).NE.1.AND.NX(3).NE.1) THEN
C             Separating 2-D slices of 3-D grid by a blank line
              WRITE(LU2,*)
            END IF
            COOR(3)=RAM(2*MX+I3)
            DO 62 I2=1,NX(2)
              COOR(2)=RAM(MX+I2)
              DO 61 I1=1,NX(1)
                COOR(1)=RAM(I1)
                CALL VAL2(ICLASS,IGROUP,NFUNCT,COOR,F,POWER)
                AUX=GROUP*POWERW/POWER(JFUNCT)
                RAM(3*MX+I1)=F(1,JFUNCT)
                IF(WHAT) THEN
                  IF(AUX.NE.1.) THEN
                    IF(RAM(3*MX+I1).GE.0.) THEN
                      RAM(3*MX+I1)=RAM(3*MX+I1)**AUX
                    ELSE
C                     
C                     MODMOD-06
                      CALL ERROR('MODMOD-06: Power of a negative value')
C                     Function with negative values is interpolated
C                     while its non-unit power should be written.
C                     Such an operation is not permitted.
                    END IF
                  END IF
                END IF
   61         CONTINUE
              CALL WARRAY(LU2,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                                NX(1),RAM(3*MX+1))
   62       CONTINUE
   63     CONTINUE
   80   CONTINUE
C       End of loop for functions
C
        N=0
   81   CONTINUE
          CALL NEWLIN(LU1,LU2,N)
          GROUP=1.
        READ(LU2,*,END=81) TEXT,GROUP
   89   CONTINUE
   90 CONTINUE
C     End of loop for groups of functions
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE NEWLIN(LU1,LU2,N)
      INTEGER LU1,LU2,N
C
C Subroutines and external functions required:
      EXTERNAL LENGTH
      INTEGER  LENGTH
C
C-----------------------------------------------------------------------
C
      CHARACTER*160 LINE
      INTEGER I
C
C.......................................................................
C
C     Returning from the position after the end of file:
      IF(N.GT.0) THEN
        BACKSPACE(LU2)
      END IF
C
C     Copying one more line:
      READ (LU1,'(A)') LINE
      WRITE(LU2,'(A)') LINE(1:LENGTH(LINE))
      N=N+1
C
C     Rewinding to the position before reading:
      DO 10 I=1,N
        BACKSPACE(LU2)
   10 CONTINUE
      RETURN
      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
C
C=======================================================================
C