C 
C Program GBOPT to calculate quantities for optimization of the shape of
C 2-D Gaussian beams for prestack depth migrations
C
C Version: 8.00
C Date: 2022, April 24
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz
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 input files with the results of ray tracing from a receiver:
C     CRTOUT='string'... File with the names of the output files of
C             program CRT.
C             For general description of file CRTOUT refer to file
C            'writ.for'.
C             Description specific to this program:
C               Only the file
C              'CRT-R'
C               with the quantities stored along rays and the file
C              'CRT-I'
C               with the quantities at the initial points of rays
C               are read by GBOPT.
C               File 'CRT-R' must contain all rays traced by CRT, not
C               only two-point rays.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C             'CRT-I'='r01i.out'
C Names of input formatted files with quantities corresponding to the
C             common source, interpolated at the nodes of input grid:
C     NUM='string'...Name of the input file with N1*N2*N3 array
C             of integer values, containing the number of arrivals at
C             each gridpoint of the target grid.
C             Default: NUM='mtt-num.out'
C     MTT='string'...Name of the input file with the array of real
C             values, containing for each gridpoint the travel times
C             of all determined arrivals.
C             Default: MTT='mtt-tt.out'
C     AMP='string'...Name of the input file with the array of real
C             values, containing for each gridpoint the amplitudes.
C             Default: AMP='mtt-ap.out'
C Input parameters:
C     NTR=positive integer... Initial number of travel times
C             in the 3-D Hamiltonian hypersurface.
C             Default: NTR=1
C     NP=positive integer... Initial number of take-off angles
C             in the 3-D Hamiltonian hypersurface.
C             Default: NP=1
C Data specifying the dimensions of the input MTT grid:
C     O1=real, O2=real, O3=real ... Coordinates of the origin of the
C             grid.
C             Default: O1=0., O2=0., O3=0.
C     D1=real... Grid interval along the X1 axis.
C             Default: D1=1.
C     D2=real... Grid interval along the X2 axis.
C             Default: D2=1.
C     D3=real... Grid interval along the X3 axis.
C             Default: D3=1.
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 Names of output files:
C     FGBR22='string'... Name of the output file containing optimum
C             value R0=-C21/C22 of the real part of parameter M0.
C             Default: FGBR22='r0s.out'
C     FGBY22='string'... Name of the output file containing optimum
C             value Y0=SQRT(C11/C22)
C             of the imaginary part of parameter M0.
C             Default: FGBY22='y0s.out'
C     FSIGMA='string'... Name of the output file containing initial
C             value SIGMA=SQRT(Y0/C22) of a standard deviation for
C             smoothing.
C             Default: FSIGMA='sigmar.out'
C     FTAU='string'... Name of the output file containing sum of
C             travel times from source to point in model and from
C             intersection of central ray of Gaussian packet with
C             profile to same point in model.
C             Default: FTAU='tau.out'
C Form of the input/output files:
C     FORM='string'... Form of the input/output files: either
C             'FORMATTED' or 'UNFORMATTED' (case insensitive).
C             Default: FORM='FORMATTED'
C     FORMR='string'... Form of the input file. If both FORM and FORMR
C             are specified, FORMR is used for output files.
C             Default: FORMR=FORM
C Value of undefined quantities:
C     UNDEF=real... The value to be used for undefined real quantities.
C             Default: UNDEF=undefined value used in forms.for
C
C=======================================================================
C
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     External functions and subroutines:
      EXTERNAL AP00,AP28,OARRAY,RARRAY,RARRAI,ERROR,RSEP1,RSEP3T,RSEP3I,
     *         RSEP3R
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL O1,O2,O3,D1,D2,D3,NORM0,NORM1,NORM2,FAK1,FAK2,FAK3,DET
      REAL FX,FFX,COR3,COR6,COR7
      REAL H11,H12,H13,H21,H22,H23,H31,H32,H33,V1,V2,E22
      REAL C11,C12,C21,C22,Y0,R0,SIGMA,TAU,TT
      REAL C(10),DC(7),SDC(7),SSDC(3),FN1(7),FFN1(3),T(8)
      INTEGER LU1,LU2,LU3,LU4,LU5,LU6,IPUF,NTR,IP,NP,NPN
      INTEGER IFN1(7),IFN2(7),IFFN1(3),IFFN2(3),IRAM(MRAM)
      INTEGER K1,K2,K3,INDX,ISH,IUF,ISM,IGRD,ITAU,N1,N2,N3,N1N2N3
      INTEGER ISRF1,NFN1,NFFN1,NFN2,NFFN2
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4,LU5=5,LU6=6)
      CHARACTER*80 FILSEP,FILCRT,RAYALL,INIALL,CH,NUM,MTT,AMP
      CHARACTER*80 FGBY22,FGBR22,FSIGMA,FTAU
      EQUIVALENCE (IRAM,RAM)
