C 
C Program GREENTI designed to read two formatted output files of program
C GREEN containing Green tensors calculated by the prevailing-frequency
C approximation of the coupling ray theory along the SH-wave and SV-wave
C rays, and to generate new file containing the more accurate Green
C tensor composed of the two input Green tensors.
C
C Version: 7.20
C Date: 2015, May 15
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Faculty of Mathematics and Physics,
C     Charles University in Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
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 Names of input files:
C     REC='string'... The name of the file with the names
C             and coordinates of the receiver points.
C             Description of file REC.
C             Default: REC='rec.dat'
C     GREEN1='string'... Name of the input formatted file with the Green
C             tensor calculated by prevailing-frequency approximation of
C             the coupling ray theory along SH-wave rays (KSWAVE=1).
C             Description of file GREEN1.
C             Default: GREEN='green1.out'
C     GREEN2='string'... Name of the input formatted file with the Green
C             tensor calculated by prevailing-frequency approximation of
C             the coupling ray theory along SV-wave rays (KSWAVE=2).
C             Description of file GREEN2.
C             Default: GREEN='green2.out'
C Names of output files:
C     GREEN='string'... Name of the output formatted file with the Green
C             tensor.
C             Description of file GREEN.
C             Default: GREEN='green.out' (mostly sufficient)
C     GREENTI='string'... Name of the output formatted file containing
C             information on input Green tensors.
C             Description of file GREENTI.
C             Default: GREENTI='greenti.out'
C Data describing sorting of the Green tensors:
C     ISORT=integer... Specifies primary sorting criterion for the
C             output files with the Green tensors:
C             ISORT=0... Sorting according to the receiver index.
C             Default: ISORT=0 (sorting according to the receiver index)
C     TIA1=real, TIA2=real, TIA3=real:  Symmetry axis of a transversely
C             isotropic or approximately transversely isotropic (TI)
C             medium.
C             Default: TIA1=0, TIA2=0, TIA3=0 (no TI symmetry specified)
C
C-----------------------------------------------------------------------
C
C     Input and output files:
      INTEGER LU1,LU2,LU3,LU4,LU5
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5)
      CHARACTER*80 FILSEP,FILREC,FILE1,FILE2,FILE3,FILWRT
      CHARACTER*1 TEXT
C
C     Temporary storage locations:
      CHARACTER*11 REC,REC1,REC2,SRC1,SRC2
      INTEGER ISORT,I
      REAL TI1,TI2,TI3,GREEN1(33),GREEN2(33)
C
C.......................................................................
C
C     Reading input SEP file with input data:
      WRITE(*,'(A)') '+GREENTI: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU1,FILSEP)
      ELSE
