C
C Program GRDBORN to calculate the Born approximation of the wavefield
C at specified receivers
C
C Version 6.60
C Date: 2012, June 5
C
C Coded by: Libor Sachl
C     Department of Geophysics, Charles University in Prague,
C     E-mail: sachl@karel.troja.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of the input files describing the medium perturbations:
C     The files have form GRD (GRiDded data).
C     VPPER='string'... Name of the file with the grid values of the
C             perturbations of the product of the density and the square
C             of P-wave velocity.
C             Default: VPPER=' '
C     VSPER='string'... Name of the file with the grid values of the
C             perturbations of the product of the density and the square
C             of S-wave velocity.
C             Default: VSPER=' '
C     RHOPER='string'... Name of the file with the grid values of the
C             density perturbations.
C             Default: RHOPER=' '
C Name of the input file specifying the receiver Green functions:
C     BORN='string'... Name of the file containing the names of the
C             receivers and the names of the corresponding files with
C             the quantities required for the Born approximation.
C             The multivalued quantities are discretized on the grid,
C             are related to the waves arriving to the receiver from the
C             gridpoints, and are stored in the files of form
C             MGRD (Multivalued GRiDded data).
C             Description of file BORN
C             No default, BORN must be specified and cannot be blank.
C     IREC=positive integer... Index of the line in file BORN
C             corresponding to the first receiver.
C             Default: IREC=1
C     NREC=positive integer... Number of receivers read from file BORN.
C             Receivers IREC to IREC+NREC-1 are thus considered.
C             Default: NREC=1
C Names of the input files describing the incident wavefield:
C     The files have form MGRD (Multivalued GRiDded data).
C     NUM='string'... Name of the file with N1*N2*N3 array of integer
C             values, containing the number of arrivals at each
C             gridpoint of the grid.
C             No default, NUM must be specified and cannot be blank.
C     MTT='string'... Name of the file with the array of values,
C             containing for each gridpoint the travel times
C             of all determined arrivals.
C             No default, MTT must be specified and cannot be blank.
C     MP1='string', MP2='string', MP3='string'... Names of the files
C             with the array of values, containing the first, second and
C             third component of the slowness vector (i.e. with the
C             derivatives of travel time with respect to the coordinate
C             of the point).
C             No default, MP1, MP2, MP3 must be specified and cannot be
C             blank.
C     MTTXX='string'... Name of the file with the grid values of the
C             second derivative of travel time with respect to the
C             coordinate perpendicular to the 2D grid.
C             This file is used for computing the Born approximation of
C             a 3D wavefield in a 2D grid covering a 2D velocity model,
C             with the source and the receivers situated in the symmetry
C             plane of the 2D model.
C             MTTXX is not used and need not be specified for a 3D grid.
C             No default, MTTXX must be specified and cannot be blank
C             for a 2D grid.
C     AUR1='string', AUI1='string',
C     AUR2='string', AUI2='string',
C     AUR3='string', AUI3='string'... Names of the files with the array
C             of values containing the real or imaginary part of the
C             amplitude of the wave field.
C             No default, AUR1, AUI1, AUR2, AUI2, AUR3, AUI3 must be
C             specified and cannot be blank.
C Input files with source and receiver coordinates:
C     SRC='string'... String with the name of the input data file
C             with the name and coordinates of the source.
C             The point names cannot be longer than 6 characters.
C             Description of file SRC
C             No default, SRC must be specified and cannot be blank.
C     ISRC=positive integer... Index of the point in file SRC
C             corresponding to the source.
C             Default: ISRC=1
C     REC='string'... String with the name of the input data file
C             with the coordinates corresponding to the receiver names.
C             The point names cannot be longer than 6 characters.
C             Description of file REC
C             No default, REC must be specified and cannot be blank.
C Output file:
C     RF='string'... Name of the output file with the real and imaginary
C             parts of the response function corresponding to the Born
C             approximation.
C             Description of file RF.
C             Default: RF='rf.out'
C Data specifying dimensions of the input grid:
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     D1=real... Grid spacing along the X1 axis.
C             Default: D1=1.
C     D2=real... Grid spacing along the X2 axis.
C             Default: D2=1.
C     D3=real... Grid spacing along the X3 axis.
C             Default: D3=1.
C Data describing the frequency domain common with program 'ss.for':
C     These data are used to generate the optimum default values of
C             input parameters DF, OF and NF compatible with data for
C             '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.
C             Default: FMIN=0 (mostly sufficient)
C     FMAX=real... The highest frequency.
C             Default: FMAX=0.5/DT (Nyquist frequency, mostly
C               sufficient)
C Optional data for alternative specification of the frequency domain:
C     DF=real... Frequency step.
C             Default: DF=1/(NFFT*DT) (recommended)
C     OF=real... Born approximation is calculated for NF frequencies
C             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 Data specifying the cosine window used in the Born approximation.
C             The length of the cosine tapper is L.  The cosine window
C             vanishes at the distance greater than half the grid
C             interval from each side of the grid.
C     CSWIN=real... The global specification of parameter named L in the
C             specification above for cosine window. If none of
C             parameters CSWIN1,CSWIN2,CSWIN3,CSWIN4,CSWIN5,CSWIN6 is
C             specified, L is set to this value.
C             Default: CSWIN=0
C     CSWIN1=real... Parameter named L in the specification above for
C             the cosine window in the negative direction of the 1-st
C             axis.
C             Default: CSWIN1=CSWIN
C     CSWIN2=real... Parameter named L in the specification above for
C             the cosine window in the positive direction of the 1-st
C             axis.
C             Default: CSWIN2=CSWIN
C     CSWIN3=real... Parameter named L in the specification above for
C             the cosine window in the negative direction of the 2-nd
C             axis.
C             Default: CSWIN3=CSWIN
C     CSWIN4=real... Parameter named L in the specification above for
C             the cosine window in the positive direction of the 2-nd
C             axis.
C             Default: CSWIN4=CSWIN
C     CSWIN5=real... Parameter named L in the specification above for
C             the cosine window in the negative direction of the 3-rd
C             axis.
C             Default: CSWIN5=CSWIN
C     CSWIN6=real... Parameter named L in the specification above for
C             the cosine window in the positive direction of the 3-rd
C             axis.
C             Default: CSWIN6=CSWIN
C
C                                                    
C Input file BORN containing names of the receivers and the names of the
C files containing the quantities necessary for the Born approximation:
C     The quantities are discretized on the grid and are related to the
C     waves arriving to the receiver from each gridpoint of the grid.
C     The files have form MGRD (Multivalued GRiDded data).
C (1) For each receiver data (1.1):
C (1.1) TXTREC,NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG,
C     AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31,
C     AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32,
C     AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33,/
C     TXTREC... Name of the receiver.
C             The name cannot be longer than 6 characters.
C     NUMG... Name of the file with N1*N2*N3 array of integer
C             values containing the number of arrivals.
C     MTTG... Name of the file containing travel times.
C     MPG1,MPG2,MPG3... 3 names of the files containg the first,
C             second and third component of the slowness vector
C             (i.e. with the derivatives of travel time with respect to
C             the coordinates of the point).
C     MTTXXG... Name of the file with the grid values of the second
C             derivative of travel time with respect to the coordinate
C             perpendicular to the 2D grid.
C             This file is used for computing the Born approximation of
C             a 3D wavefield in a 2D grid covering a 2D velocity model,
C             with the source and the receivers situated in the symmetry
C             plane of the 2D model.
C             MTTXXG is not used and may be blank for a 3D grid.
C             MTTXXG cannot be blank for a 2D grid.
C     AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31,
C     AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32,
C     AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33... 18 names of the
C             files containing the real or imaginary part of the
C             tensorial amplitude of the Green function.
C     The filenames have no defaults. They must be specified.
C Input file BORN may contain the above data only, without any
C additional line.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FNAME,FILREC,FILSRC
      CHARACTER*80 NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG
      CHARACTER*80 AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31,
     &             AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32,
     &             AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33
      INTEGER LU1,LU2,LU3,LU4
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4)
      REAL PI,P2I,GIANT
      PARAMETER (PI=3.14159265,P2I=6.28318530,GIANT=2000000000.)
