C
C Program GREENTC to multiply Green function by factor exp(2*Pi*f*t).
C The program reads given formatted file GREEN containing
C the ray-theory elastodynamic Green function in the form of the output
C file GREEN
C generated by program GREEN. This program multiplies the amplitudes of
C the Green function by factor exp(2*Pi*f*t), and writes the components
C of the resulting propagator matrix into the file GREENTC.
C
C Version: 7.40
C Date: 2017, May 18
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.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 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 Green function:
C     GREEN='string'... Name of the input formatted file with the
C             Green tensor.
C             Description of file
C        GREEN.
C             Default: GREEN='green.out'
C Name of the output file:
C     GREENTC='string'... Name of the output formatted file with the
C             propagator matrix.
C             Description of file GREENTC.
C             Default: GREENTC='greentc.out'
C Data describing the initial travel time in the source:
C     OT=real ... Initial travel time.
C             Default: OT=0.
C Data describing the frequency domain:
C     DF=real ... Frequency step.
C             Default: DF=0.
C     OF=real ... The elementary Green functions are calculated for NF
C             frequencies OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF.
C             Default: OF=0.
C     NF=integer ... Number of frequencies.
C             Default: NF=1
C Value of undefined quantities:
C     UNDEF=real... The value to be used for undefined real quantities.
C             Default: UNDEF=undefined value used in forms.for
C
C                                                 
C Output formatted file GREENTC:
C   For each frequency following line:
C     F,REU11,IMU11,REU21,IMU21,REU31,IMU31,
C       REU12,IMU12,REU22,IMU22,REU32,IMU32,
C       REU13,IMU13,REU23,IMU23,REU33,IMU33,
C   where F is the frequency and REUii,IMUii are real and imaginary
C   parts of the 3x3 ray propagator matrix.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,RSEP3R,UARRAY
      REAL UARRAY
C     ERROR ... File
C     error.for.
C     UARRAY ... File
C     forms.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C     File sep.for.
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      REAL GREEN(MRAM)
      EQUIVALENCE (GREEN,RAM)
C
C Auxiliary storage locations:
      REAL PI,UNDEF
      PARAMETER (PI=3.1415926)
      REAL REG11,IMG11,REG21,IMG21,REG31,IMG31,REG12,IMG12,REG22,IMG22,
     *     REG32,IMG32,REG13,IMG13,REG23,IMG23,REG33,IMG33
      REAL REU11,IMU11,REU21,IMU21,REU31,IMU31,
     *     REU12,IMU12,REU22,IMU22,REU32,IMU32,
     *     REU13,IMU13,REU23,IMU23,REU33,IMU33
      REAL F,AUX,REAUX,IMAUX,TT
      CHARACTER*80 FILSEP,FILIN,FILOUT
      CHARACTER*80 TEXT,TXTSRC,TXTREC
      INTEGER NF,NGREEN
      REAL OF,DF,OT
      INTEGER LU1,LU2,I1,I2,I
      PARAMETER (LU1=1,LU2=2)
C
C-----------------------------------------------------------------------
C
C     Main input data:
      FILSEP=' '
      WRITE(*,'(A)') '+GREENTC: Enter input filename: '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       GREENTC-01
        CALL ERROR('GREENTC-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)') '+GREENTC: Working...            '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
      UNDEF=UARRAY()
C
C     Name of the input file:
      CALL RSEP3T('GREEN',FILIN,'green.out')
C     Name of the output file:
      CALL RSEP3T('GREENTC',FILOUT,'greentc.out')
C     Number of frequencies:
      CALL RSEP3I('NF',NF,1)
C     Origin of frequencies:
      CALL RSEP3R('OF',OF,0.)
C     Step in frequencies:
      CALL RSEP3R('DF',DF,0.)
C     Origin of travel time:
      CALL RSEP3R('OT',OT,0.)
C
C     Opening input file with the Green function:
      OPEN(LU1,FILE=FILIN,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
C
C     Opening the output file:
      OPEN(LU2,FILE=FILOUT)
C
C     Loop over the records in file GREEN(MUL):
  10  CONTINUE
        NGREEN=14+18*NF
        DO 12, I1=1,NGREEN
          GREEN(I1)=0.
  12    CONTINUE
        GREEN(33)=UNDEF
C       Reading:
        TXTREC='$'
        READ(LU1,*) TXTREC,TXTSRC,(GREEN(I),I=1,NGREEN)
        IF (TXTREC.EQ.'$') GOTO 20
        IF(GREEN(33).EQ.UNDEF) THEN
          DO 13, I1=33,NGREEN
            GREEN(I1)=GREEN(MOD(I1-15,18)+15)
  13      CONTINUE
        END IF
        TT=GREEN(1)-OT
        I2=15
        DO 14, I1=1,NF
          F=OF+(I1-1)*DF
          REG11=GREEN(I2   )/1000000.
          IMG11=GREEN(I2+ 1)/1000000.
          REG21=GREEN(I2+ 2)/1000000.
          IMG21=GREEN(I2+ 3)/1000000.
          REG31=GREEN(I2+ 4)/1000000.
          IMG31=GREEN(I2+ 5)/1000000.
          REG12=GREEN(I2+ 6)/1000000.
          IMG12=GREEN(I2+ 7)/1000000.
          REG22=GREEN(I2+ 8)/1000000.
          IMG22=GREEN(I2+ 9)/1000000.
          REG32=GREEN(I2+10)/1000000.
          IMG32=GREEN(I2+11)/1000000.
          REG13=GREEN(I2+12)/1000000.
          IMG13=GREEN(I2+13)/1000000.
          REG23=GREEN(I2+14)/1000000.
          IMG23=GREEN(I2+15)/1000000.
          REG33=GREEN(I2+16)/1000000.
          IMG33=GREEN(I2+17)/1000000.
          AUX=2.*PI*F*TT
          REAUX=COS(AUX)
          IMAUX=SIN(AUX)
          REU11=REAUX*REG11-IMAUX*IMG11
          IMU11=REAUX*IMG11+IMAUX*REG11
          REU21=REAUX*REG21-IMAUX*IMG21
          IMU21=REAUX*IMG21+IMAUX*REG21
          REU31=REAUX*REG31-IMAUX*IMG31
          IMU31=REAUX*IMG31+IMAUX*REG31
          REU12=REAUX*REG12-IMAUX*IMG12
          IMU12=REAUX*IMG12+IMAUX*REG12
          REU22=REAUX*REG22-IMAUX*IMG22
          IMU22=REAUX*IMG22+IMAUX*REG22
          REU32=REAUX*REG32-IMAUX*IMG32
          IMU32=REAUX*IMG32+IMAUX*REG32
          REU13=REAUX*REG13-IMAUX*IMG13
          IMU13=REAUX*IMG13+IMAUX*REG13
          REU23=REAUX*REG23-IMAUX*IMG23
          IMU23=REAUX*IMG23+IMAUX*REG23
          REU33=REAUX*REG33-IMAUX*IMG33
          IMU33=REAUX*IMG33+IMAUX*REG33
          WRITE(LU2,'(19(1F16.6,1X))')
     *          F,REU11,IMU11,REU21,IMU21,REU31,IMU31,
     *            REU12,IMU12,REU22,IMU22,REU32,IMU32,
     *            REU13,IMU13,REU23,IMU23,REU33,IMU33
          I2=I2+18
  14    CONTINUE
      GOTO 10
C
  20  CONTINUE
C
      CLOSE(LU1)
      CLOSE(LU2)
      WRITE(*,'(A)') '+GREENTC: Done.                 '
      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
C
C=======================================================================
C