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.  For S waves, second-order perturbations of anisotropic travel
C times may be optionally calculated.
C The Green function can be calculated only for a point source.
C
C Version: 7.40
C Date: 2017, May 19
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     http://sw3d.cz/staff/bulant.htm
C     http://sw3d.cz/staff/klimes.htm
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     No default, 'SEP' must be specified and cannot be blank.
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 '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 are considered.
C             Description of file CRTOUT.
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 names
C             and coordinates of the source points.  The names are used
C             within the strings describing the two-point rays and the
C             coordinates are used for the linear interpolation of
C             travel times 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 to be used for the coupling
C             ray theory Green function for S waves.  If this parameter
C             is specified and the model is not isotropic, the coupling
C             ray theory is applied.  If this parameter is blank, the
C             model is assumed isotropic, the data for the model are not
C             read, and the coupling ray theory is not applied to S
C             waves.  Similarly, if this parameter is specified and the
C             model is isotropic, the coupling ray theory is not
C             applied to S waves.
C             This parameter has no effect to P waves.
C             Description of file 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 Prevailing-frequency approximation of the coupling ray theory:
C     SINGLF=positive real... Prevailing frequency for the
C             prevailing-frequency approximation of the coupling ray
C             theory.  This parameter is used only if MODEL is not blank
C             and the model is anisotropic.
C             Default: SINGLF=0. (no prevailing-frequency approximation)
C Optional data describing sorting of the Green functions:
C     ISORT=integer... Specifies primary sorting criterion for the
C             output files with the Green functions:
C             ISORT=0 ... Sorting according to the receiver index
C             ISORT=1 ... Sorting according to the source index
C             Default: ISORT=0 (sorting according to the receiver index)
C     ISORT1=integer... Specifies secondary sorting criterion for the
C             output files with the Green functions. The Green functions
C             are first sorted according to receiver index or source
C             index, as specified by ISORT, and then according to the
C             quantity specified by ISORT1:
C             ISORT1=0 ... receiver index
C             ISORT1=1 ... source index
C             ISORT1=2 ... travel time
C             ISORT1=3 ... imaginary part of complex travel time
C             ISORT1=4 to 6 ... receiver coordinates
C             ISORT1=7 to 9 ... source coordinates
C             ISORT1=10 to 15 .. derivatives of travel time with respect
C               to the coordinates of the receiver and of the source
C             ISORT1=16 to 33 .. amplitude of the Green function
C             Default: ISORT1=2 (sorting according to the travel time)
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, and SINGLF=0. (no prevailing-frequency
C             approximation), i.e., when the Green function is frequency
C             dependent.  These data common with 'ss.for' are used to
C             generate the optimum default values of input parameters
C             DF, OF and NF compatible with data for 'ss.for'.
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, and SINGLF=0. (no prevailing-frequency
C             approximation), i.e., when the Green function is frequency
C             dependent.
C     DF=real... Frequency step.
C             Default: DF=1/(NFFT*DT) (recommended)
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=DF*NINT(FMIN/DF) (i.e. FMIN rounded to the
C               nearest multiple of the frequency step, recommended)
C     NF=integer...   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 and leave the default values of DF, OF and NF
C     if you are going to generate synthetic seismograms by 'ss.for'.
C     Specify explicitly DF, OF and NF only if you generate the Green
C     function for another purpose than the generation of synthetic
C     seismograms by programs 'greenss.for' and 'ss.for'.
C Auxiliary parameters for anisotropic models, need not be specified:
C     ERRWAN=real... Maximum error in propagator matrix along whole ray.
C             Default: ERRWAN=0.001
C     ERRQI=real... Maximum error of the quasi-isotropic projection of
C             the polarization vectors.  If ERRQI is exceeded, the
C             warning message is generated.
C             This parameter has no meaning for QIPV=0.
C             Default: ERRQI=ERRWAN
C     QIPV=integer... Optional quasi-isotropic projection of
C             the polarization vectors.
C             QIPV=0: No modification of the polarization vectors of the
C               coupling ray theory.  Eigenvectors of the Christoffel
C               matrix are used as the polarization vectors in
C               anisotropic models (default).
C             QIPV=1: Quasi-isotropic projection of the polarization
C               vectors onto the basis vectors of the ray-centred
C               coordinate system.  This option corresponds to versions
C               5.20 to 5.50 of this program, and may also be useful for
C               comparison of the results with package ANRAY.
C             Both options should yield equivalent results if QICHM=1.
C             Default: QIPV=0
C     QICHM=integer... Optional quasi-isotropic modification of the
C             Christoffel matrix.
C             QICHM=0: No modification of the Christoffel matrix
C               (default).
C             QICHM=1: Quasi-isotropic modification of the Christoffel
C               matrix.  Mixed components Gamma_12, Gamma_21, Gamma_13
C               and Gamma_31 of the Christoffel matrix in the
C               ray-centred coordinates are put to zero, in order to
C               keep the reference isotropic polarization vectors
C               unperturbed.  This modification is designed to compare
C               the results with package ANRAY.
C             Default: QICHM=0
C     QITT=integer... Optional quasi-isotropic approximation of the
C             average anisotropic travel time.
C             QITT=0: No modification of the average anisotropic travel
C               time calculated along the reference ray.  Anisotropic
C               travel times are linearized with respect to the
C               Christoffel matrix powered to -1/2, which represents
C               the slowness (default).
C               QITT=0 corresponds to the minus first degree
C               of the homogeneous Hamiltonian used for the calculation
C               of perturbations along the reference ray.
C             QITT=1: Quasi-isotropic modification of the average
C               anisotropic travel time calculated along the reference
C               ray.  The anisotropic travel times are linearized with
C               respect to the Christoffel matrix, which represents
C               the squared velocity.  This modification is designed to
C               compare the results with package ANRAY, for the default
C               value of parameter INEWB of program 'weakan.for'.
C               QITT=1 corresponds to INEWB=0, QITT=0 to INEWB=1.
C               QITT=1 corresponds to the second degree
C               of the homogeneous Hamiltonian used for the calculation
C               of perturbations along the reference ray.
C             Default: QITT=0
C     QIDT=integer... Optional quasi-isotropic approximation of the
C             difference between the anisotropic travel times.
C             QIDT=0: No modification of the difference between the
C               anisotropic travel times (default).
C             QIDT=1: Quasi-isotropic modification of the difference
C               between the anisotropic travel times.  The anisotropic
C               travel times are linearized with respect to the
C               Christoffel matrix instead of its inverse square root.
C               This modifications is designed to compare the results
C               with package ANRAY, for the default value of parameter
C               INEWB of program 'weakan.for'.
C             Default: QIDT=0
C     EIGLST='string' ... Name of the optional output file containing
C             eigenvectors of Christoffel matrix.
C             Specification of EIGLST does not influence the output
C             file with the Green function.
C             Description of file EIGLST.
C             Default: EIGLST=' ' (the file is not generated)
C Parameters for optional second-order perturbations of anisotropic
C travel times of S waves. The perturbations may be calculated only
C along the isotropic or anisotropic rays which do not cross any
C interface.
C The perturbations can not be calculated for P waves.
C Furthermore, in more complex models in the vicinity of caustics,
C the second-order perturbations may be infinitely large and should
C be avoided.
C     QIRAY=integer ... Optional second-order perturbations
C             of anisotropic travel times of S waves for the
C             calculation of the output Green function.
C             QIRAY=0: No second-order perturbations (default).
C             QIRAY=1: Second-order perturbations of anisotropic
C               travel times of S waves are calculated and are taken
C               into account when calculating the output Green
C               function.
C             Default: QIRAY=0
C     QILST='string' ... Name of the optional output file
C             containing second-order perturbations of anisotropic
C             travel times of S waves.
C             Specification of QILST does not influence the output
C             file with the Green function, the second-order
C             perturbations are calculated and written to file QILST.
C             Description of file QILST.
C             Default: QILST=' ' (the file is not generated)
C Parameters used for ray tracing. The values of the
C parameters must be the same as used for ray tracing:
C     CRTANI=real... Switch to anisotropic ray tracing.
C             CRTANI=0... Isotropic ray tracing.
C             CRTANI.NE.0... Anisotropic ray tracing.
C             See the description in the file crtin.for.
C             Default: CRTANI=0.
C     DEGREE=real... Degree of the homogeneous Hamiltonian
C             arithmetically averaged during common S-wave ray tracing.
C             See the description in the file crtin.for.
C             Default: DEGREE=-1.
C     KSWAVE=integer... Selection of the anisotropic-ray-theory
C             S-wave rays instead of common S-wave rays. Affects only
C             S waves.
C             See the description in the file crtin.for.
C Optional parameters specifying the form of the real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
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 Example of file SEP.
C
C                                                  
C Input formatted file CRTOUT:
C     For general description of file CRTOUT refer to the file writ.for.
C Description specific to this program:
C     'CRT-R'... The file is not read if MODEL=' ' or if the model
C             is isotropic.  If the file is read, only two-point rays
C             stored in both files 'CRT-R' and 'CRT-S' are selected.
C     'CRT-S'... Only two-point rays (rays incident to the vicinities of
C             the receivers situated at the reference surface) are read
C             from file 'CRT-S'.  The output elementary Green functions
C             are then evaluated at the receivers.
C     'CRT-I'... If file 'CRT-R' is not read, file 'CRT-I' must contain
C             at least initial points of all rays contained in file
C             'CRT-S'.  If file 'CRT-R' is read, file 'CRT-I' must
C             contain at least the initial points of all rays common
C             to files 'CRT-R' and 'CRT-S'.
C     'CRT-T'... Not used.
C     Only the first data line is read in the current version.
C     Defaults for the first data line:
C             'CRT-R'='r01.out', 'CRT-S'='s01.out', 'CRT-I'='s01i.out'
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 Output formatted file QILST with S-wave travel-time perturbations:
C (1) / (a slash).
C (2) For each two-point ray one line (2.1):
C     For isotropic two-point ray:
C (2.1) IREC,TI,DTLI1,DTLI2,DTQI1,DTQI2,T1,T1,T2,T2,
C       DTLIA,DTQIA,TA,TA,DTQA1,DTQA2,DTQ12 /
C     IREC... Number of the receiver corresponding to the two-point ray.
C     TI...   Isotropic travel time of S wave.
C     DTLI1,DTLI2... First-order perturbations of travel time,
C             from the isotropic S-wave reference ray to the
C             anisotropic ray theory S1 and S2 rays.
C     DTQI1,DTQI2... Second-order perturbations of travel time,
C             from the isotropic S-wave reference ray to the
C             anisotropic ray theory S1 and S2 rays.
C     T1,T2...Approximation of the anisotropic travel times by the
C             second-order perturbation expansion of travel time,
C             from the isotropic S-wave reference ray to the
C             anisotropic ray theory S1 and S2 rays.
C     DTLIA...First-order perturbations of travel time,
C             from the isotropic S-wave reference ray to the
C             anisotropic S-wave common reference ray.
C     DTQIA...Second-order perturbations of travel time,
C             from the isotropic S-wave reference ray to the
C             anisotropic S-wave common reference ray.
C     TA...   Approximation of the anisotropic S-wave common
C             reference ray travel time by the second-order
C             perturbation expansion of travel time,
C             from the isotropic S-wave reference ray to the
C             anisotropic S-wave common reference ray.
C     DTQA1,DTQA2... Estimated second-order perturbations of travel
C             time, from the anisotropic S-wave common reference ray
C             to the anisotropic ray theory S1 and S2 rays.
C     DTQ12...Estimated second-order perturbations of travel
C             time, from the anisotropic ray theory S1 ray
C             to the anisotropic ray theory S2 ray.
C     /...    A slash at the end of line.
C     For common anisotropic two-point ray:
C (2.1) IREC,TA,DTLA1,DTLA2,DTQA1,DTQA2,T1,T1,T2,T2,
C       DTLAI,DTQAI,TI,TI,DTQI1,DTQI2,DTQ12 /
C     IREC... Number of the receiver corresponding to the two-point ray.
C     TA...   Common anisotropic travel time of S wave.
C     DTLA1,DTLA2... First-order perturbations of travel time,
C             from the anisotropic common S-wave reference ray to the
C             anisotropic ray theory S1 and S2 rays.
C     DTQA1,DTQA2... Second-order perturbations of travel time,
C             from the anisotropic common S-wave reference ray to the
C             anisotropic ray theory S1 and S2 rays.
C     T1,T2...Approximation of the anisotropic travel times by the
C             second-order perturbation expansion of travel time,
C             from the anisotropic common S-wave reference ray to the
C             anisotropic ray theory S1 and S2 rays.
C     DTLAI...First-order perturbations of travel time,
C             from the anisotropic common S-wave reference ray to the
C             isotropic ray theory S wave ray.
C     DTQAI...Second-order perturbations of travel time,
C             from the anisotropic common S-wave reference ray to the
C             isotropic ray theory S wave ray.
C     TI...   Approximation of the isotropic S wave ray travel time
C             by the second-order perturbation expansion of travel time,
C             from the anisotropic common S-wave reference ray to the
C             isotropic ray theory S wave ray.
C     DTQI1,DTQI2... Estimated second-order perturbations of travel
C             time, from the isotropic S wave ray
C             to the anisotropic ray theory S1 and S2 rays.
C     DTQ12...Estimated second-order perturbations of travel
C             time, from the anisotropic ray theory S1 ray
C             to the anisotropic ray theory S2 ray.
C     /...    A slash at the end of line.
C (3) / (a slash).
C
C                                                  
C Output formatted file EIGLST with eigenvectors of Christoffel matrix:
C (1) / (a slash).
C (2) For each two-point ray data (2.1) and (2.2):
C (2.1) IREC,/
C     IREC... Number of receiver corresponding to the two-point ray.
C     /...    A slash at the end of line.
C (2.2) For each point along the two-point ray data (2.2.1)
C (2.2.1) IPTS,TT,G11,G12,G13,G21,G22,G23,E1,E2,/
C
C     IPTS... Index of the point, zero for points added to the ray
C             by interpolation.
C     TT...   Travel time in the point.
C     G11,G12,G13,G21,G22,G23... Eigenvectors of qS1 and qS2 wave.
C     E1...   Eigenvalue of Christoffel matrix
C             corresponding to qS1 (faster) wave.
C     E2...   Eigenvalue of Christoffel matrix
C             corresponding to qS2 (slower) wave.
C     /...    A slash at the end of line.
C (3) / (a slash).
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL KOOR,NSRFC,AP00,AP21,TXT1,TXT2,REC,SRC,TTSORT,
     *         FORM2,LENGTH,WAN
      INTEGER  KOOR,NSRFC,LENGTH
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     FORM2...File 'forms.for'.
C     FORM1...File 'forms.for' (called by FORM2).
C     LENGTH..File 'length.for'.
C     TTSORT..File 'ttsort.for'.
C     INDEXX..File 'indexx.for' (called by TTSORT).
C     WAN...  File 'wan.for'.
C     WAPTS,WAREAD,WACHRI,WAMAT,WACHAN,WAPROJ,WAPERT,WASUM,WASUM4,WASUM5
C             ... File 'wan.for' (called by WAN).
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 LU0,LU1,LU2,LU3,LU4,LUQI,LUEIG
      PARAMETER (LU0=1,LU1=2,LU2=3,LU3=4,LU4=5)
      CHARACTER*80 FILSEP,FILMOD,FILCRT,FILREC,FILSRC,FGREEN,FQILST,FEIG
      CHARACTER*80 FILE1,FILE2,FILE3,FILO1,FILO2,FILO3
      CHARACTER*260 FORMAT
      CHARACTER*80 RAYTXT,TEXT
      INTEGER MQ,MW,MTPR,NTPR,IW,IGREEN,I,K
      INTEGER NFFT,NF,NQ,KQIPV,ISORT,ISORT1,ISINGL,NSINGL
      REAL DT,FMIN,FMAX,DF,OF,ERRWAN,ERRQI,ERR(2),SINGLF
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             IRAM(MQ*NTPR+3)... Index of the source.
C             RAM(MQ*NTPR+4) to RAM(MQ*NTPR+17)... Travel time,
C               coordinates, gradient of travel time.
C             RAM(MQ*NTPR+18) to RAM(MQ*NTPR+IRAM(MQ*NTPR+1))...
C               Amplitudes.
C     MW...   MQ for 'wan.for'
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     IW...   Loop variable: index of the elementary wave.
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GREEN: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       GREEN-06
        CALL ERROR('GREEN-06: 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)') '+GREEN: Working ...           '
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     Parameters describing sorting of the Green functions
      CALL RSEP3I('ISORT', ISORT, 0)
      CALL RSEP3I('ISORT1',ISORT1,2)
      ISORT1=ISORT1+2
      IF (.NOT.(ISORT.EQ.0.OR.ISORT.EQ.1).OR.
     *    .NOT.(ISORT1.GE.0.AND.ISORT1.LE.33)) THEN