C
C     Input data:
      INTEGER N1,N2,N3
      INTEGER IREC,NREC
      INTEGER NPTS,NF
      REAL D1,D2,D3
      REAL TD,FMIN,FMAX,FSTEP
C
C     Data needed for output file
      INTEGER ISRC
      CHARACTER*6 TXTREC
      REAL XS1,XS2,XS3
      REAL XR1,XR2,XR3
      REAL TTMIN,TTMAX
C
C     Memory allocation:
      INTEGER INUMU,ILA,IMU,IRHO
      INTEGER IMTTU,IPU1,IPU2,IPU3,ITTXX,IUR1,IUI1,IUR2,IUI2,IUR3,IUI3
      INTEGER IBORN
      INTEGER INUMG,IMTTG,IPG1,IPG2,IPG3,ITTXXG
      INTEGER IGR11,IGI11,IGR21,IGI21,IGR31,IGI31,IGR12,IGI12,IGR22
      INTEGER IGI22,IGR32,IGI32,IGR13,IGI13,IGR23,IGI23,IGR33,IGI33
C
C     Auxiliary names for
C       perturbations of elastic moduli and density:
      REAL RLA,RMU,RHO
C       travel times and slowness vectors of the incident wavefield:
      REAL TTU,PU1,PU2,PU3
C       second derivative of travel times of the incident wavefield in
C       the direction perpendicular to the symmetry plane of the
C       2D model.
      REAL TTXX
C       real and imaginary parts of amplitudes of incident wavefield:
      REAL UR1,UR2,UR3
      REAL UI1,UI2,UI3
C       travel times and slowness vectors of the Green function:
      REAL TTG,PG1,PG2,PG3
C       second derivative of travel times of the Green function in the
C       direction perpendicular to the symmetry plane of the 2D model.
      REAL TTXXG
C       real and imaginary parts of amplitudes of the Green function:
      REAL GR11,GR21,GR31,GR12,GR22,GR32,GR13,GR23,GR33
      REAL GI11,GI21,GI31,GI12,GI22,GI32,GI13,GI23,GI33
C       travel time of the incident wavefield plus travel time of
C       the Green function
      REAL TT
C
C     Other variables:
      LOGICAL LINPUT,LRLA,LRMU,LRHO,LVPPER,LVSPER,LLRLA,LLRMU,LLRHO,L2D
      CHARACTER*80 MSSG,NAMPTS,MRECN
      CHARACTER*1  TEXT
      INTEGER I,II,J,JJ,K,IR,JF,ICOUNT
      INTEGER N12,N23,N123,N,NN,N6F,NWIN1,NWIN2,NWIN3,NWIN4,NWIN5,NWIN6
      INTEGER J1,J2,J3,J4,J5,J6,J6F
      INTEGER I01,I02,I03,I04,I05,I06,I07,I08,I09,I10,I11,I12
      REAL D123,ARG,ARGINC,RL,RL1,RL2,RL3,RL4,RL5,RL6,WINDOW
      REAL AUAG1R,AUAG1I,AUAG2R,AUAG2I,AUAG3R,AUAG3I
      REAL PGAG1R,PGAG1I,PGAG2R,PGAG2I,PGAG3R,PGAG3I
      REAL PUAUR,PUAUI
      REAL PUAG1R,PUAG1I,PUAG2R,PUAG2I,PUAG3R,PUAG3I
      REAL PGAUR,PGAUI,PUPG
      REAL DER,DEI,B1R,B1I,B2R,B2I,B3R,B3I,BR,BI
      REAL F,DOM2
      REAL AMAX
C
C     Functions
      INTEGER LENGTH
C
C-----------------------------------------------------------------------
C
C     Reading input data
C
C     Reading input SEP parameter file:
      WRITE(*,'(A)') '+GRDBORN: Enter input filename: '
      FSEP=' '
      READ(*,*) FSEP
      IF(FSEP.EQ.' ') THEN
C       GRDBORN-01
        CALL ERROR('GRDBORN-01: SEP file not given.')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      END IF
      CALL RSEP1(LU1,FSEP)
      WRITE(*,'(A)') '+GRDBORN: Working ...           '
C
C     Reading input data specifying dimensions of the input grid:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      N12=N1*N2
      N23=N2*N3
      N123=N1*N2*N3
C
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      L2D=.FALSE.
      IF(N2.EQ.1.AND.N3.EQ.1) THEN
C       D123=ABS(D1)
C       WRITE(*,'(A)') '+GRDBORN: 1D Born approximation on a 1D grid.'
C       GRDBORN-14
        CALL ERROR('GRDBORN-14: 1D Born approximation is not coded.')
C       The Born approximation of 1D, 2D or 3D wavefield on a 1D grid
C       has not been coded yet.
      ELSE IF(N1.EQ.1.AND.N3.EQ.1) THEN
C       D123=ABS(D2)
C       WRITE(*,'(A)') '+GRDBORN: 1D Born approximation on a 1D grid.'
C       GRDBORN-15
        CALL ERROR('GRDBORN-15: 1D Born approximation is not coded.')
C       The Born approximation of 1D, 2D or 3D wavefield on a 1D grid
C       has not been coded yet.
      ELSE IF(N1.EQ.1.AND.N2.EQ.1) THEN
C       D123=ABS(D2)
C       WRITE(*,'(A)') '+GRDBORN: 1D Born approximation on a 1D grid.'
C       GRDBORN-16
        CALL ERROR('GRDBORN-16: 1D Born approximation is not coded.')
C       The Born approximation of 1D, 2D or 3D wavefield on a 1D grid
C       has not been coded yet.
      ELSE IF(N1.EQ.1) THEN
        D123=ABS(D2*D3)
        WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 2D grid.'
        L2D=.TRUE.
C       Note that the Born approximation of 2D wavefield on a 2D grid
C       has not been coded yet.
      ELSE IF(N2.EQ.1) THEN
        D123=ABS(D1*D3)
        WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 2D grid.'
        L2D=.TRUE.
C       Note that the Born approximation of 2D wavefield on a 2D grid
C       has not been coded yet.
      ELSE IF(N3.EQ.1) THEN
        D123=ABS(D1*D2)
        WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 2D grid.'
        L2D=.TRUE.
C       Note that the Born approximation of 2D wavefield on a 2D grid
C       has not been coded yet.
      ELSE
        D123=ABS(D1*D2*D3)
        WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 3D grid.'
      END IF
C
C     Checking the size of RAM:
C     Is array RAM(MRAM) allocated in include file 'ram.inc' big enough
C     to store content of file specified by NUM and to read medium
C     perturbations?
C     You may wish to increase the dimension MRAM in file
C     ram.inc.
      CALL CHSP(MRAM,4*N123+1,'perturbations of density')
C
C     Reading data specifying dimensions of the output grid:
      CALL RSEP3I('IREC',IREC,1)
      CALL RSEP3I('NREC',NREC,1)
C
C     Reading data specifying the cosine windows.
      CALL RSEP3R('CSWIN',RL,0.)
      CALL RSEP3R('CSWIN1',RL1,RL)
      CALL RSEP3R('CSWIN2',RL2,RL)
      CALL RSEP3R('CSWIN3',RL3,RL)
      CALL RSEP3R('CSWIN4',RL4,RL)
      CALL RSEP3R('CSWIN5',RL5,RL)
      CALL RSEP3R('CSWIN6',RL6,RL)
C
C     Reading data describing the frequency domain:
      CALL RSEP3I('NFFT ',NPTS  ,1)
      CALL RSEP3R('DT   ',TD    ,1.)
      CALL RSEP3R('FMIN ',FMIN  ,0.)
      CALL RSEP3R('FMAX ',FMAX  ,0.5/TD)
      CALL RSEP3R('DF   ',FSTEP ,1./(FLOAT(NPTS)*TD))
      CALL RSEP3R('OF   ',FMIN  ,FSTEP*NINT(FMIN/FSTEP))
      CALL RSEP3I('NF   ',NF    ,NINT((FMAX-FMIN)/FSTEP)+1)
