C 
C Program GREEN converting the unformatted output of program CRT into a
C formatted file containing the ray-theory elastodynamic Green function.
C Program GREEN can consider coupling ray-theory in weakly anisotropic
C models (without structural interfaces in this version).
C
C Version: 5.20
C Date: 1998, November 11
C
C Coded by: Petr Bulant & Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mails: bulant@seis.karlov.mff.cuni.cz
C              klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'SEP'/
C     'SEP'...Name of the file with input parameters.
C             Description of file SEP
C             If blank, default values of the corresponding data are
C             considered.
C Default: 'SEP'=' '.
C
C                                                     
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file
C     'sep.for'.
C Names of optional input files (related to ray tracing):
C     CRTOUT='string'...File with the names of the output files of
C             program CRT.
C             If blank, default names of the output files of program CRT
C             are considered.
C             Description of file CRT
C             Default: CRTOUT=' ' (usually sufficient for the first
C               elementary wave)
C     REC='string'... If non-blank, the name of the file with the names
C             and coordinates of the receiver points.  The names are
C             used within the strings describing the points of two-point
C             rays and the coordinates for the linear interpolation of
C             travel times from ray endpoints to the receivers.
C             If blank, the two-point rays are denoted by the receiver
C             index and the interpolation is not performed.
C             Description of file REC
C             Default: REC=' '
C     SRC='string'... If non-blank, the name of the file with the name
C             and coordinates of the source point.  The name is used
C             within the strings describing the two-point rays and the
C             coordinates for the linear interpolation of travel times
C             from ray initial points to the source.
C             If blank, the interpolation is not performed.
C             Description of file SRC
C             Default: SRC=' ' (mostly sufficient)
C     MODEL='string'... Name of the input formatted file with the input
C             data for the anisotropic model.  If blank, the model is
C             assumed isotropic and the data for the model are not read.
C             Description of MODEL
C             Default: MODEL=' ' (sufficient in isotropic models)
C The name of the output file:
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 Data describing the frequency domain common with program 'ss.for':
C     These data are used only if MODEL is not blank and the model is
C             anisotropic.
C     DT=real... Time step to digitize the source time function and
C             seismograms.
C             Default: DT=1.
C     NFFT=integer... Number of the time samples for the fast Fourier
C             transform.  Will be used to convert the time step to the
C             frequency step.
C             NFFT must be an integer power of 2.
C             Default: NFFT=1
C     FMIN=real... The lowest frequency of the cosine band-pass filter.
C             Default: FMIN=0 (mostly sufficient)
C     FMAX=real... The highest frequency of the cosine band-pass filter.
C             Default: FMAX=0.5/DT (Nyquist frequency, mostly
C               sufficient)
C Data describing the frequency domain specific to this program:
C     These data are used only if MODEL is not blank and the model is
C             anisotropic.
C     DF...   Frequency step.
C             Default: DF=1/(NPTS*DT) (recommended)
C     OF...   The elementary Green functions are calculated for NF
C             frequencies OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF.
C             Default: OF=DF*NINT(FMIN/DF) (i.e. FMIN rounded to the
C               nearest multiple of the frequency step, recommended)
C     NF...   Number of frequencies.
C             Default: NF=NINT((FMAX-OF)/DF)+1 (recommended)
C     Hint: Do not specify DF, OF and NF.  Use the data for program
C     'ss.for' instead.
C Example of file 'SEP'
C
C Input formatted files 'MODEL', 'REC' and 'SRC', if specified, must
C     correspond to the receiver and source files used during Complete
C     Ray Tracing.
C
C Input formatted file 'MODEL':
C     See the description within source code file 'model.for'.
C     Description of file MODEL
C
C                                                     
C Input formatted file 'CRT':
C     The data consist of character strings read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'CRT-R','CRT-S','CRT-I',/
C     'CRT-R'.. Name of the input unformatted file with the quantities
C             stored along rays (see C.R.T.5.5.1).  The file is not read
C             if SEP is blank or if MODEL=' ' in input data SEP or if
C             the model is isotropic.
C             Description of file CRT-R
C     'CRT-S'.. File with the quantities stored at the points of
C             intersections of rays with the specified reference surface
C             (see C.R.T.5.5.2).  Only two-point rays incident to the
C             vicinities of the receivers situated at the reference
C             surface are considered.  The output elementary Green
C             functions are then evaluated at the receivers.
C             Description of file CRT-S
C     'CRT-I'.. File with the quantities at the initial points of rays,
C             corresponding to files 'CRT-R' and 'CRT-S'
C             (see C.R.T.6.1).
C             Description of file CRT-I
C Default: 'CRT-R'='r01.out', 'CRT-S'='s01.out', 'CRT-I'='s01i.out'.
C
C Input unformatted file 'CRT-R':
C     See the description within source code file 'writ.for'.
C     Description of file CRT-R
C
C Input unformatted file 'CRT-S':
C     See the description within source code file 'writ.for'.
C     Description of file CRT-S
C
C Input unformatted file 'CRT-I':
C     See the description within source code file 'writ.for'.
C     Description of file CRT-I
C
C                                                   
C Output formatted file 'GREEN':
C (1) / (a slash).
C (2) For each two-point ray (2.1):
C (2.1) 'R','S',(GREEN(I),I=1,NRAYPT),/
C             Here NRAYPT=32 except for S waves in anisotropic model,
C             where NRAYPT=14+18*NF
C     'R'...  String in apostrophes describing the receiver.
C     'S'...  String in apostrophes describing the source.
C     GREEN(1)... Travel time between the receiver and the source.
C     GREEN(2)... Imaginary part of the complex-valued travel time
C             between receiver and source due to attenuation.
C     GREEN(3:8)... Coordinates of the receiver and coordinates of the
C             source.
C     GREEN(9:14)... Derivatives of the travel time with respect to the
C             coordinates of the receiver and coordinates of the source.
C     GREEN(15:32)... 1000000 times enlarged amplitude of the Green
C             function: contravariant components of the complex-valued
C             3*3 matrix Gij in model coordinates, where the first
C             subscript corresponds to the receiver and the second
C             subscript corresponds to the source.  The components are
C             ordered as
C             Re(G11),Im(G11),Re(G21),Im(G21),Re(G31),Im(G31),
C             Re(G12),Im(G12),Re(G22),Im(G22),Re(G32),Im(G32),
C             Re(G13),Im(G13),Re(G23),Im(G23),Re(G33),Im(G33).
C     GREEN(33:*)... Analogy of GREEN(15:32) for other frequencies.
C             Written just for frequency-dependent elementary Green
C             functions.
C     /...    An obligatory slash at the end of line.  It may also stand
C             for zeros at the end of line.
C (3) / (a slash).
C File form GREEN is, to some extent, an extension of file form
C Travel Times.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL AP00,AP21,TXT1,TXT2,REC,SRC,TTSORT,FORM2
C     AP00,AP21... File 'ap.for'.
C     AP03,AP03A... File 'ap.for' (called by AP21,AP03).
C     TXT1,TXT2,REC,SRC... File 'crtout.for'.
C     LENGTH..File 'length.for' (called by TXT2).
C     TTSORT..File 'ttsort.for'.
C     INDEXX..File 'indexx.for' (called by TTSORT).
C     FORM2...File 'forms.for'.
C     FORM1...File 'forms.for' (called by FORM2).
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C-----------------------------------------------------------------------
C
C     Allocation of working arrays:
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C     Auxiliary storage locations:
      INTEGER LU1,LU2,LU3,LU4
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4)
      CHARACTER*80 FILSEP,FILMOD,FILCRT,FILREC,FILSRC,FGREEN
      CHARACTER*80 FILE1,FILE2,FILE3
      CHARACTER*260 FORMAT
      CHARACTER*80 RAYTXT
      INTEGER MQ,MTPR,NTPR,IGREEN,I,K
      INTEGER NFFT,NF,NQ
      REAL DT,FMIN,FMAX,DF,OF