C       GREEN-09
        CALL ERROR('GREEN-09: Wrong value of ISORT or ISORT1')
C       Refer to the description of input parameters ISORT an ISORT1
C       above.
      ENDIF
C
C     Reading the data for a possibly anisotropic model:
      IF(FILMOD.NE.' ') THEN
        OPEN(LU0,FILE=FILMOD,STATUS='OLD')
        CALL MODEL1(LU0)
        CLOSE(LU0)
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(KOOR().NE.0) THEN
C         GREEN-07
          CALL ERROR
     *        ('GREEN-07: Anisotropic model in curvilinear coordinates')
C         Current version of subroutine file 'wan.for' for weakly
C         anisotropic models is coded only for Cartesian coordinates.
        END IF
        IF(NSRFC().NE.0) THEN
C         GREEN-51
          CALL WARN('GREEN-51: Anisotropic model with interfaces')
C         This version of program GREEN is not debugged for anisotropic
C         models with structural interfaces.
        END IF
      END IF
C
C     Number of frequencies, maximum number of values per two-point ray:
      CALL RSEP3R('SINGLF',SINGLF,0.)
      IF(FILMOD.EQ.' ') THEN
        CALL RSEP3I('QIPV',KQIPV,0)
        NF=1
        MQ=35
      ELSE