C     Auxiliary variable
      N6F=6*NF
C
C     Positions of the first elements of appropriate quantities
C     in array  RAM minus 1.
C     Quantities:number of arrivals for incident wave field.
      INUMU=1
      IRAM(INUMU)=0
C
C     Reading the numbers of arrivals for incident wave field
      CALL CHLI('NUM',' ',LU1,N123,.TRUE.,IRAM(INUMU+1),LINPUT)
C
C     Sequential number of the lastly read arrival at the I-th point:
      DO 10 I=INUMU+1,INUMU+N123
        IRAM(I)=IRAM(I)+IRAM(I-1)
   10 CONTINUE
C
      N=IRAM(INUMU+N123)
C     Reading input data specifying model and computation of positions
C     of the first elements of appropriate quantities in array RAM
C     minus 1.
C     Quantities:perturbations of rho*(v_s)^2,
C                perturbations of rho*(v_p)^2,density perturbations,
C                travel times,slowness vectors and amplitudes
C                of the incident wavefield,Born approximation,
C                number of arrivals and travel times of the Green
C                function
      IMU=INUMU+N123
C
      CALL CHLR('VSPER',' ',LU1,N123,.FALSE.,RAM(IMU+1),LVSPER)
      LRMU=LVSPER
      IF(LRMU) THEN
        ILA=IMU+N123
      ELSE
        ILA=IMU
      END IF
C
      CALL CHLR('VPPER',' ',LU1,N123,.FALSE.,RAM(ILA+1),LVPPER)
      IF(LVPPER.OR.LVSPER) THEN
        LRLA=.TRUE.
      ELSE
        LRLA=.FALSE.
      END IF
      IF(LRLA) THEN
        IRHO=ILA+N123
      ELSE
        IRHO=ILA
      END IF
C
      CALL CHLR('RHOPER',' ',LU1,N123,.FALSE.,RAM(IRHO+1),LRHO)
      IF(LRHO) THEN
        IMTTU=IRHO+N123
      ELSE
        IMTTU=IRHO
      END IF
C
      IPU1=IMTTU+N
      IPU2=IPU1+N
      IPU3=IPU2+N
C
      IF(L2D) THEN
        ITTXX=IPU3+N
      ELSE
        ITTXX=IPU3
      END IF
C
      IUR1=ITTXX+N
      IUI1=IUR1+N
      IUR2=IUI1+N
      IUI2=IUR2+N
      IUR3=IUI2+N
      IUI3=IUR3+N
C
      IBORN=IUI3+N
C
      INUMG=IBORN+N6F+1
      IRAM(INUMG)=0
C
      IMTTG=INUMG+N123
C
C     Checking the size of RAM:
C     Is array RAM(MRAM) allocated in include file 'ram.inc' big enough
C     to contain both input and output grids (without the Green
C     function)?
C     You may wish to increase the dimension MRAM in file
C     ram.inc.
      CALL CHSP(MRAM,INUMG+N123,'NUMG')
C
C     Reading input data specifying incident wavefield.
      CALL CHLR('MTT',' ',LU1,N,.TRUE.,RAM(IMTTU+1),LINPUT)
C
      CALL CHLR('MP1',' ',LU1,N,.TRUE.,RAM(IPU1+1),LINPUT)
      CALL CHLR('MP2',' ',LU1,N,.TRUE.,RAM(IPU2+1),LINPUT)
      CALL CHLR('MP3',' ',LU1,N,.TRUE.,RAM(IPU3+1),LINPUT)
C
      IF(L2D) THEN
        CALL CHLR('MTTXX',' ',LU1,N,.TRUE.,RAM(ITTXX+1),LINPUT)
        DO 17 I=1,N
          IF(RAM(ITTXX+I).LT.0.) THEN