C       GREENTI-01
        CALL ERROR('GREENTI-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.
      END IF
C
      WRITE(*,'(A)') '+GREENTI: Working...            '
C
C     Names of the input and output files:
      CALL RSEP3T('REC'   ,FILREC,'rec.dat')
      CALL RSEP3T('GREEN1',FILE1 ,'green1.out')
      CALL RSEP3T('GREEN2',FILE2 ,'green2.out')
      CALL RSEP3T('GREEN' ,FILE3 ,'green.out')
      CALL RSEP3T('GREENTI',FILWRT,'greenti.out')
C
C     Parameters describing sorting of the Green tensors:
      CALL RSEP3I('ISORT', ISORT, 0)
      IF(ISORT.NE.0) THEN
C       GREENTI-02
        CALL ERROR('GREENTI-02: Nonzero value of ISORT')
C       Program GREENTI assumes input SEP parameter ISORT=0.
      ENDIF
      CALL RSEP3R('TIA1',TI1,0.)
      CALL RSEP3R('TIA2',TI2,0.)
      CALL RSEP3R('TIA3',TI3,0.)
      IF(TI1.EQ.0..AND.TI2.EQ.0..AND.TI3.EQ.0.) THEN
C       GREENTI-03
        CALL ERROR('GREENTI-03: No TI symmetry axis specified')
C       Input SEP parameters TIA1,TIA2,TIA3 are probably not specified.
      ENDIF
C
C     Opening input and output files
      OPEN(LU1,FILE=FILREC,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
      OPEN(LU2,FILE=FILE1 ,STATUS='OLD')
      READ(LU2,*) (TEXT,I=1,20)
      OPEN(LU3,FILE=FILE2 ,STATUS='OLD')
      READ(LU3,*) (TEXT,I=1,20)
      OPEN(LU4,FILE=FILE3)
      WRITE(LU4,'(A)') '/'
      OPEN(LU5,FILE=FILWRT)
C
C     Writing the header of the report
      WRITE(LU5,'(2A)') 'receiver  rays',
     *                  ' SV-1rec  SV-1ini  SV-2rec  SV-2ini'
C
C     Reading the first elementary Green tensors of the input files:
      REC1='$'
      READ(LU2,*,ERR=90,END=90) REC1,SRC1,(GREEN1(I),I=1,33)
      REC2='$'
      READ(LU3,*,ERR=90,END=90) REC2,SRC2,(GREEN2(I),I=1,33)
C
C     Loop over receivers:
   10 CONTINUE
        REC='$'
        READ(LU1,*,ERR=90,END=90) REC
        IF(REC.EQ.'$') THEN
C         End of receivers
          GO TO 90
        END IF
C
C       Loop over the pairs of the SH elementary Green tensors
   20   CONTINUE
          IF(REC1.EQ.REC) THEN
            CALL TIWRIT(LU2,LU4,LU5,REC,'SH',GREEN1,TI1,TI2,TI3)
          END IF
          CALL TIREAD(LU2,REC1,SRC1,GREEN1)
        IF(REC1.EQ.REC) GO TO 20
C
C       Loop over the pairs of the SV elementary Green tensors
   30   CONTINUE
          IF(REC2.EQ.REC) THEN
            CALL TIWRIT(LU3,LU4,LU5,REC,'SV',GREEN2,TI1,TI2,TI3)
          END IF
          CALL TIREAD(LU3,REC2,SRC2,GREEN2)
        IF(REC2.EQ.REC) GO TO 30
      GO TO 10
C
   90 CONTINUE
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
      WRITE(LU4,'(A)') '/'
      CLOSE(LU4)
      CLOSE(LU5)
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE TIREAD(LU1,REC1,SRC1,GREEN1)
      INTEGER LU1
      CHARACTER*(*) REC1,SRC1
      REAL GREEN1(33)
C
      REAL UARRAY
      EXTERNAL UARRAY
C
C-----------------------------------------------------------------------
C
      INTEGER I
C
C     Undefined value:
      REAL UNDEF
C
      UNDEF=UARRAY()
C
      REC1='$'
      GREEN1(33)=UNDEF
      READ(LU1,*) REC1,SRC1,(GREEN1(I),I=1,33)
      IF(GREEN1(33).NE.UNDEF) THEN
C       Frequency-dependent Green tensor:
C       GREENTI-04
        CALL ERROR('GREENTI-04: Frequency dependent Green tensor')
C       Frequency-dependent Green tensor has been identified in
C       the input file instead of the prevailing-frequency
C       approximation of the Green tensor.
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TIWRIT(LU2,LU4,LU5,REC,WAVE,GREEN1,TI1,TI2,TI3)
      INTEGER LU2,LU4,LU5
      CHARACTER*(*) REC,WAVE
      REAL GREEN1(33),TI1,TI2,TI3
C
      INTEGER LENGTH
      EXTERNAL LENGTH
C
C-----------------------------------------------------------------------
C
      CHARACTER*260 FORMAT
      CHARACTER*11 REC3,SRC3
      INTEGER I
      REAL GREEN3(33),C1C1,S1S1,C2C2,S2S2
C
C     Reading the even elementary Green tensor
      REC3='$'
      READ(LU2,*,END=90) REC3,SRC3,(GREEN3(I),I=1,33)
      IF(REC3.EQ.REC) THEN
        CALL TIPRD(GREEN1,TI1,TI2,TI3,C1C1,C2C2)
        CALL TIPRD(GREEN3,TI1,TI2,TI3,S1S1,S2S2)
        C1C1=C1C1/(C1C1+S1S1)
        C2C2=C2C2/(C2C2+S2S2)
C
C       Writing the output Green tensor
        FORMAT(1:4)='(6A,'
        IF((WAVE.EQ.'SH'.AND.C1C1+C2C2.LE.1.).OR.
     *     (WAVE.EQ.'SV'.AND.C1C1+C2C2.GT.1.)) THEN
          CALL FORM2(32,GREEN1,GREEN1,FORMAT(5:260))
          WRITE(LU4,FORMAT) '''',REC(1:LENGTH(REC)),''' ''',
     *                           SRC3(1:LENGTH(SRC3)),'''',
     *                      (' ',GREEN1(I),I=1,32),' /'
        ELSE
          CALL FORM2(32,GREEN3,GREEN3,FORMAT(5:260))
          WRITE(LU4,FORMAT) '''',REC(1:LENGTH(REC)),''' ''',
     *                           SRC3(1:LENGTH(SRC3)),'''',
     *                      (' ',GREEN3(I),I=1,32),' /'
        END IF
C
C       Writing the report
        WRITE(LU5,'(3A,4F9.6)') REC,' ',WAVE,C2C2,C1C1,1.-C2C2,1.-C1C1
        RETURN
      END IF
C
C     GREENTI-05
      CALL ERROR('GREENTI-05: Input Green tensors are not paired')
C     Input prevaling-frequency elementary Green tensors must be written
C     in pairs for each receiver.
      RETURN
C
   90 CONTINUE
C     GREENTI-06
      CALL ERROR('GREENTI-06: End of input Green tensor')
C     Input prevaling-frequency elementary Green tensors must be written
C     in pairs for each receiver.
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TIPRD(GREEN,TI1,TI2,TI3,C1C1,C2C2)
      REAL GREEN(33),TI1,TI2,TI3,C1C1,C2C2,AMPL
C
C-----------------------------------------------------------------------
C
      AMPL=GREEN(15)**2+GREEN(21)**2+GREEN(27)**2
     *    +GREEN(16)**2+GREEN(22)**2+GREEN(28)**2
     *    +GREEN(17)**2+GREEN(23)**2+GREEN(29)**2
     *    +GREEN(18)**2+GREEN(24)**2+GREEN(30)**2
     *    +GREEN(19)**2+GREEN(25)**2+GREEN(31)**2
     *    +GREEN(20)**2+GREEN(26)**2+GREEN(32)**2
      C1C1=(GREEN(15)*TI1+GREEN(21)*TI2+GREEN(27)*TI3)**2
     *    +(GREEN(16)*TI1+GREEN(22)*TI2+GREEN(28)*TI3)**2
     *    +(GREEN(17)*TI1+GREEN(23)*TI2+GREEN(29)*TI3)**2
     *    +(GREEN(18)*TI1+GREEN(24)*TI2+GREEN(30)*TI3)**2
     *    +(GREEN(19)*TI1+GREEN(25)*TI2+GREEN(31)*TI3)**2
     *    +(GREEN(20)*TI1+GREEN(26)*TI2+GREEN(32)*TI3)**2
      C2C2=(GREEN(15)*TI1+GREEN(17)*TI2+GREEN(19)*TI3)**2
     *    +(GREEN(16)*TI1+GREEN(18)*TI2+GREEN(20)*TI3)**2
     *    +(GREEN(21)*TI1+GREEN(23)*TI2+GREEN(25)*TI3)**2
     *    +(GREEN(22)*TI1+GREEN(24)*TI2+GREEN(26)*TI3)**2
     *    +(GREEN(27)*TI1+GREEN(29)*TI2+GREEN(31)*TI3)**2
     *    +(GREEN(28)*TI1+GREEN(30)*TI2+GREEN(32)*TI3)**2
      C1C1=C1C1/AMPL
      C2C2=C2C2/AMPL
      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
C
C=======================================================================
C