C       Coupling ray theory
        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=17+18*NF
        MW=MQ
        IF (SINGLF.NE.0.) THEN
C         Prevailing-frequency approximation
          MQ=35
          MW=70
        ENDIF
        CALL RSEP3R('ERRWAN',ERRWAN,0.001)
        CALL RSEP3R('ERRQI' ,ERRQI ,ERRWAN)
        ERR(1)=0.
        ERR(2)=0.
        CALL RSEP3T('QILST' ,FQILST ,' ')
        IF (FQILST.EQ.' ') THEN
          LUQI=0
        ELSE
          LUQI=6
          OPEN(LUQI,FILE=FQILST)
          WRITE(LUQI,'(A)') '/'
        ENDIF
        CALL RSEP3T('EIGLST',FEIG,' ')
        IF (FEIG.EQ.' ') THEN
          LUEIG=0
        ELSE
          LUEIG=7
          OPEN(LUEIG,FILE=FEIG)
          WRITE(LUEIG,'(A)') '/'
        ENDIF
      END IF
C
C     Reading the source and receiver names and preparing them for TXT2:
      CALL TXT1(LU0,FILSRC,FILREC)
C
C     Initializing the output file:
      OPEN(LU4,FILE=FGREEN)
      WRITE(LU4,'(A)') '/'