C           GRDBORN-12
            CALL ERROR('GRDBORN-12: Negative second derivative of '
     &              // 'travel time in the direction perpendicular '
     &              // 'to the symmetry plane.')
          END IF
   17   CONTINUE
      END IF
C
      CALL CHLR('AUR1',' ',LU1,N,.TRUE.,RAM(IUR1+1),LINPUT)
      CALL CHLR('AUI1',' ',LU1,N,.TRUE.,RAM(IUI1+1),LINPUT)
      CALL CHLR('AUR2',' ',LU1,N,.TRUE.,RAM(IUR2+1),LINPUT)
      CALL CHLR('AUI2',' ',LU1,N,.TRUE.,RAM(IUI2+1),LINPUT)
      CALL CHLR('AUR3',' ',LU1,N,.TRUE.,RAM(IUR3+1),LINPUT)
      CALL CHLR('AUI3',' ',LU1,N,.TRUE.,RAM(IUI3+1),LINPUT)
C
C.......................................................................
C
C     Adjustment. Computation of perturbations of elastic moduli lambda
C     from perturbations of rho*(v_p)^2 and rho*(v_s)^2 :
C
      IF(LVPPER.AND.LVSPER) THEN
        DO 11 I=1,N123
          RAM(ILA+I)=RAM(ILA+I)-2.*RAM(IMU+I)
   11   CONTINUE
      ELSE IF(LVSPER) THEN
        DO 16 I=1,N123
          RAM(ILA+I)=-2.*RAM(IMU+I)
   16   CONTINUE
      END IF
C
C.......................................................................
C
C     Application of the cosine window on the model boundaries.
C     The affected gridpoints form a set of gridpoint planes
C     (two sets of gridpoint planes per an axis).
C     Number of gridpoint planes in each set:
      IF(N1.NE.1) THEN
        NWIN1=INT((RL1+D1/2)/ABS(D1))
        NWIN2=INT((RL2+D1/2)/ABS(D1))
      ELSE
        NWIN1=0
        NWIN2=0
      END IF
      IF(N2.NE.1) THEN
        NWIN3=INT((RL3+D2/2)/ABS(D2))
        NWIN4=INT((RL4+D2/2)/ABS(D2))
      ELSE
        NWIN3=0
        NWIN4=0
      END IF
      IF(N3.NE.1) THEN
        NWIN5=INT((RL5+D3/2)/ABS(D3))
        NWIN6=INT((RL6+D3/2)/ABS(D3))
      ELSE
        NWIN5=0
        NWIN6=0
      END IF
C
C     1-st axis
C     Variables II,J3 and JJ,J4 are pointers belonging to the first or
C     the second set of gridpoint planes.
C       Negative direction
C     Setting initial values of cosine argument, argument increment,
C     and pointer II.
      IF (NWIN1.GE.1) THEN
        ARG=PI*(D1/2.)/RL1
        ARGINC=2.*ARG
        II=1
      END IF
C     Loop over gridpoint values of the cosine window.
      DO 61 I=1,NWIN1
        WINDOW=(1.-COS(ARG))/2.
C       Pointer J3 is set to the value of pointer II. In the next loop
C       we will use and alter J3.
        J3=II
C       Loop over points with the same value of the cosine window
        DO 60 J=1,N23
          IF(LRHO) THEN
            J1=IRHO+J3
            RAM(J1)=RAM(J1)*WINDOW
          END IF
          IF(LRLA) THEN
            J1=ILA+J3
            RAM(J1)=RAM(J1)*WINDOW
          END IF
          IF(LRMU) THEN
            J1=IMU+J3
            RAM(J1)=RAM(J1)*WINDOW
          END IF
C         Pointer modification.
          J3=J3+N1
   60   CONTINUE
C       Argument increment
        ARG=ARG+ARGINC
C       Pointer modification
        II=II+1
   61 CONTINUE
C       Positive direction
C     Setting initial values of cosine argument, argument increment
C     and pointer JJ.
      IF (NWIN2.GE.1) THEN
        ARG=PI*(D1/2.)/RL2
        ARGINC=2.*ARG
        JJ=N1
      END IF
C     Loop over gridpoint values of the cosine window.
      DO 63 I=1,NWIN2
        WINDOW=(1.-COS(ARG))/2.
C       Pointer J4 is set to the value of pointer JJ. In the next loop
C       we will use and alter J4.
        J4=JJ
C       Loop over points with the same value of the cosine window
        DO 62 J=1,N23
          IF(LRHO) THEN
            J2=IRHO+J4
            RAM(J2)=RAM(J2)*WINDOW
          END IF
          IF(LRLA) THEN
            J2=ILA+J4
            RAM(J2)=RAM(J2)*WINDOW
          END IF
          IF(LRMU) THEN
            J2=IMU+J4
            RAM(J2)=RAM(J2)*WINDOW
          END IF
C         Pointer modification.
          J4=J4+N1
   62   CONTINUE
C       Argument increment
        ARG=ARG+ARGINC
C       Pointer modification
        JJ=JJ-1
   63 CONTINUE
C     The structure of the code in the following two cases is the same
C     (3-rd axis) or very similar (2-nd axis), therefore commentary is
C     briefer
C
C     2-nd axis
C     Variables II,J5,J3 and JJ,J6,J4 are pointers belonging to
C     the first or the second set of gridpoint planes.
C       Negative direction
      IF (NWIN3.GE.1) THEN
        ARG=PI*(D2/2.)/RL3
        ARGINC=2.*ARG
        II=1
      END IF
C     Loop over gridpoint values of the cosine window.
      DO 66 I=1,NWIN3
        WINDOW=(1.-COS(ARG))/2.
        J5=II
C       Two loops over points with the same value of the cosine window.
        DO 65 K=1,N3
          J3=J5
          DO 64 J=1,N1
            IF(LRHO) THEN
              J1=IRHO+J3
              RAM(J1)=RAM(J1)*WINDOW
            END IF
            IF(LRLA) THEN
              J1=ILA+J3
              RAM(J1)=RAM(J1)*WINDOW
            END IF
            IF(LRMU) THEN
              J1=IMU+J3
              RAM(J1)=RAM(J1)*WINDOW
            END IF
            J3=J3+1
   64     CONTINUE
          J5=J5+N12
   65   CONTINUE
        ARG=ARG+ARGINC
        II=II+N1
   66 CONTINUE
C       Positive direction
      IF (NWIN4.GE.1) THEN
        ARG=PI*(D2/2.)/RL4
        ARGINC=2.*ARG
        JJ=N12
      END IF
C     Loop over gridpoint values of the cosine window.
      DO 69 I=1,NWIN4
        WINDOW=(1.-COS(ARG))/2.
        J6=JJ
C       Two loops over points with the same value of the cosine window.
        DO 68 K=1,N3
          J4=J6
          DO 67 J=1,N1
            IF(LRHO) THEN
              J2=IRHO+J4
              RAM(J2)=RAM(J2)*WINDOW
            END IF
            IF(LRLA) THEN
              J2=ILA+J4
              RAM(J2)=RAM(J2)*WINDOW
            END IF
            IF(LRMU) THEN
              J2=IMU+J4
              RAM(J2)=RAM(J2)*WINDOW
            END IF
            J4=J4-1
   67     CONTINUE
          J6=J6+N12
   68   CONTINUE
        ARG=ARG+ARGINC
        JJ=JJ-N1
   69 CONTINUE
C
C     3-rd axis
C     Variables II,J3 and JJ,J4 are pointers belonging to the first or
C     the second set of gridpoint planes.
C       Negative direction
      IF (NWIN5.GE.1) THEN
        ARG=PI*(D3/2.)/RL5
        ARGINC=2.*ARG
        II=1
      END IF
C     Loop over gridpoint values of the cosine window.
      DO 71 I=1,NWIN5
        WINDOW=(1.-COS(ARG))/2.
        J3=II
C       Loop over points with the same value of the cosine window
        DO 70 J=1,N12
          IF(LRHO) THEN
            J1=IRHO+J3
            RAM(J1)=RAM(J1)*WINDOW
          END IF
          IF(LRLA) THEN
            J1=ILA+J3
            RAM(J1)=RAM(J1)*WINDOW
          END IF
          IF(LRMU) THEN
            J1=IMU+J3
            RAM(J1)=RAM(J1)*WINDOW
          END IF
          J3=J3+1
   70   CONTINUE
        ARG=ARG+ARGINC
        II=II+N12
   71 CONTINUE
C       Positive direction
      IF (NWIN6.GE.1) THEN
        ARG=PI*(D3/2.)/RL6
        ARGINC=2.*ARG
        JJ=N123
      END IF
C     Loop over gridpoint values of the cosine window.
      DO 73 I=1,NWIN6
        WINDOW=(1.-COS(ARG))/2.
        J4=JJ
C       Loop over points with the same value of the cosine window
        DO 72 J=1,N12
          IF(LRHO) THEN
            J2=IRHO+J4
            RAM(J2)=RAM(J2)*WINDOW
          END IF
          IF(LRLA) THEN
            J2=ILA+J4
            RAM(J2)=RAM(J2)*WINDOW
          END IF
          IF(LRMU) THEN
            J2=IMU+J4
            RAM(J2)=RAM(J2)*WINDOW
          END IF
          J4=J4-1
   72   CONTINUE
        ARG=ARG+ARGINC
        JJ=JJ-N12
   73 CONTINUE
C
C.......................................................................
C
C     Opening input file with the Green function:
C
      CALL RSEP3T('BORN',FNAME ,' ')
C
C     Checking
      IF(FNAME.EQ.' ') THEN
C       GRDBORN-02
        CALL ERROR('GRDBORN-02: Input file BORN not given.')
      END IF
C
      OPEN(LU1,FILE=FNAME,STATUS='OLD')
      DO 12 I=1,IREC-1
        READ(LU1,*) TXTREC,
     &              NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG,
     &              AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31,
     &              AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32,
     &              AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33
   12 CONTINUE
C
C.......................................................................
C
C     Opening the output file and writing parts (1),(2),(3)
C
C     Reading the name of the file with coordinates of the source
      CALL RSEP3T('SRC',FILSRC,' ')
      IF(FILSRC.EQ.' ') THEN
C       GRDBORN-03
        CALL ERROR('GRDBORN-03: File with coordinates of source was '
     &          // 'not specified.')
      END IF
C     Reading the position of desired source.
      CALL RSEP3I('ISRC',ISRC,1)
C     Reading the coordinates of the source
      OPEN(LU4,FILE=FILSRC,STATUS='OLD')
      READ(LU4,*) (TEXT,I=1,20)
      ICOUNT=0
   13 CONTINUE
        ICOUNT=ICOUNT+1
        NAMPTS='$$$$$$$'
        XS1=0.
        XS2=0.
        XS3=0.
        READ(LU4,*,END=14) NAMPTS,XS1,XS2,XS3
        IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
          GO TO 14
        END IF
        IF(NAMPTS(7:).NE.' ') THEN
C         GRDBORN-04
          CALL ERROR('GRDBORN-04: Source name exceeds 6 characters.')
C         Names of points in file given by input parameter PTS cannot
C         be longer than 6 characters.  This limitation is imposed by
C         the GSE standard.
        END IF
        IF(ICOUNT.EQ.ISRC) THEN
          GO TO 15
        END IF
      GO TO 13
   14 CONTINUE
C     GRDBORN-05
      CALL ERROR('GRDBORN-05: Source name and coordinates were not'
     &        // 'specified.')
   15 CONTINUE
      CLOSE(LU4)
C     Reading the name of the output file and opening that file
      CALL RSEP3T('RF',FNAME,'rf.out')
      OPEN(LU3,FILE=FNAME)
C     Writing.
      WRITE(LU3,'(A)') '/'
      WRITE(LU3,'(3(F9.3,1X),A)') XS1,XS2,XS3,'/'
      WRITE(LU3,'(2(E11.5,1X),I8,1X,A)') FMIN,FSTEP,NF,'/'
C
C.......................................................................
C
C     Computation
C
C     Reading the name of the file with coordinates of receivers
      CALL RSEP3T('REC',FILREC,' ')
      IF(FILREC.EQ.' ') THEN
C       GRDBORN-06
        CALL ERROR('GRDBORN-06: File with coordinates of receivers was '
     &          // 'not specified.')
      END IF
C     Opening the file with coordinates of receivers
      OPEN(LU4,FILE=FILREC,STATUS='OLD')
      READ(LU4,*) (TEXT,I=1,20)
C
C     Loop over receivers.
      DO 50 IR=1,NREC
        TTMIN=GIANT
        TTMAX=-GIANT
C       First of all,elements of array RAM which will contain Born
C       approximation have to be nullified.
        DO 20 I=IBORN+1,IBORN+N6F
          RAM(I)=0.
   20   CONTINUE
C
        TXTREC=' '
        READ(LU1,*,END=21) TXTREC,
     &                     NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG,
     &                     AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31,
     &                     AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32,
     &                     AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33
C
C       Checking
   21   CONTINUE
        IF(TXTREC.EQ.' ') THEN
          WRITE(MRECN,*) IREC+IR-1
C         GRDBORN-07
          WRITE(MSSG,'(3A)') 'GRDBORN-07: Name of the receiver number ',
     &                       MRECN(1:LENGTH(MRECN)),
     &                       ' was not specified.'
          CALL ERROR(MSSG)
C         Note: NUMG.EQ.' ' is caught by subroutine CHLI2 further on.
        END IF
C
C       Reading the coordinates of the receiver
   22   CONTINUE
          NAMPTS='$$$$$$$'
          XR1=0.
          XR2=0.
          XR3=0.
          READ(LU4,*,END=23) NAMPTS,XR1,XR2,XR3
          IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
            GO TO 23
          END IF
          IF(NAMPTS(7:).NE.' ') THEN
C           GRDBORN-08
            CALL ERROR('GRDBORN-08: Receiver name exceeds '
     &              // '6 characters.')
C           Names of points in file given by input parameter PTS cannot
C           be longer than 6 characters.  This limitation is imposed by
C           the GSE standard.
          END IF
          IF(TXTREC.EQ.NAMPTS) THEN
            GO TO 24
          END IF
        GO TO 22
   23   CONTINUE
C       GRDBORN-09
        CALL ERROR('GRDBORN-09: No point name matching the input '
     &          // 'receiver name.')
   24   CONTINUE
        REWIND(UNIT=LU4)
C
C       Reading the data specifying the number of arrivals for
C       the Green function.
        CALL CHLI2('NUMG',NUMG,LU2,N123,.TRUE.,IRAM(INUMG+1),LINPUT)
C
C       Sequential number of the lastly read arrival at the J-th
C       point:
        DO 25 J=INUMG+1,INUMG+N123
          IRAM(J)=IRAM(J)+IRAM(J-1)
   25   CONTINUE
C
        NN=IRAM(INUMG+N123)
C
C       Positions of the first elements of appropriate quantities
C       in array  RAM minus 1.
C       Quantities:slowness vectors and amplitudes of the Green function
C
        IPG1=IMTTG+NN
        IPG2=IPG1+NN
        IPG3=IPG2+NN
C
        IF(L2D) THEN
          ITTXXG=IPG3+NN
        ELSE
          ITTXXG=IPG3
        END IF
C
        IGR11=ITTXXG+NN
        IGI11=IGR11+NN
        IGR21=IGI11+NN
        IGI21=IGR21+NN
        IGR31=IGI21+NN
        IGI31=IGR31+NN
        IGR12=IGI31+NN
        IGI12=IGR12+NN
        IGR22=IGI12+NN
        IGI22=IGR22+NN
        IGR32=IGI22+NN
        IGI32=IGR32+NN
        IGR13=IGI32+NN
        IGI13=IGR13+NN
        IGR23=IGI13+NN
        IGI23=IGR23+NN
        IGR33=IGI23+NN
        IGI33=IGR33+NN
C
C       Checking
        CALL CHSP(MRAM,IGI33+NN,'Im(G_33)')
C
C       Reading the Green function
        CALL CHLR2('MTTG',MTTG,LU2,NN,.TRUE.,RAM(IMTTG+1),LINPUT)
C
        CALL CHLR2('MPG1',MPG1,LU2,NN,.TRUE.,RAM(IPG1+1),LINPUT)
        CALL CHLR2('MPG2',MPG2,LU2,NN,.TRUE.,RAM(IPG2+1),LINPUT)
        CALL CHLR2('MPG3',MPG3,LU2,NN,.TRUE.,RAM(IPG3+1),LINPUT)
C
        IF(L2D) THEN
          CALL CHLR2('MTTXXG',MTTXXG,LU2,NN,.TRUE.,RAM(ITTXXG+1),LINPUT)
          DO 26 I=1,NN
            IF(RAM(ITTXXG+I).LT.0.) THEN
C             GRDBORN-13
              CALL ERROR('GRDBORN-13: Negative second derivative of '
     &                // 'travel time in the direction perpendicular '
     &                // 'to the symmetry plane.')
            END IF
   26     CONTINUE
        END IF
C
        CALL CHLR2('AMPR11',AMPR11,LU2,NN,.TRUE.,RAM(IGR11+1),LINPUT)
        CALL CHLR2('AMPI11',AMPI11,LU2,NN,.TRUE.,RAM(IGI11+1),LINPUT)
C
        CALL CHLR2('AMPR21',AMPR21,LU2,NN,.TRUE.,RAM(IGR21+1),LINPUT)
        CALL CHLR2('AMPI21',AMPI21,LU2,NN,.TRUE.,RAM(IGI21+1),LINPUT)
C
        CALL CHLR2('AMPR31',AMPR31,LU2,NN,.TRUE.,RAM(IGR31+1),LINPUT)
        CALL CHLR2('AMPI31',AMPI31,LU2,NN,.TRUE.,RAM(IGI31+1),LINPUT)
C
        CALL CHLR2('AMPR12',AMPR12,LU2,NN,.TRUE.,RAM(IGR12+1),LINPUT)
        CALL CHLR2('AMPI12',AMPI12,LU2,NN,.TRUE.,RAM(IGI12+1),LINPUT)
C
        CALL CHLR2('AMPR22',AMPR22,LU2,NN,.TRUE.,RAM(IGR22+1),LINPUT)
        CALL CHLR2('AMPI22',AMPI22,LU2,NN,.TRUE.,RAM(IGI22+1),LINPUT)
C
        CALL CHLR2('AMPR32',AMPR32,LU2,NN,.TRUE.,RAM(IGR32+1),LINPUT)
        CALL CHLR2('AMPI32',AMPI32,LU2,NN,.TRUE.,RAM(IGI32+1),LINPUT)
C
        CALL CHLR2('AMPR13',AMPR13,LU2,NN,.TRUE.,RAM(IGR13+1),LINPUT)
        CALL CHLR2('AMPI13',AMPI13,LU2,NN,.TRUE.,RAM(IGI13+1),LINPUT)
C
        CALL CHLR2('AMPR23',AMPR23,LU2,NN,.TRUE.,RAM(IGR23+1),LINPUT)
        CALL CHLR2('AMPI23',AMPI23,LU2,NN,.TRUE.,RAM(IGI23+1),LINPUT)
C
        CALL CHLR2('AMPR33',AMPR33,LU2,NN,.TRUE.,RAM(IGR33+1),LINPUT)
        CALL CHLR2('AMPI33',AMPI33,LU2,NN,.TRUE.,RAM(IGI33+1),LINPUT)
C
        WRITE(*,'(A)') '+GRDBORN: Working.                            '
C       Loop over gridpoints.
        DO 30 I=1,N123
C         Auxiliary names and setting local logical variables, which
C         determine whether perturbations of medium parameters are zero
C         or not in the current gridpoint.
C
          IF(LRLA) THEN
            RLA=RAM(ILA+I)
            IF(RLA.NE.0.) THEN
              LLRLA=.TRUE.
            ELSE
              LLRLA=.FALSE.
            END IF
          ELSE
            LLRLA=.FALSE.
          END IF
C
          IF(LRMU) THEN
            RMU=RAM(IMU+I)
            IF(RMU.NE.0.) THEN
              LLRMU=.TRUE.
            ELSE
              LLRMU=.FALSE.
            END IF
          ELSE
            LLRMU=.FALSE.
          END IF
C
          IF(LRHO) THEN
            RHO=RAM(IRHO+I)
            IF(RHO.NE.0.) THEN
              LLRHO=.TRUE.
            ELSE
              LLRHO=.FALSE.
            END IF
          ELSE
            LLRHO=.FALSE.
          END IF
C
          IF(LLRHO.OR.LLRLA.OR.LLRMU) THEN
C         Density and/or elastic perturbations are not zero
C         in this gridpoint. We will compute its contribution
C         to Born approximation.
C
C           Loop over multivalued quantities of the Green function
            DO 29 JJ=IRAM(INUMG+I-1)+1,IRAM(INUMG+I)
C             Auxiliary names.
              TTG=RAM(IMTTG+JJ)
C
              PG1=RAM(IPG1+JJ)
              PG2=RAM(IPG2+JJ)
              PG3=RAM(IPG3+JJ)
C
              IF(L2D) THEN
                TTXXG=RAM(ITTXXG+JJ)
              END IF
C
              GR11=RAM(IGR11+JJ)
              GR21=RAM(IGR21+JJ)
              GR31=RAM(IGR31+JJ)
              GR12=RAM(IGR12+JJ)
              GR22=RAM(IGR22+JJ)
              GR32=RAM(IGR32+JJ)
              GR13=RAM(IGR13+JJ)
              GR23=RAM(IGR23+JJ)
              GR33=RAM(IGR33+JJ)
C
              GI11=RAM(IGI11+JJ)
              GI21=RAM(IGI21+JJ)
              GI31=RAM(IGI31+JJ)
              GI12=RAM(IGI12+JJ)
              GI22=RAM(IGI22+JJ)
              GI32=RAM(IGI32+JJ)
              GI13=RAM(IGI13+JJ)
              GI23=RAM(IGI23+JJ)
              GI33=RAM(IGI33+JJ)
C             Loop over multivalued quantities of the incident
C             wave field.
              DO 28 II=IRAM(INUMU+I-1)+1,IRAM(INUMU+I)
C               Auxiliary names.
                TTU=RAM(IMTTU+II)
C
                PU1=RAM(IPU1+II)
                PU2=RAM(IPU2+II)
                PU3=RAM(IPU3+II)
C
                IF(L2D) THEN
                  TTXX=RAM(ITTXX+II)
                END IF
C
                UR1=RAM(IUR1+II)
                UR2=RAM(IUR2+II)
                UR3=RAM(IUR3+II)
C
                UI1=RAM(IUI1+II)
                UI2=RAM(IUI2+II)
                UI3=RAM(IUI3+II)
C
C               Computing travel time from the source to the receiver.
                TT=TTU+TTG
C               And finding minimum and maximum travel time.
                TTMIN=AMIN1(TTMIN,TT)
                TTMAX=AMAX1(TTMAX,TT)
C
C               Computation of auxiliary variables and computation of
C               frequency independent parts of the Born approximation:
C
                B1R=0.
                B1I=0.
                B2R=0.
                B2I=0.
                B3R=0.
                B3I=0.
C
C               ........................................................
C
C               1-st term (and auxiliary variables for the 3-rd)
                IF(LLRHO.OR.LLRMU) THEN
C                 Auxiliary variables for the 1-st and the 3-rd term
                  AUAG1R=UR1*GR11+UR2*GR21+UR3*GR31
     &                  -UI1*GI11-UI2*GI21-UI3*GI31
                  AUAG1I=UI1*GR11+UI2*GR21+UI3*GR31
     &                  +UR1*GI11+UR2*GI21+UR3*GI31
C
                  AUAG2R=UR1*GR12+UR2*GR22+UR3*GR32
     &                  -UI1*GI12-UI2*GI22-UI3*GI32
                  AUAG2I=UI1*GR12+UI2*GR22+UI3*GR32
     &                  +UR1*GI12+UR2*GI22+UR3*GI32
C
                  AUAG3R=UR1*GR13+UR2*GR23+UR3*GR33
     &                  -UI1*GI13-UI2*GI23-UI3*GI33
                  AUAG3I=UI1*GR13+UI2*GR23+UI3*GR33
     &                  +UR1*GI13+UR2*GI23+UR3*GI33
C
                  IF(LLRHO) THEN
C                   Computation of the 1-st term
                    B1R=RHO*AUAG1R
                    B1I=RHO*AUAG1I
                    B2R=RHO*AUAG2R
                    B2I=RHO*AUAG2I
                    B3R=RHO*AUAG3R
                    B3I=RHO*AUAG3I
                  END IF
                END IF
C
C               ........................................................
C
C               2-nd term
                IF(LLRLA) THEN
C                 Auxiliary variables for the 2-nd term
                  PGAG1R=PG1*GR11+PG2*GR21+PG3*GR31
                  PGAG1I=PG1*GI11+PG2*GI21+PG3*GI31
C
                  PGAG2R=PG1*GR12+PG2*GR22+PG3*GR32
                  PGAG2I=PG1*GI12+PG2*GI22+PG3*GI32
C
                  PGAG3R=PG1*GR13+PG2*GR23+PG3*GR33
                  PGAG3I=PG1*GI13+PG2*GI23+PG3*GI33
C
                  PUAUR=PU1*UR1+PU2*UR2+PU3*UR3
                  PUAUI=PU1*UI1+PU2*UI2+PU3*UI3
C
C                 Computation of the 2-st term
                  B1R=B1R+RLA*(PGAG1R*PUAUR-PGAG1I*PUAUI)
                  B1I=B1I+RLA*(PGAG1I*PUAUR+PGAG1R*PUAUI)
                  B2R=B2R+RLA*(PGAG2R*PUAUR-PGAG2I*PUAUI)
                  B2I=B2I+RLA*(PGAG2I*PUAUR+PGAG2R*PUAUI)
                  B3R=B3R+RLA*(PGAG3R*PUAUR-PGAG3I*PUAUI)
                  B3I=B3I+RLA*(PGAG3I*PUAUR+PGAG3R*PUAUI)
                END IF
C
C               ........................................................
C
C               3-rd term
                IF(LLRMU) THEN
C                 Auxiliary variables for the 3-rd term
                  PUAG1R=PU1*GR11+PU2*GR21+PU3*GR31
                  PUAG1I=PU1*GI11+PU2*GI21+PU3*GI31
C
                  PUAG2R=PU1*GR12+PU2*GR22+PU3*GR32
                  PUAG2I=PU1*GI12+PU2*GI22+PU3*GI32
C
                  PUAG3R=PU1*GR13+PU2*GR23+PU3*GR33
                  PUAG3I=PU1*GI13+PU2*GI23+PU3*GI33
C
C
                  PGAUR=PG1*UR1+PG2*UR2+PG3*UR3
                  PGAUI=PG1*UI1+PG2*UI2+PG3*UI3
C
C
                  PUPG=PU1*PG1+PU2*PG2+PU3*PG3
C
C                 Computation of the 3-rd term
                  B1R=B1R+RMU*(PUPG*AUAG1R
     &                        +PUAG1R*PGAUR-PUAG1I*PGAUI)
                  B1I=B1I+RMU*(PUPG*AUAG1I
     &                        +PUAG1R*PGAUI+PUAG1I*PGAUR)
                  B2R=B2R+RMU*(PUPG*AUAG2R
     &                        +PUAG2R*PGAUR-PUAG2I*PGAUI)
                  B2I=B2I+RMU*(PUPG*AUAG2I
     &                        +PUAG2R*PGAUI+PUAG2I*PGAUR)
                  B3R=B3R+RMU*(PUPG*AUAG3R
     &                        +PUAG3R*PGAUR-PUAG3I*PGAUI)
                  B3I=B3I+RMU*(PUPG*AUAG3I
     &                        +PUAG3R*PGAUI+PUAG3I*PGAUR)
                END IF
C
C               ........................................................
C
C               Computation of term
C               (frequency independent part)*exp(i*2*pi*fmin(tauU+tauG))
C               Warning. DER here means real, DEI imaginary part of this
C               exponential. But DER and DEI will further have
C               a different meaning.
                DER=COS(P2I*FMIN*TT)
                DEI=SIN(P2I*FMIN*TT)
                IF(L2D) THEN
                  DER=DER/SQRT(TTXX+TTXXG)
                  DEI=DEI/SQRT(TTXX+TTXXG)
                END IF
C
                BR=B1R*DER-B1I*DEI
                BI=B1R*DEI+B1I*DER
                B1R=BR
                B1I=BI
C
                BR=B2R*DER-B2I*DEI
                BI=B2R*DEI+B2I*DER
                B2R=BR
                B2I=BI
C
                BR=B3R*DER-B3I*DEI
                BI=B3R*DEI+B3I*DER
                B3R=BR
                B3I=BI
C
                DER=COS(P2I*FSTEP*TT)
                DEI=SIN(P2I*FSTEP*TT)
C               The innermost loop over frequencies:
                DO 27 JF=IBORN+1,IBORN+N6F,6
C                 Computation of Born approximation (not multiplied
C                 by term frequency squared and term D123) using
C                 previously computed auxiliary variables:
C                 Variables J2-J6 are auxiliary for better performance.
                  J2=JF+1
                  J3=JF+2
                  J4=JF+3
                  J5=JF+4
                  J6=JF+5
C
                  RAM(JF)=RAM(JF)+B1R
                  RAM(J2)=RAM(J2)+B1I
                  RAM(J3)=RAM(J3)+B2R
                  RAM(J4)=RAM(J4)+B2I
                  RAM(J5)=RAM(J5)+B3R
                  RAM(J6)=RAM(J6)+B3I
C
                  BR=B1R*DER-B1I*DEI
                  BI=B1R*DEI+B1I*DER
                  B1R=BR
                  B1I=BI
C
                  BR=B2R*DER-B2I*DEI
                  BI=B2R*DEI+B2I*DER
                  B2R=BR
                  B2I=BI
C
                  BR=B3R*DER-B3I*DEI
                  BI=B3R*DEI+B3I*DER
                  B3R=BR
                  B3I=BI
   27           CONTINUE
   28         CONTINUE
   29       CONTINUE
          END IF
   30   CONTINUE
C
C       Adjustment of Born approximation.
C       Multiplication by term
C       1D, 3D: Frequency squared and D1*D2*D3 (in the code denoted by
C       D123}
C       2D: Frequency to the power of three halves, D1*D3 (in the code
C       denoted by D123), square root of pi and complex number 1+i.
        DO 31 JF=0,NF-1
          F=FMIN+FLOAT(JF)*FSTEP
C         Variables J1-J6,J6F,DOM2 are auxiliary for better performance.
          IF(L2D) THEN
            DOM2=((SQRT(P2I*F))**3)*D123*SQRT(PI)
          ELSE
            DOM2=D123*(P2I*F)**2.
          END IF
          J6F=IBORN+6*JF
          J1 = 1 + J6F
          J2 = 2 + J6F
          J3 = 3 + J6F
          J4 = 4 + J6F
          J5 = 5 + J6F
          J6 = 6 + J6F
C
          IF(L2D) THEN
            BR=DOM2*(RAM(J1)-RAM(J2))
            BI=DOM2*(RAM(J1)+RAM(J2))
            RAM(J1)=BR
            RAM(J2)=BI
            BR=DOM2*(RAM(J3)-RAM(J4))
            BI=DOM2*(RAM(J3)+RAM(J4))
            RAM(J3)=BR
            RAM(J4)=BI
            BR=DOM2*(RAM(J5)-RAM(J6))
            BI=DOM2*(RAM(J5)+RAM(J6))
            RAM(J5)=BR
            RAM(J6)=BI
          ELSE
            RAM(J1)=DOM2*RAM(J1)
            RAM(J2)=DOM2*RAM(J2)
            RAM(J3)=DOM2*RAM(J3)
            RAM(J4)=DOM2*RAM(J4)
            RAM(J5)=DOM2*RAM(J5)
            RAM(J6)=DOM2*RAM(J6)
          END IF
   31   CONTINUE
C
C.......................................................................
C
C       Saving results into the output file.
C
        AMAX=0.
        DO 40 I=IBORN+1,IBORN+N6F,6
          AMAX=AMAX1(ABS(RAM(I))  ,ABS(RAM(I+1)),
     *               ABS(RAM(I+2)),ABS(RAM(I+3)),
     *               ABS(RAM(I+4)),ABS(RAM(I+5)),
     *                                                     AMAX)
   40   CONTINUE
C
        IF(AMAX.LE.0.) THEN
          IF(TTMIN.GT.0.) THEN
            TTMAX=0.
          ELSE
            TTMAX=TTMIN-1.
          END IF
        END IF
C
        IF(XR1.LT.-99.999.OR.999.999.LT.XR1.OR.
     *     XR2.LT.-99.999.OR.999.999.LT.XR2.OR.
     *     XR3.LT.-99.999.OR.999.999.LT.XR3) THEN
          WRITE(LU3,'(3A,3(F7.0,1X),3(E11.5,1X),A)')
     *      '''',TXTREC,''' ',XR1,XR2,XR3,TTMIN,TTMAX,AMAX,'/'
        ELSE
          WRITE(LU3,'(3A,3(F7.3,1X),3(E11.5,1X),A)')
     *    '''',TXTREC,''' ',XR1,XR2,XR3,TTMIN,TTMAX,AMAX,'/'
        END IF
        IF(TTMIN.LE.TTMAX) THEN
          DO 41 I=IBORN+1,IBORN+N6F,12
            I01=IFIX(99999.1*RAM(I)/AMAX)
            I02=IFIX(99999.1*RAM(I+1)/AMAX)
            I03=IFIX(99999.1*RAM(I+2)/AMAX)
            I04=IFIX(99999.1*RAM(I+3)/AMAX)
            I05=IFIX(99999.1*RAM(I+4)/AMAX)
            I06=IFIX(99999.1*RAM(I+5)/AMAX)
            IF(I.LT.(IBORN+1+6*(NF-1))) THEN
              I07=IFIX(99999.1*RAM(I+6)/AMAX)
              I08=IFIX(99999.1*RAM(I+7)/AMAX)
              I09=IFIX(99999.1*RAM(I+8)/AMAX)
              I10=IFIX(99999.1*RAM(I+9)/AMAX)
              I11=IFIX(99999.1*RAM(I+10)/AMAX)
              I12=IFIX(99999.1*RAM(I+11)/AMAX)
              WRITE(LU3,'(12I6)') I01,I02,I03,I04,I05,I06,
     *                            I07,I08,I09,I10,I11,I12
            ELSE
              WRITE(LU3,'(12I6)') I01,I02,I03,I04,I05,I06
            END IF
   41     CONTINUE
        END IF
   50 CONTINUE
C
      WRITE(LU3,'(A)') '/'
C
      CLOSE(LU1)
      CLOSE(LU3)
      CLOSE(LU4)
      WRITE(*,'(A)') '+GRDBORN: Done.                            '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE CHLR(FNAME,TDEF,LU,N,LMSSG,ARRAY,INPUT)
      CHARACTER*80 MSSG
      CHARACTER*(*) FNAME,TDEF
      LOGICAL INPUT,LMSSG
      INTEGER LU,N,LENGTH
      REAL ARRAY(N)
C
C Subroutine CHeck and Load is designed to check whether there is any
C input filename in the input SEP parameter file. Filename is specified
C by parameter FNAME. If there is no input, variable INPUT will be false
C otherwise true.
C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is
C false, nothing will happen.
C If there is any input filename, real values contained in this file
C will be read from the file.  Otherwise no reading will be performed.
C
C Input:
C     FNAME...String containing the name of the parameter.  Except for
C             its case, it should match the parameter name in the input
C             SEP parameter file.
C     TDEF... Default value of the parameter.
C     LU...   Logical unit number to be used.
C     N...    Array dimension (number of elements to read).
C     LMSSG...Logical variable. Determines whether absence of input
C             file will create error or not.
C
C Output:
C     ARRAY.. Array having been read.
C     INPUT...Logical variable.  TRUE if input file has been specified,
C             FALSE otherwise.
C
C-----------------------------------------------------------------------
C
C  Other variables:
      CHARACTER*80 TOUT
C     TOUT... Value of the parameter specified by FNAME.
C
C-----------------------------------------------------------------------
C
      INPUT=.TRUE.
C
      CALL RSEP3T(FNAME,TOUT,TDEF)
C
      IF(TOUT.EQ.' ') THEN
        INPUT=.FALSE.
      END IF
C
      IF(INPUT) THEN
        CALL RARRAY(LU,TOUT,'FORMATTED',.TRUE.,0.,N,ARRAY)
      ELSE
        IF(LMSSG) THEN
C         GRDBORN-10
          WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ',
     &                        FNAME(1:LENGTH(FNAME)),
     &                       ' not specified.'
          CALL ERROR(MSSG)
        END IF
      END IF
      END
C
C=======================================================================
C
      SUBROUTINE CHLI(FNAME,TDEF,LU,N,LMSSG,IARRAY,INPUT)
      CHARACTER *80 MSSG
      CHARACTER*(*) FNAME,TDEF
      logical INPUT,LMSSG
      INTEGER LU,N,LENGTH
      INTEGER IARRAY(N)
C
C Subroutine CHeck and Load is designed to check whether there is any
C input filename in the input SEP parameter file. Filename is specified
C by parameter FNAME. If there is no input,variable INPUT will be false
C otherwise true.
C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is
C false, nothing will happen.
C If there is any input filename, integer values contained in this file
C will be read from the file. Otherwise no reading will be performed
C
C Input:
C     FNAME...String containing the name of the parameter. Except for
C             its case, it should match the parameter name in the input
C             SEP parameter file.
C     TDEF... Default value of the parameter.
C     LU...   Logical unit number to be used.
C     N...    Array dimension (number of elements to read).
C     LMSSG...Logical variable. Determines whether absence of input
C             file will create error or not.
C
C Output:
C     IARRAY..Array having been read.
C     INPUT...Logical variable.  TRUE if input file has been specified,
C             FALSE otherwise.
C
C-----------------------------------------------------------------------
C
C  Other variables:
      CHARACTER*80 TOUT
C     TOUT... Value of the parameter specified by FNAME.
C
C-----------------------------------------------------------------------
C
      INPUT=.TRUE.
C
      CALL RSEP3T(FNAME,TOUT,TDEF)
C
      IF(TOUT.EQ.' ') THEN
        INPUT=.FALSE.
      END IF
C
      IF(INPUT) THEN
        CALL RARRAI(LU,TOUT,'FORMATTED',.TRUE.,0,N,IARRAY)
      ELSE
        IF(LMSSG) THEN
          WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ',
     &                        FNAME(1:LENGTH(FNAME)),
     &                       ' not specified.'
          CALL ERROR(MSSG)
        END IF
      END IF
      END
C
C=======================================================================
C
      SUBROUTINE CHLR2(FNAME,TOUT,LU,N,LMSSG,ARRAY,INPUT)
      CHARACTER *80 MSSG
      CHARACTER*(*) FNAME,TOUT
      LOGICAL INPUT,LMSSG
      INTEGER LU,N,LENGTH
      REAL ARRAY(N)
C
C Subroutine CHeck and Load is designed to check whether there is any
C input filename. If not, variable INPUT will be false otherwise true.
C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is
C false, nothing will happen.
C If there is any input filename, real values contained in this file
C will be read from the file. Otherwise no reading will be performed
C
C Input:
C     FNAME...String containing the name of the parameter specifying
C             the file
C     TOUT... String containg the filename.
C     LU...   Logical unit number to be used.
C     N...    Array dimension (number of elements to read).
C     LMSSG...Logical variable. Determines whether absence of input
C             file will create error or not.
C
C Output:
C     ARRAY.. Array having been read.
C     INPUT...Logical variable.  TRUE if input file has been specified,
C             FALSE otherwise.
C
C-----------------------------------------------------------------------
C
      INPUT=.TRUE.
C
      IF(TOUT.EQ.' ') THEN
        INPUT=.FALSE.
      END IF
C
      IF(INPUT) THEN
        CALL RARRAY(LU,TOUT,'FORMATTED',.TRUE.,0.,N,ARRAY)
      ELSE
        IF(LMSSG) THEN
          WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ',
     &                        FNAME(1:LENGTH(FNAME)),
     &                       ' not specified.'
          CALL ERROR(MSSG)
        END IF
      END IF
      END
C
C=======================================================================
C
      SUBROUTINE CHLI2(FNAME,TOUT,LU,N,LMSSG,IARRAY,INPUT)
      CHARACTER *80 MSSG
      CHARACTER*(*) FNAME,TOUT
      LOGICAL INPUT,LMSSG
      INTEGER LU,N,LENGTH
      INTEGER IARRAY(N)
C
C Subroutine CHeck and Load is designed to check whether there is any
C input filename. If not, variable INPUT will be false otherwise true.
C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is
C false, nothing will happen.
C If there is any input filename, integer values contained in this file
C will be read from the file. Otherwise no reading will be performed
C
C Input:
C     FNAME...String containing the name of the parameter specifying
C             the file
C     TOUT... String containg the filename.
C     LU...   Logical unit number to be used.
C     N...    Array dimension (number of elements to read).
C     LMSSG...Logical variable. Determines whether absence of input file
C             will create error or not.
C
C Output:
C     IARRAY..Array having been read.
C     INPUT...Logical variable.  TRUE if input file has been specified,
C             FALSE otherwise.
C
C-----------------------------------------------------------------------
C
      INPUT=.TRUE.
C
      IF(TOUT.EQ.' ') THEN
        INPUT=.FALSE.
      END IF
C
      IF(INPUT) THEN
        CALL RARRAI(LU,TOUT,'FORMATTED',.TRUE.,0,N,IARRAY)
      ELSE
        IF(LMSSG) THEN
          WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ',
     &                        FNAME(1:LENGTH(FNAME)),
     &                       ' not specified.'
          CALL ERROR(MSSG)
        END IF
      END IF
      END
C
C=======================================================================
C
      SUBROUTINE CHSP(MEMORY,NEEDED,QNAME)
      INTEGER MEMORY,NEEDED,LENGTH
      CHARACTER *(*) QNAME
C
C Subroutine designed to CHeck SPace in memory which we would like
C to use to store data or as a computational space.
C
C Input:
C     MEMORY..Reserved memory space.
C     NEEDED..Needed memory space.
C     QNAME ..Name of the last quantity in the memory.
C
C-----------------------------------------------------------------------
C
C  Other variables:
      CHARACTER *80 MMEM,MNEED
      CHARACTER *140 MSSG
C
C-----------------------------------------------------------------------
C
      IF(MEMORY.LT.NEEDED) THEN
        WRITE(MMEM,*) MEMORY
        WRITE(MNEED,*) NEEDED
C       GRDBORN-11
        WRITE(MSSG,'(7A)') 'GRDBORN-11: Not enough space for ',
     &    QNAME(1:LENGTH(QNAME)),' in array RAM.',
     &    ' Available:',MMEM(1:LENGTH(MMEM)),
     &    ' Needed:',MNEED(1:LENGTH(MNEED))
        CALL ERROR(MSSG)
      END IF
      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