C     MQ...   Maximum number of quantities stored for a two-point ray:
C             IRAM(MQ*NTPR+1)... Number of quantities stored for the
C               (NTPR+1)-th two-point ray.
C             IRAM(MQ*NTPR+2)... Index of the receiver.
C             RAM(MQ*NTPR+3) to RAM(MQ*NTPR+16) ... Travel time,
C               coordinates, gradient of travel time.
C             RAM(MQ*NTPR+17) to RAM(MQ*NTPR+IRAM(MQ*NTPR+1))...
C               Amplitudes.
C     MTPR... Maximum number of two-point rays.
C             Note: MTPR=1 indicates no sorting according to receivers.
C     NTPR... Number of two-point rays.
C
C-----------------------------------------------------------------------
C
C     Opening input and output files:
C
C     Main input data:
      FILSEP=' '
      WRITE(*,'(A)') '+GREEN: Enter input filename: '
      READ(*,*) FILSEP
      WRITE(*,'(A)') '+GREEN: Working...            '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Name of the model file
      CALL RSEP3T('MODEL',FILMOD,' ')
C     Name of the receiver file
      CALL RSEP3T('REC',FILREC,' ')
C     Name of the source file
      CALL RSEP3T('SRC',FILSRC,' ')
C     Name of the output file
      CALL RSEP3T('GREEN',FGREEN,'green.out')