C
C     Maximum number of two-point rays:
      MTPR=MRAM/(MQ+2)
      IF (SINGLF.NE.0.) THEN
C       Prevailing-frequency approximation
        MTPR=MTPR-1
      ENDIF
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     Number of stored two-point rays:
      NTPR=0
C
C.......................................................................
C
C     Loop for elementary waves:
      FILO1=' '
      FILO2=' '
      FILO3=' '
      DO 70 IW=1,1
C
C.......................................................................
C
C     Reading output filenames of the CRT program:
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      IF(IW.EQ.1) THEN
        FILE1='r01.out'
        FILE2='s01.out'
        FILE3='s01i.out'
      ELSE
        FILE1=' '
        FILE2=' '
        FILE3=' '
      END IF
      IF(FILCRT.NE.' ') THEN
        IF(IW.EQ.1) THEN
          OPEN(LU0,FILE=FILCRT,STATUS='OLD')
        END IF
        READ(LU0,*,END=90) FILE1,FILE2,FILE3
      END IF
      IF(FILE1.EQ.' '.AND.FILE2.EQ.' '.AND.FILE3.EQ.' ') THEN
C       No new elementary wave
        GO TO 90
      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
        IF(FILE1.NE.FILO1) THEN
          IF(IW.GT.1) THEN
            CLOSE(LU1)
          END IF
          OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD',ERR=2)
        END IF
      END IF
      IF(FILE2.NE.FILO2) THEN
        IF(IW.GT.1) THEN
          CLOSE(LU2)
        END IF
        OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD',ERR=3)
      END IF
      IF(FILE3.NE.FILO3) THEN
        IF(IW.GT.1) THEN
          CLOSE(LU3)
        END IF
        OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD',ERR=4)
      END IF
      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.......................................................................