C
C     Storage in array (I)RAM:
C     N1N2N3=N1*N2*N3
C     ISH=2*N1N2N3+1
C     IRAM(1:N1N2N3)... N1*N2*N3 array of integer values, containing
C             the number of arrivals at each gridpoint of the target
C             grid (NUM input file).
C     IRAM(N1N2N3+1:2*N1N2N3)... Elementary sums of all arrivals at
C             gridpoints of the target grid (counted from the origin
C             of the grid).
C     IRAM(ISH)... Sum of all arrivals at all gridpoints of
C             the target grid.
C     RAM(ISH+1:ISH+1+IRAM(ISH))... IRAM(ISH) array of values,
C             containing for each gridpoint the travel times of all
C             determined arrivals (MTT input file).
C     RAM(ISH+IRAM(ISH)+1:ISH+IRAM(ISH)+1+IRAM(ISH))... IRAM(ISH) grid
C             values of the norm of the vectorial amplitude of
C             the Green function (AMP input file).
C     RAM(ISH+2*IRAM(ISH)+1:ISH+2*IRAM(ISH)+1+N1N2N3)... N1*N2*N3 grid
C             of normalized values of multivalued travel times.
C             The value of travel time is multiplied by the norm
C             and the sum of normalized travel times is divided by
C             the number of travel times at the gridpoint.
C
C-----------------------------------------------------------------------
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GBOPT: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GBOPT-01
        CALL ERROR('GBOPT-01: No input file specified')
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
      WRITE(*,'(A)')'+GBOPT: Reading, calculating... '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading filenames of the input files with computed rays
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      RAYALL='r01.out'
      INIALL='r01i.out'
      IF(FILCRT.NE.' ') THEN
        OPEN (LU2,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU2,*) RAYALL,CH,INIALL
        CLOSE (LU2)
      ENDIF
C
C     Reading filenames of the input files with multivalued
C     travel times
      CALL RSEP3T('NUM',NUM,'mtt-num.out')
      CALL RSEP3T('MTT',MTT,'mtt-tt.out')
      CALL RSEP3T('AMP',AMP,'mtt-ap.out')
C
C     Input parameters
      CALL RSEP3I('NTR',NTR,1)
      CALL RSEP3I('NP',NP,1)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
C
C     Reading filenames of the output files
      CALL RSEP3T('FGBR22',FGBR22,'r0s.out')
      CALL RSEP3T('FGBY22',FGBY22,'y0s.out')
      CALL RSEP3T('FSIGMA',FSIGMA,'sigmar.out')
      CALL RSEP3T('FTAU',FTAU,'tau.out')
C
      N1N2N3=N1*N2*N3
      ISH=2*N1N2N3+1
      FX=0.
      FFX=0.
      IPUF=0
      NPN=1
      ISRCO=1
C
C     Read N1*N2*N3 array of integer values, containing the number
C     of arrivals at each grid point of the target grid
      CALL OARRAY(LU1,NUM,'R')
      CALL RARRAI(LU1,' ',.FALSE.,0,N1N2N3,IRAM)
      CLOSE(LU1)
C     Elementary sums of all arrivals at gridpoints of the target grid
C     (counted from the origin of the grid)
      IRAM(N1N2N3+1)=0
      DO 10 IGRD=1,N1N2N3
        IRAM(N1N2N3+IGRD+1)=IRAM(N1N2N3+IGRD)+IRAM(IGRD)
   10 CONTINUE
