C 
C Program TCCOMP to read two files with propagator matrices  U2 and U3
C written in the form of file GREENTC,
C to compute matrix DU=(U2-U3), and to compute its 'norm' as specified
C by input parameter INORM.
C
C Version: 5.60
C Date: 2002, May 23
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
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 a SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Name of the input file with the propagation matrices:
C     GREENTC2='string'... Name of the file with matrix U2.
C         Description of file GREENTC2
C             Default: GREENTC2=' '
C     GREENTC='string'... Name of the file with matrix U3.
C         Description of file GREENTC3
C             Default: GREENTC=' '
C Name of the output file:
C     TCCOMP='string'... Name of the output formatted file with
C             the norm of matrix DU specified by INORM.
C             Default: TCCOMP='tccomp.out'
C Specification of the norm to be used:
C     INORM=integer
C             INORM=1   ... output quantity ADU is the norm
C               of the matrix DU divided by square root of 2.
C             otherwise ... output quantity ADU is calculated as
C               ADU=SQRT(2.*|DU|/(|U2|+|U3|)), where |DU| stands for
C               the square of the norm of the matrix DU.
C             Default: INORM=1
C Output formatted file TCCOMP:
C   For each frequency following line:
C     F,ADU,TEXT
C   where F is the frequency, ADU is the 'norm' specified by INORM,
C   and text is a character enabling the output file to be included
C   into PostScript file. Thus the first line terminates with
C   character M, and other lines terminate with character L.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T
C     ERROR ... File error.for.
C     RSEP1,RSEP3T ... File sep.for.
C
C-----------------------------------------------------------------------
C
C     Input and output data files:
      CHARACTER*80 FSEP,FIN2,FIN3,FOUT
      CHARACTER*2 TEXT
      INTEGER LU,LU2,LU3,LU4,INORM
      PARAMETER (LU=5,LU2=2,LU3=3,LU4=4)
C     Auxiliary storage locations:
      REAL RV11,RV21,RV31,RV12,RV22,RV32,RV13,RV23,RV33
      REAL RW11,RW21,RW31,RW12,RW22,RW32,RW13,RW23,RW33
      REAL RX11,RX21,RX31,RX12,RX22,RX32,RX13,RX23,RX33
      REAL IV11,IV21,IV31,IV12,IV22,IV32,IV13,IV23,IV33
      REAL IW11,IW21,IW31,IW12,IW22,IW32,IW13,IW23,IW33
      REAL IX11,IX21,IX31,IX12,IX22,IX32,IX13,IX23,IX33
      REAL F,ADU
C
C.......................................................................
C
C     Main input data:
      FSEP=' '
      WRITE(*,'(A)') '+TCCOMP: Enter input filename: '
      READ(*,*) FSEP
      IF (FSEP.EQ.' ') THEN
C       TCCOMP-01
        CALL ERROR('TCCOMP-01: SEP file not given')
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.
      ENDIF
C
      WRITE(*,'(A)') '+TCCOMP: Working...            '
C
C     Reading all the data from file FSEP to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
C
C     Recalling the data:
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3T('GREENTC2',FIN2,' ')
      CALL RSEP3T('GREENTC' ,FIN3,' ')
      CALL RSEP3T('TCCOMP',FOUT,'tccomp.out')
      CALL RSEP3I('INORM',INORM,1)
      OPEN(LU2,FILE=FIN2,STATUS='OLD')
      OPEN(LU3,FILE=FIN3,STATUS='OLD')
      OPEN(LU4,FILE=FOUT)
      TEXT=' M'
  10  CONTINUE
        READ(LU2,*,END=100) F,RV11,IV11,RV21,IV21,RV31,IV31,RV12,IV12,
     *  RV22,IV22,RV32,IV32,RV13,IV13,RV23,IV23,RV33,IV33
        READ(LU3,*,END=100) F,RW11,IW11,RW21,IW21,RW31,IW31,RW12,IW12,
     *  RW22,IW22,RW32,IW32,RW13,IW13,RW23,IW23,RW33,IW33
        RX11=RV11-RW11
        IX11=IV11-IW11
        RX21=RV21-RW21
        IX21=IV21-IW21
        RX31=RV31-RW31
        IX31=IV31-IW31
        RX12=RV12-RW12
        IX12=IV12-IW12
        RX22=RV22-RW22
        IX22=IV22-IW22
        RX32=RV32-RW32
        IX32=IV32-IW32
        RX13=RV13-RW13
        IX13=IV13-IW13
        RX23=RV23-RW23
        IX23=IV23-IW23
        RX33=RV33-RW33
        IX33=IV33-IW33
        ADU=RX11**2+IX11**2+RX21**2+IX21**2+RX31**2+IX31**2+
     *      RX12**2+IX12**2+RX22**2+IX22**2+RX32**2+IX32**2+
     *      RX13**2+IX13**2+RX23**2+IX23**2+RX33**2+IX33**2
        IF (INORM.EQ.1) THEN
          ADU=SQRT(0.5*ADU)
        ELSE
          ADV=RV11**2+IV11**2+RV21**2+IV21**2+RV31**2+IV31**2+
     *        RV12**2+IV12**2+RV22**2+IV22**2+RV32**2+IV32**2+
     *        RV13**2+IV13**2+RV23**2+IV23**2+RV33**2+IV33**2
          ADW=RW11**2+IW11**2+RW21**2+IW21**2+RW31**2+IW31**2+
     *        RW12**2+IW12**2+RW22**2+IW22**2+RW32**2+IW32**2+
     *        RW13**2+IW13**2+RW23**2+IW23**2+RW33**2+IW33**2
          ADU=SQRT(2.*ADU/(ADV+ADW))
        ENDIF
        WRITE(LU4,'(1(F9.6,1X),F9.6,A)')
     *    F,ADU,TEXT
        TEXT=' L'
      GOTO 10
  100 CONTINUE
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      WRITE(*,'(A)') '+TCCOMP: Done.                  '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
C=======================================================================
C