C
C     Loop for the points of rays:
   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 35 storage locations for
C         each two-point ray in an isotropic model and 17+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(KQIPV,RAM(IGREEN+4))
          NQ=32
        ELSE
C         Anisotropic model
          CALL WAN(LU1,LU2,LU3,RAM(IGREEN+4),MW-3,RAM(IGREEN+MW+1),
     *             IRAM(IGREEN+MW+1),MRAM-IGREEN-MW,NQ,ERR,LUQI,LUEIG)
          IF(NQ.LE.0) THEN
C           End of rays
            GO TO 60
          END IF
        END IF
C
C       Check on a point source:
        IF(YI(12).NE.0..OR.YI(16).NE.0..OR.
     *     YI(13).NE.0..OR.YI(17).NE.0.) THEN
C         GREEN-08
          CALL ERROR('GREEN-08: Non-point initial conditions')
C         The Green function can be calculated only for a point source.
C         This restriction has been introduced in order to avoid
C         incorrect synthetic seismograms calculated by a mistake.
C         You may overcome this restriction by disabling the above
C         check.
        END IF
C
        IF(SINGLF.NE.0.) THEN
C         Prevailing-frequency approximation
          NQ=32
          NSINGL=1
          DO 21 I=IGREEN+35,IGREEN+4,-1
            RAM(I+MQ)=RAM(I+32)
   21     CONTINUE
        ELSE
          NSINGL=0
        END IF
        DO 40 ISINGL=0,NSINGL
          IRAM(IGREEN+1)=NQ+3
          IRAM(IGREEN+2)=IREC
          IRAM(IGREEN+3)=ISRC