C
      IF(3*N1N2N3+2*IRAM(ISH)+1.GT.MRAM) THEN
C       GBOPT-02
        CALL ERROR('GBOPT-02: Too small array RAM(MRAM)')
C       Too small array RAM(MRAM).  If possible, increase dimension MRAM
C       in include file
C       ram.inc.
      END IF
C
C     Read array of values, containing for each grid point the travel
C     times of all determined arrivals
      CALL OARRAY(LU1,MTT,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,IRAM(ISH),
     *            RAM(ISH+1))
      CLOSE(LU1)
C
C     Read the grid values of the norm of the vectorial amplitude
C     of the Green function.
      CALL OARRAY(LU1,AMP,'R')
      CALL RARRAY(LU1,' ',.FALSE.,1.,IRAM(ISH),
     *            RAM(ISH+IRAM(ISH)+1))
      CLOSE(LU1)
C
C     Calculation of normalized travel times
      DO 30 IGRD=1,N1N2N3
        IF(IRAM(IGRD).GE.1) THEN
          NORM1=0.
          NORM2=0.
          DO 20 ITAU=1,IRAM(IGRD)
            NORM0=RAM(ISH+IRAM(ISH)+IRAM(N1N2N3+IGRD)+ITAU)**2
            NORM1=NORM1+RAM(ISH+IRAM(ISH)+IRAM(N1N2N3+IGRD)+ITAU)**2
            NORM2=NORM2+NORM0*RAM(ISH+IRAM(N1N2N3+IGRD)+ITAU)
   20     CONTINUE
          IF(NORM1.EQ.0.) THEN
            RAM(ISH+2*IRAM(ISH)+IGRD)=-1.
          ELSE
            RAM(ISH+2*IRAM(ISH)+IGRD)=NORM2/NORM1
          END IF
        ELSE
          RAM(ISH+2*IRAM(ISH)+IGRD)=-1.
        END IF
   30 CONTINUE
C
C     Opening the input files with computed rays:
      OPEN(LU1,FILE=RAYALL,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=INIALL,FORM='UNFORMATTED',STATUS='OLD')
C     Opening the output files:
      OPEN(LU3,FILE=FGBR22,FORM='FORMATTED')
      OPEN(LU4,FILE=FGBY22,FORM='FORMATTED')
      OPEN(LU5,FILE=FSIGMA,FORM='FORMATTED')
      OPEN(LU6,FILE=FTAU,FORM='FORMATTED')
C
C     Loop for the points of rays
   40 CONTINUE
C
C       Reading the results of the complete ray tracing
C
C       Subroutine AP00 reads from the CRT output files the quantities
C       describing a point on a ray
        CALL AP00(0,LU1,LU2)
C
        IF(IWAVE.LT.1) THEN
C         End of rays
          IF(IPUF.LT.NTR) THEN
C           Check number of travel times
            IPUF=NTR-IPUF
            DO 41 IP=1,IPUF
              WRITE(LU3,'(I7,A)') 1,'*'
              WRITE(LU4,'(I7,A)') 1,'*'
              WRITE(LU5,'(I7,A)') 1,'*'
              WRITE(LU6,'(I7,A)') 1,'*'
   41       CONTINUE
          ELSE
C           GBOPT-03
            CALL ERROR('GBOPT-03: Too low parameter NTR - End of rays')
          END IF
          NPN=NPN-1
          IF(NPN.LT.NP) THEN
C           Check number of slowness vectors
            NPN=NTR*(NP-NPN)
            DO 42 IP=1,NPN
              WRITE(LU3,'(I7,A)') 1,'*'
              WRITE(LU4,'(I7,A)') 1,'*'
              WRITE(LU5,'(I7,A)') 1,'*'
              WRITE(LU6,'(I7,A)') 1,'*'
   42       CONTINUE
          END IF
          GO TO 100
        END IF
C
        IF(ISRC.NE.ISRCO) THEN
C         New source
          IF(IPUF.LT.NTR) THEN
