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.
C
C Version: 5.50
C Date: 2001, February 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     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 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 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, 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.
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, 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 QIERR is exceeded, the
C             warning message is generated.
C             Default: ERRQI=ERRWAN
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 modifications 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=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 modifications 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             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 Example of file 'SEP'
C
C                                                  
C Input formatted file CRTOUT:
C     For general description of file CRTOUT refer to the file
C     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
C Subroutines and external functions required:
      EXTERNAL KOOR,NSRFC,AP00,AP21,TXT1,TXT2,REC,SRC,TTSORT,FORM2
      INTEGER  KOOR,NSRFC
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 LU0,LU1,LU2,LU3,LU4
      PARAMETER (LU0=1,LU1=2,LU2=3,LU3=4,LU4=5)
      CHARACTER*80 FILSEP,FILMOD,FILCRT,FILREC,FILSRC,FGREEN
      CHARACTER*80 FILE1,FILE2,FILE3,FILO1,FILO2,FILO3
      CHARACTER*260 FORMAT
      CHARACTER*80 RAYTXT,TEXT
      INTEGER MQ,MTPR,NTPR,IW,IGREEN,I,K
      INTEGER NFFT,NF,NQ
      REAL DT,FMIN,FMAX,DF,OF,ERRWAN,ERRQI,ERR(2)
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     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
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:
      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
        CALL RSEP3R('ERRWAN',ERRWAN,0.001)
        CALL RSEP3R('ERRQI' ,ERRQI ,ERRWAN)
        ERR(1)=0.
        ERR(2)=0.
      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)
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 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,ERR)
          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
      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 for sorting:
        DO 91 K=1,NTPR
          IRAM(MQ*NTPR+K)=IRAM(MQ*(K-1)+2)
   91   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 92 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))
   92   CONTINUE
      END IF
C
C.......................................................................
C
      WRITE(LU4,'(A)') '/'
      CLOSE(LU4)
      CLOSE(LU3)
      CLOSE(LU2)
      IF(FILMOD.NE.' ') THEN
        CLOSE(LU1)
      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(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