C         Linear Taylor expansion of travel time
          CALL RECSRC(FILREC,FILSRC,RAM(IGREEN+4))
C         Rescaling the amplitudes of the Green function
          DO 30 I=IGREEN+18,IGREEN+3+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(IGREEN+1)-3,IRAM(IGREEN+2),
     *                      IRAM(IGREEN+3),RAM(IGREEN+4))
            NTPR=0
          END IF
C
          IGREEN=IGREEN+MQ
   40   CONTINUE
      GO TO 10
   60 CONTINUE
C
C.......................................................................
C
      FILO1=FILE1
      FILO2=FILE2
      FILO3=FILE3
   70 CONTINUE
C     End of loop for elementary waves.
C
C.......................................................................
C
C     Sorting and writing two-point rays:
   90 CONTINUE
      IF(MTPR.GT.1) THEN
C       Collecting the indices of the receivers or sources for sorting:
        DO 91 K=1,NTPR
          IRAM(MQ*NTPR+K)=IRAM(MQ*(K-1)+2+ISORT)
   91   CONTINUE
C
C       Sorting elementary Green functions according to the receivers:
        IF (ISORT1.LE.3) THEN
          DO 93, K=1,NTPR
            RAM(MQ*(K-1)+ISORT1)=FLOAT(IRAM(MQ*(K-1)+ISORT1))
   93     CONTINUE
        ENDIF
        CALL TTSORT(MQ,NTPR,ISORT1,RAM,IRAM(MQ*NTPR+1),
     *                             RAM(MQ*NTPR+1),IRAM((MQ+1)*NTPR+1))
        IF (ISORT1.LE.3) THEN
          DO 94, K=1,NTPR
            IRAM(MQ*(K-1)+ISORT1)=NINT(RAM(MQ*(K-1)+ISORT1)+0.001)
   94     CONTINUE
        ENDIF