C           Check number of travel times
            IPUF=NTR-IPUF
            DO 43 IP=1,IPUF
              WRITE(LU3,'(I7,A)') 1,'*'
              WRITE(LU4,'(I7,A)') 1,'*'
              WRITE(LU5,'(I7,A)') 1,'*'
              WRITE(LU6,'(I7,A)') 1,'*'
   43       CONTINUE
          ELSE
C           GBOPT-03
            CALL ERROR('GBOPT-03: Too low parameter NTR - New source')
          END IF
          NPN=NPN-1
          IF(NPN.LT.NP) THEN
C           Check number of slowness vectors
            NPN=NTR*(NP-NPN)
            DO 44 IP=1,NPN
              WRITE(LU3,'(I7,A)') 1,'*'
              WRITE(LU4,'(I7,A)') 1,'*'
              WRITE(LU5,'(I7,A)') 1,'*'
              WRITE(LU6,'(I7,A)') 1,'*'
   44       CONTINUE
          END IF
          ISRCO=ISRC
          IPUF=0
        END IF
C
        IF(IPT.LE.1) THEN
C         New ray
          IF(IRAY.GT.NPN) THEN
C           Check number of slowness vectors
            NPN=NTR*(IRAY-NPN)
            DO 45 IP=1,NPN
              WRITE(LU3,'(I7,A)') 1,'*'
              WRITE(LU4,'(I7,A)') 1,'*'
              WRITE(LU5,'(I7,A)') 1,'*'
              WRITE(LU6,'(I7,A)') 1,'*'
   45       CONTINUE
          END IF
          NPN=IRAY+1
          IF(IRAY.GT.1) THEN
            IF(IPUF.LT.NTR) THEN
C           Check number of travel times
              IPUF=NTR-IPUF
              DO 48 IP=1,IPUF
                WRITE(LU3,'(I7,A)') 1,'*'
                WRITE(LU4,'(I7,A)') 1,'*'
                WRITE(LU5,'(I7,A)') 1,'*'
                WRITE(LU6,'(I7,A)') 1,'*'
   48         CONTINUE
            ELSE
C             GBOPT-04
              CALL ERROR('GBOPT-04: Too low parameter NTR')
            END IF
          END IF
C
C         Initial values for numerical quadrature at the initial point
          IPUF=0
          NFN1=7
          NFFN1=3
          NFN2=7
          NFFN2=3
          DO 50 ISM=1,7
           IFN1(ISM)=ISM
           FN1(ISM)=0.
           IFN2(ISM)=ISM
   50     CONTINUE
          DO 55 ISM=1,3
           IFFN1(ISM)=ISM
           FFN1(ISM)=0.
           IFFN2(ISM)=ISM
   55     CONTINUE
          FFN1(1)=0.25
          FFN1(3)=0.25
        END IF
C
C       Numerical quadrature of the differential equation (31)
        DC(1)=Y(12)*Y(20)+Y(13)*Y(21)
        DC(2)=Y(16)*Y(20)+Y(17)*Y(21)
        DC(3)=Y(20)*Y(20)+Y(21)*Y(21)
        DC(4)=Y(12)*Y(24)+Y(13)*Y(25)
        DC(5)=Y(16)*Y(24)+Y(17)*Y(25)
        DC(6)=Y(20)*Y(24)+Y(21)*Y(25)
        DC(7)=Y(24)*Y(24)+Y(25)*Y(25)
C
C       Subroutine AP28 performs the numerical quadrature of the set
C       of given functions along a ray.  It has to be called once at
C       each point along the ray in which the computed quantities
C       are stored
        CALL AP28(7,SDC,1,1,0.01,FX,ISRF1,NFN1,IFN1,FN1,NFN2,IFN2,DC)
C
C       Correction for trapezoidal quadrature of a quadratic function
        IF(IPT.LE.1) THEN
          COR3=SDC(3)/3.
          COR6=SDC(6)/3.
          COR7=SDC(7)/3.
        END IF
        SDC(3)=SDC(3)-COR3
        SDC(6)=SDC(6)-COR6
        SDC(7)=SDC(7)-COR7