C
C     Reading the data for a possibly anisotropic model:
      IF(FILMOD.NE.' ') THEN
        OPEN(LU1,FILE=FILMOD,STATUS='OLD')
        CALL MODEL1(LU1)
        CLOSE(LU1)
C       Checking whether the model is anisotropic
        CALL PARM4(I)
        IF(I.GT.0) THEN
C         Isotropic model
          FILMOD=' '
        END IF
      END IF
      IF(FILMOD.NE.' ') THEN
        IF(NSRFC().NE.0) THEN
C         GREEN-06
          CALL ERROR('GREEN-06: Anisotropic model with interfaces')
C         This version of program GREEN cannot consider anisotropic
C         models with structural interfaces.
        END IF
      END IF
C
C     Reading the source and receiver names and preparing them for TXT2:
      CALL TXT1(LU2,FILSRC,FILREC)
C
C     Reading output filenames of the CRT program:
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      FILE1='r01.out'
      FILE2='s01.out'
      FILE3='s01i.out'
      IF(FILCRT.NE.' ') THEN
        OPEN(LU1,FILE=FILE2,STATUS='OLD')
        READ(LU1,*) FILE1,FILE2,FILE3
        CLOSE(LU1)
      END IF
C
C     Opening output files of the CRT program for input:
      IF(FILMOD.NE.' ') THEN
        IF(FILE1.EQ.' ') THEN
C         GREEN-01
          CALL ERROR('GREEN-01: No input file with rays')
C         Calculation of the coupling ray theory Green function
C         in the anisotropic model, specified by input parameter MODEL,
C         requires the two-point rays stored by program 'crt.for' into
C         the file which name is taken from the file given by parameter
C         CRTOUT.
        END IF
        OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD',ERR=2)
      END IF
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD',ERR=3)
      OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD',ERR=4)
      GO TO 9
C       Error messages:
    2   CONTINUE
C       GREEN-02
        CALL ERROR('GREEN-02: Cannot open the input file with rays')
C       The default name of the file is 'r01.out'.
    3   CONTINUE
C       GREEN-03
        CALL ERROR('GREEN-03: Cannot open the input file with points')
C       The default name of the file is 's01.out'.
    4   CONTINUE
C       GREEN-04
        CALL ERROR
     *      ('GREEN-04: Cannot open the input file with initial points')
C       The default name of the file is 's01i.out'.
    9 CONTINUE
C
C     Initializing the output file:
      OPEN(LU4,FILE=FGREEN)
      WRITE(LU4,'(A)') '/'
C
C     Number of frequencies, maximum number of values per two-point ray:
      IF(FILMOD.EQ.' ') THEN
        NF=1
        MQ=34
      ELSE
        CALL RSEP3I('NFFT',NFFT,1)
        CALL RSEP3R('DT  ',DT  ,1.)
        CALL RSEP3R('FMIN',FMIN,0.)
        CALL RSEP3R('FMAX',FMAX,0.5/DT)
        CALL RSEP3R('DF  ',DF  ,1./(FLOAT(NFFT)*DT))
        CALL RSEP3R('OF  ',OF  ,DF*NINT(FMIN/DF))
        CALL RSEP3I('NF  ',NF  ,NINT((FMAX-OF)/DF)+1)
        MQ=16+18*NF
      END IF