C
C       Writing elementary Green functions:
        DO 92 K=(MQ+1)*NTPR+1,(MQ+2)*NTPR
          IGREEN=MQ*(IRAM(K)-1)
          CALL WGREEN(LU4,IRAM(IGREEN+1)-3,IRAM(IGREEN+2),
     *                    IRAM(IGREEN+3),RAM(IGREEN+4))
   92   CONTINUE
      END IF
C
C.......................................................................
C
      WRITE(LU4,'(A)') '/'
      CLOSE(LU4)
      CLOSE(LU3)
      CLOSE(LU2)
      IF(FILMOD.NE.' ') THEN
        CLOSE(LU1)
        IF (FQILST.NE.' ') THEN
          WRITE(LUQI,'(A)') '/'
          CLOSE(LUQI)
        ENDIF
        IF (FEIG.NE.' ') THEN
          WRITE(LUEIG,'(A)') '/'
          CLOSE(LUEIG)
        ENDIF
      END IF
      IF(FILCRT.NE.' ') THEN
        CLOSE(LU0)
      END IF
      IF(FILMOD.NE.' ') THEN
        IF(ERR(1).GT.ERRQI.OR.ERR(2).GT.ERRQI) THEN
C         GREEN-52
          WRITE(TEXT,'(2(A,F8.6))')
     *                  'GREEN-52: QI projection error: source ',ERR(1),
     *                                             ', receiver ',ERR(2)
          CALL WARN(TEXT(1:LENGTH(TEXT)))
        END IF
      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(ISRC,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,ISRC,GREEN)
      INTEGER LU4,NGREEN,IREC,ISRC
      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,ISRC,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 'eigen.for'
C     eigen.for
      INCLUDE 'indexx.for'
C     indexx.for of Numerical Recipes
      INCLUDE 'trans.for'
C     trans.for
      INCLUDE 'coef.for'
C     coef.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'an.for'
C     an.for
      INCLUDE 'crtout.for'
C     crtout.for
      INCLUDE 'ttsort.for'
C     ttsort.for
      INCLUDE 'wan.for'
C     wan.for
      INCLUDE 'wanpfa.for'
C     wanpfa.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
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'hder.for'
C     hder.for
C
C=======================================================================
C