C
        C(4)=SDC(1)
        C(5)=SDC(2)
        C(6)=SDC(3)
        C(7)=SDC(4)
        C(8)=SDC(5)
        C(9)=SDC(6)
        C(10)=SDC(7)
        DC(4)=Y(12)
        DC(5)=Y(13)
        DC(6)=Y(16)
        DC(7)=Y(17)
        DET=C(6)*C(10)-C(9)*C(9)
        IF(DET.EQ.0.) THEN
C             GBOPT-05
              CALL ERROR('GBOPT-05: Determinant is equal to zero')
        END IF
        DC(4)=(Y(20)*C(10)-Y(24)*C(9))*C(4)
        DC(4)=DC(4)+(Y(24)*C(6)-Y(20)*C(9))*C(7)
        DC(4)=DC(4)/DET
        DC(4)=Y(12)-DC(4)
C
        DC(5)=(Y(21)*C(10)-Y(25)*C(9))*C(4)
        DC(5)=DC(5)+(Y(25)*C(6)-Y(21)*C(9))*C(7)
        DC(5)=DC(5)/DET
        DC(5)=Y(13)-DC(5)
C
        DC(6)=(Y(20)*C(10)-Y(24)*C(9))*C(5)
        DC(6)=DC(6)+(Y(24)*C(6)-Y(20)*C(9))*C(8)
        DC(6)=DC(6)/DET
        DC(6)=Y(16)-DC(6)
C
        DC(7)=(Y(21)*C(10)-Y(25)*C(9))*C(5)
        DC(7)=DC(7)+(Y(25)*C(6)-Y(21)*C(9))*C(8)
        DC(7)=DC(7)/DET
        DC(7)=Y(17)-DC(7)
C
        DC(1)=DC(4)**2+DC(5)**2
        DC(2)=DC(4)*DC(6)+DC(5)*DC(7)
        DC(3)=DC(6)**2+DC(7)**2
C
C       Subroutine AP28 performs the numerical quadrature of the set
C       of given functions along a ray.  It has to be called once at
C       each point along the ray in which the computed quantities
C       are stored
        CALL AP28(3,SSDC,1,1,0.01,FFX,ISRF1,
     *                                  NFFN1,IFFN1,FFN1,NFFN2,IFFN2,DC)
C
        C(1)=SSDC(1)
        C(2)=SSDC(2)
        C(3)=SSDC(3)
C
C       Calculation of the mean arrival time at the reciever,
C       eight grid points are used
        IUF=1
        K1=INT((Y(3)-O1)/D1)
        K2=INT((Y(4)-O2)/D2)
        K3=INT((Y(5)-O3)/D3)
C       Check position of gridpoints
        IF(K1.LT.0.OR.(N1.GT.1.AND.K1.GE.(N1-1))) THEN
          IUF=0
        END IF
        IF(K2.LT.0.OR.(N2.GT.1.AND.K2.GE.(N2-1))) THEN
          IUF=0
        END IF
        IF(K3.LT.0.OR.(N3.GT.1.AND.K3.GE.(N3-1))) THEN
          IUF=0
        END IF
C
        IF(IUF.EQ.1) THEN
          INDX=K3*N1*N2+K2*N1+K1+1
          T(1)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF(N1.GT.1) THEN
            INDX=K3*N1*N2+K2*N1+K1+2
          ELSE
            INDX=K3*N1*N2+K2*N1+1
          END IF
          T(2)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF(N2.GT.1) THEN
            INDX=K3*N1*N2+(K2+1)*N1+K1+1
          ELSE
            INDX=K3*N1*N2+K1+1
          END IF
          T(3)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF((N1-1)*(N2-1).GT.0) THEN
            INDX=K3*N1*N2+(K2+1)*N1+K1+2
          ELSE
            IF(N1.GT.1)THEN
              INDX=K3*N1*N2+K1+2
            ELSE
              INDX=K3*N1*N2+(K2+1)*N1+1
            END IF
          END IF
          T(4)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF(N3.GT.1) THEN
            INDX=(K3+1)*N1*N2+K2*N1+K1+1
          ELSE
            INDX=K2*N1+K1+1
          END IF
          T(5)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF((N1-1)*(N3-1).GT.0) THEN
            INDX=(K3+1)*N1*N2+K2*N1+K1+2
          ELSE
            IF(N1.GT.1)THEN
              INDX=K2*N1+K1+2
            ELSE
              INDX=(K3+1)*N1*N2+K2*N1+1
            END IF
          END IF
          T(6)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF((N2-1)*(N3-1).GT.0) THEN
            INDX=(K3+1)*N1*N2+(K2+1)*N1+K1+1
          ELSE
            IF(N2.GT.1)THEN
              INDX=(K2+1)*N1+K1+1
            ELSE
              INDX=(K3+1)*N1*N2+K1+1
            END IF
          END IF
          T(7)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          IF((N1-1)*(N2-1)*(N3-1).GT.0) THEN
            INDX=(K3+1)*N1*N2+(K2+1)*N1+K1+2
          ELSE
            IF(N1.GT.1) THEN
              IF(N2.GT.1)THEN
                INDX=(K2+1)*N1+K1+2
              ELSE
                INDX=(K3+1)*N1*N2+K1+2
              END IF
            ELSE
              INDX=(K3+1)*N1*N2+(K2+1)*N1+1
            END IF
          END IF
          T(8)=RAM(ISH+2*IRAM(ISH)+INDX)