C
C     Maximum number of two-point rays:
      MTPR=MRAM/(MQ+2)
C     Note: MQ is the maximum number of quantities per a two-point ray
C           (including two integers), and the space for 2 working arrays
C           (required for sorting) is reserved at the end of array RAM.
C     Note: If MTPR=1, the rays are not sorted according to receivers.
C
C.......................................................................
C
C     Loop for the points of rays:
      NTPR=0
   10 CONTINUE
C       Checking free memory for the next two-point ray:
        IF (NTPR.GE.MTPR) THEN
C         GREEN-05
          CALL ERROR('GREEN-05: Too many Elementary Green functions')
C         Elementary Green functions require 36 storage locations for
C         each two-point ray in an isotropic model and 18+18*NF storage
C         locations for each two-point ray in an anisotropic model.
C         The Green functions are stored in array RAM declared in
C         include file ram.inc.
        END IF
C
C       Reading the results of the complete ray tracing
        IGREEN=MQ*NTPR
        IF(FILMOD.EQ.' ') THEN
C         Isotropic model
   20     CONTINUE
            CALL AP00(0,LU2,LU3)
            IF(IWAVE.LT.1) THEN
C             End of rays
              GO TO 60
            END IF
          IF (IREC.LE.0) GO TO 20
C         Determination of the Green function
          CALL AP21(RAM(IGREEN+3))
          NQ=32
        ELSE
C         Anisotropic model
          CALL WAN(LU1,LU2,LU3,RAM(IGREEN+3),MQ-2,RAM(IGREEN+MQ+1),
     *                              IRAM(IGREEN+MQ+1),MRAM-IGREEN-MQ,NQ)
          IF(NQ.LE.0) THEN
C           End of rays
            GO TO 60
          END IF
        END IF
C
        IRAM(IGREEN+1)=NQ+2
        IRAM(IGREEN+2)=IREC
C       Linear Taylor expansion of travel time
        CALL RECSRC(FILREC,FILSRC,RAM(IGREEN+3))
C       Rescaling the amplitudes of the Green function
        DO 30 I=IGREEN+17,IGREEN+2+NQ
          RAM(I)=1000000.*RAM(I)
   30   CONTINUE
C
C       Writing unsorted two-point rays if MTPR=1:
        NTPR=NTPR+1
        IF(MTPR.EQ.1) THEN
          CALL WGREEN(LU4,IRAM(1)-2,IRAM(2),RAM(3))
          NTPR=0
        END IF
      GO TO 10
   60 CONTINUE
C
C.......................................................................
C
C     Sorting and writing two-point rays:
      IF(MTPR.GT.1) THEN
C       Collecting the indices of the receivers for sorting:
        DO 70 K=1,NTPR
          IRAM(MQ*NTPR+K)=IRAM(MQ*(K-1)+2)
   70   CONTINUE
C
C       Sorting elementary Green functions according to the receivers:
        CALL TTSORT(MQ,NTPR,3,RAM,IRAM(MQ*NTPR+1),
     *                             RAM(MQ*NTPR+1),IRAM((MQ+1)*NTPR+1))
C
C       Writing elementary Green functions:
        DO 80 K=(MQ+1)*NTPR+1,(MQ+2)*NTPR
          IGREEN=MQ*(IRAM(K)-1)
          CALL WGREEN(LU4,IRAM(IGREEN+1)-2,IRAM(IGREEN+2),RAM(IGREEN+3))
   80   CONTINUE
      END IF
C
C.......................................................................
C
      WRITE(LU4,'(A)') '/'
      CLOSE(LU4)
      CLOSE(LU3)
      CLOSE(LU2)
      IF(FILMOD.NE.' ') THEN
        CLOSE(LU1)
      END IF
      WRITE(*,'(A)') '+GREEN: Done.                                    '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE RECSRC(FILREC,FILSRC,GREEN)
      CHARACTER*(*) FILREC,FILSRC
      REAL GREEN(14)