C
          DO 70 ISM=1,8
            IF(T(ISM).LT.0.) THEN
              IUF=0
            END IF
   70     CONTINUE
          DO 75 ISM=3,5
            IF(Y(ISM).LT.0.) THEN
              IUF=0
            END IF
   75     CONTINUE
C
          IF(IUF.EQ.1) THEN
C           Trilinear interpolation of travel time
            FAK1=(Y(3)-O1)/D1-K1
            FAK2=(Y(4)-O2)/D2-K2
            FAK3=(Y(5)-O3)/D3-K3
            TAU=T(1)+(T(2)-T(1))*FAK1+(T(3)-T(1)+
     *          (T(1)-T(2)-T(3)+T(4))*FAK1)*FAK2
            TT= T(5)+(T(6)-T(5))*FAK1+(T(7)-T(5)+
     *          (T(5)-T(6)-T(7)+T(8))*FAK1)*FAK2
            TAU=TAU+(TT-TAU)*FAK3
            TAU=Y(1)+TAU
C
            IPUF=IPUF+1
            H11=YI(9)
            H21=YI(10)
            H31=YI(11)
            H13=YLI(1)*YI(6)
            H23=YLI(1)*YI(7)
            H33=YLI(1)*YI(8)
            H12=H23*H31-H33*H21
            H22=H33*H11-H13*H31
            H32=H13*H21-H23*H11
            V1=YLI(4)/YLI(1)**2
            V2=YLI(5)/YLI(1)**2
            E22=-2*H23*V2+H23*H23*(H13*V1+H23*V2)
C
            IF(H21.LT.0.1) THEN
              H21=0.1
            END IF
            C11=C(1)*H21**2
            C22=C(6)/H21**2
            C21=C(4)-C(6)*E22/H21**2
C
C           Writing the output data
            Y0=SQRT(C11/C22)
            R0=-C21/C22
            SIGMA=SQRT(Y0/C22)
            WRITE(LU3,'(1(G16.6,1X))') R0
            WRITE(LU4,'(1(G16.6,1X))') Y0
            WRITE(LU5,'(1(G16.6,1X))') SIGMA
            WRITE(LU6,'(1(G16.6,1X))') TAU
          ELSE
            IPUF=IPUF+1
            WRITE(LU3,'(I7,A)') 1,'*'
            WRITE(LU4,'(I7,A)') 1,'*'
            WRITE(LU5,'(I7,A)') 1,'*'
            WRITE(LU6,'(I7,A)') 1,'*'
          END IF
        ELSE
          IPUF=IPUF+1
          WRITE(LU3,'(I7,A)') 1,'*'
          WRITE(LU4,'(I7,A)') 1,'*'
          WRITE(LU5,'(I7,A)') 1,'*'
          WRITE(LU6,'(I7,A)') 1,'*'
        END IF
C
        ISRCO=ISRC
      GO TO 40
C
  100 CONTINUE
C
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'forms.for'
C     forms.for
C
C=======================================================================
C