C
C-----------------------------------------------------------------------
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C-----------------------------------------------------------------------
C
      REAL COOR(3)
C
C.......................................................................
C
C     Linear Taylor expansion of travel time:
C     Receiver:
      IF(FILREC.NE.' ') THEN
C       Receiver:
        COOR(1)=GREEN(3)
        COOR(2)=GREEN(4)
        COOR(3)=GREEN(5)
        CALL REC(IREC,COOR(1),COOR(2),COOR(3))
        GREEN(1)=GREEN(1)+(COOR(1)-GREEN(3))*GREEN( 9)
     *                   +(COOR(2)-GREEN(4))*GREEN(10)
     *                   +(COOR(3)-GREEN(5))*GREEN(11)
        GREEN(3)=COOR(1)
        GREEN(4)=COOR(2)
        GREEN(5)=COOR(3)
      END IF
      IF(FILSRC.NE.' ') THEN
C       Source:
        COOR(1)=GREEN(6)
        COOR(2)=GREEN(7)
        COOR(3)=GREEN(8)
        CALL SRC(1,COOR(1),COOR(2),COOR(3))
        GREEN(1)=GREEN(1)+(COOR(1)-GREEN(6))*GREEN(12)
     *                   +(COOR(2)-GREEN(7))*GREEN(13)
     *                   +(COOR(3)-GREEN(8))*GREEN(14)
        GREEN(6)=COOR(1)
        GREEN(7)=COOR(2)
        GREEN(8)=COOR(3)
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WGREEN(LU4,NGREEN,IREC,GREEN)
      INTEGER LU4,NGREEN,IREC
      REAL GREEN(NGREEN)
C
C-----------------------------------------------------------------------
C
      CHARACTER*260 FORMAT
      CHARACTER*80 FILREC,RAYTXT
      INTEGER LENTXT,I,L
C
C.......................................................................
C
C     Text strings:
      CALL TXT2(0,1,1,1,IREC,LENTXT,RAYTXT)
      L=INDEX(RAYTXT(1:LENTXT),'''')
      L=INDEX(RAYTXT(L+1:LENTXT),'''')+L
      CALL RSEP3T('REC',FILREC,' ')
      IF(FILREC.EQ.' ') THEN
C       Shortening the receiver string part from 8 to 6 characters
        LENTXT=LENTXT-2
        DO 71 I=L+4,LENTXT
          RAYTXT(I:I)=RAYTXT(I+2:I+2)
   71   CONTINUE
      END IF
C
C     Writing:
      FORMAT(1:4)='(4A,'
      IF(NGREEN.LE.32) THEN
        CALL FORM2(32,GREEN,GREEN,FORMAT(5:260))
        WRITE(LU4,FORMAT) RAYTXT(L+2:LENTXT),' ',RAYTXT(1:L),
     *                      (' ',GREEN(I),I=1,NGREEN),' /'
      ELSE
        CALL FORM2(14,GREEN,GREEN,FORMAT(5:260))
        WRITE(LU4,FORMAT) RAYTXT(L+2:LENTXT),' ',RAYTXT(1:L),
     *                      (' ',GREEN(I),I=1,14)
        L=15
        DO 80 L=L,NGREEN-18,18
          FORMAT(1:4)='(1A,'
          CALL FORM2(18,GREEN(L),GREEN(L),FORMAT(5:260))
          WRITE(LU4,FORMAT) (' ',GREEN(I),I=L,L+17)
   80   CONTINUE
        FORMAT(1:4)='(1A,'
        CALL FORM2(NGREEN-L+1,GREEN(L),GREEN(L),FORMAT(5:260))
        WRITE(LU4,FORMAT) (' ',GREEN(I),I=L,NGREEN),' /'
      END IF
C
      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 'indexx.for'
C     indexx.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'crtout.for'
C     crtout.for
      INCLUDE 'ttsort.for'
C     ttsort.for
      INCLUDE 'wan.for'
C     wan.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