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: 5.70
C Date: 2003, April 4
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: zacek@karel.troja.mff.cuni.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     RAYALL='string'... Name of the input file CRT-R with the
C             quantities stored along rays, see
C             C.R.T.5.5.1.
C             Default: RAYALL='r01.out' (usually sufficient)
C     INIALL='string'... Name of the input file CRT-I with the
C             quantities at the initial points of rays, corresponding
C             to file RAYALL, see
C             C.R.T.6.1.
C             Default: INIALL='r01i.out' (usually sufficient)
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 Names of output files:
C     GBOUT1='string'... Name of the output file containing optimum
C             value Y0=SQRT(C11/C22)
C             of the imaginary part of parameter M0.
C             Default: GBOUT1='y0.out'
C     GBOUT2='string'... Name of the output file containing optimum
C             value R0=-C21/C22 of the real part of parameter M0.
C             Default: GBOUT2='r0.out'
C     GBOUT3='string'... Name of the output file containing initial
C             value SIGMA=SQRT(Y0/C22) of a standard deviation for
C             smoothing.
C             Default: GBOUT3='sigma.out'
C Data specifying the dimensions of the input 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
C Input unformatted file RAYALL:
C     See the description within source code file 'writ.for'.
C     Description of files CRT-R
C
C Input unformatted file INIALL:
C     See the description within source code file 'writ.for'.
C     Description of files CRT-I
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL AP00
C     AP00... File 'ap.for'.
      EXTERNAL AP28
C     AP28... File 'ap.for'.
      EXTERNAL RARRAY
C     RARRAY...File 'forms.for'.
      EXTERNAL RARRAI
C     RARRAI...File 'forms.for'.
      EXTERNAL ERROR
C     ERROR...File 'error.for'.
      EXTERNAL RSEP1
C     RSEP1...File 'sep.for'.
      EXTERNAL RSEP3T
C     RSEP3T...File 'sep.for'.
      EXTERNAL RSEP3I
C     RSEP3I...File 'sep.for'.
      EXTERNAL RSEP3R
C     RSEP3R...File 'sep.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
      INCLUDE 'sep.inc'
C     sep.inc
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL O1,O2,O3,D1,D2,D3,NORM0,NORM1,NORM2,FAK1,FAK2,FAK3,DET
      REAL TAU,TT,FX,FFX
      REAL C(10),DC(7),SDC(7),SSDC(3),FN1(7),FFN1(3),T(8)
      REAL H11,H12,H13,H21,H22,H23,H31,H32,H33,V1,V2,E11,E12,E22
      REAL CA,CB,CC11,CC12,CC22,C11,C12,C21,C22,Y0,R0,SIGMA
      INTEGER LU1,LU2,LU3,LU4,LU5,LU6,LU7,LU8,LU9,LU10,LU11,
     *        IPUF,INR,IP,IRR,IRRN
      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=21,LU2=22,LU3=23,LU4=24,LU5=25,LU6=26,LU7=27)
      PARAMETER (LU8=28,LU9=29,LU10=30,LU11=31)
      CHARACTER*80 FILSEP,RAYALL,INIALL,NUM,MTT,AMP
      CHARACTER*80 GBOUT1,GBOUT2,GBOUT3,GBOUT4,GBOUT5
      EQUIVALENCE (IRAM,RAM)
      REAL UNDEF,YMIN
      PARAMETER (UNDEF=-9.999E9,YMIN=1E-9)
      PARAMETER (INR=26)
      PARAMETER (IRR=181)
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 files with computed rays and
C     multivalued travel times
      CALL RSEP3T('RAYALL',RAYALL,'r01.out')
      CALL RSEP3T('INIALL',INIALL,'r01i.out')
      CALL RSEP3T('NUM',NUM,'mtt-num.out')
      CALL RSEP3T('MTT',MTT,'mtt-tt.out')
      CALL RSEP3T('AMP',AMP,'mtt-ap.out')
C
C     Reading filenames of the output files
      CALL RSEP3T('GBOUT1',GBOUT1,'y0.out')
      CALL RSEP3T('GBOUT2',GBOUT2,'r0.out')
      CALL RSEP3T('GBOUT3',GBOUT3,'sigma.out')
C      CALL RSEP3T('GBOUT4',GBOUT4,'time.out')
C      CALL RSEP3T('GBOUT5',GBOUT5,'h21.out')
C
      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.)
      OPEN(LU2,FILE=RAYALL,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=INIALL,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU4,FILE=NUM,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU5,FILE=MTT,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU6,FILE=AMP,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU7,FILE=GBOUT1,FORM='FORMATTED')
      OPEN(LU8,FILE=GBOUT2,FORM='FORMATTED')
      OPEN(LU9,FILE=GBOUT3,FORM='FORMATTED')
C      OPEN(LU10,FILE=GBOUT4,FORM='FORMATTED')
C      OPEN(LU11,FILE=GBOUT5,FORM='FORMATTED')
      N1N2N3=N1*N2*N3
      ISH=2*N1N2N3+1
      FX=0.
      FFX=0.
      IPUF=0
      IRRN=1
C
      CALL RARRAI(LU4,' ','FORMATTED',.FALSE.,0,N1N2N3,IRAM)
      IRAM(N1N2N3+1)=0
      DO 10 IGRD=1,N1N2N3
        IRAM(N1N2N3+IGRD+1)=IRAM(N1N2N3+IGRD)+IRAM(IGRD)
   10 CONTINUE
      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 ram.inc.
      END IF
      CALL RARRAY(LU5,' ','FORMATTED',.FALSE.,0.,IRAM(ISH),
     *            RAM(ISH+1))
      CALL RARRAY(LU6,' ','FORMATTED',.FALSE.,1.,IRAM(ISH),
     *            RAM(ISH+IRAM(ISH)+1))
      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     Loop for the points of rays
   40 CONTINUE
C       Reading the results of the complete ray tracing
        CALL AP00(0,LU2,LU3)
        IF(IWAVE.LT.1)THEN
C         End of rays
          IRRN=IRRN-1
          IF (IPUF.LT.INR) THEN
            IPUF=INR-IPUF
            DO 41 IP=1,IPUF
C              WRITE(LU7,'(1(G16.6,1X))') UNDEF
C              WRITE(LU8,'(1(G16.6,1X))') UNDEF
C              WRITE(LU9,'(1(G16.6,1X))') UNDEF
C              WRITE(LU10,'(1(G16.6,1X))') UNDEF
C              WRITE(LU11,'(1(G16.6,1X))') UNDEF
              WRITE(LU7,'(I7,A)') 1,'*'
              WRITE(LU8,'(I7,A)') 1,'*'
              WRITE(LU9,'(I7,A)') 1,'*'
   41       CONTINUE
          ELSE
C           GBOPT-03
            CALL ERROR('GBOPT-03: Too low parameter INR')
          END IF
          IF (IRRN.LT.IRR) THEN
            IRRN=INR*(IRR-IRRN)
            DO 42 IP=1,IRRN
C              WRITE(LU7,'(1(G16.6,1X))') UNDEF
C              WRITE(LU8,'(1(G16.6,1X))') UNDEF
C              WRITE(LU9,'(1(G16.6,1X))') UNDEF
C              WRITE(LU10,'(1(G16.6,1X))') UNDEF
C              WRITE(LU11,'(1(G16.6,1X))') UNDEF
              WRITE(LU7,'(I7,A)') 1,'*'
              WRITE(LU8,'(I7,A)') 1,'*'
              WRITE(LU9,'(I7,A)') 1,'*'
   42       CONTINUE
          END IF
          GO TO 100
        END IF
        IF (IPT.LE.1) THEN
C         New ray
          IF (IRAY.GT.IRRN) THEN
            IRRN=INR*(IRAY-IRRN)
            DO 45 IP=1,IRRN
C              WRITE(LU7,'(1(G16.6,1X))') UNDEF
C              WRITE(LU8,'(1(G16.6,1X))') UNDEF
C              WRITE(LU9,'(1(G16.6,1X))') UNDEF
C              WRITE(LU10,'(1(G16.6,1X))') UNDEF
C              WRITE(LU11,'(1(G16.6,1X))') UNDEF
              WRITE(LU7,'(I7,A)') 1,'*'
              WRITE(LU8,'(I7,A)') 1,'*'
              WRITE(LU9,'(I7,A)') 1,'*'
   45       CONTINUE
          END IF
          IRRN=IRAY+1
          IF (IRAY.GT.1) THEN
            IF (IPUF.LT.INR) THEN
              IPUF=INR-IPUF
              DO 48 IP=1,IPUF
C                WRITE(LU7,'(1(G16.6,1X))') UNDEF
C                WRITE(LU8,'(1(G16.6,1X))') UNDEF
C                WRITE(LU9,'(1(G16.6,1X))') UNDEF
C                WRITE(LU10,'(1(G16.6,1X))') UNDEF
C                WRITE(LU11,'(1(G16.6,1X))') UNDEF
                WRITE(LU7,'(I7,A)') 1,'*'
                WRITE(LU8,'(I7,A)') 1,'*'
                WRITE(LU9,'(I7,A)') 1,'*'
   48         CONTINUE
            ELSE
C             GBOPT-04
              CALL ERROR('GBOPT-04: Too low parameter INR')
            END IF
          END IF
          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
        END IF
C       Numerical quadrature:
        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)
        CALL AP28(7,SDC,1,1,0.01,FX,ISRF1,NFN1,IFN1,FN1,NFN2,IFN2,DC)
        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.NE.0.) THEN
C          DC(4)=Y(12)-((Y(20)*C(10)-Y(24)*C(9))*C(4)+(Y(24)*C(6)-
C     *          Y(20)*C(9))*C(7))/DET
C          DC(5)=Y(13)-((Y(21)*C(10)-Y(25)*C(9))*C(4)+(Y(25)*C(6)-
C     *          Y(21)*C(9))*C(7))/DET
C          DC(6)=Y(16)-((Y(20)*C(10)-Y(24)*C(9))*C(5)+(Y(24)*C(6)-
C     *          Y(20)*C(9))*C(8))/DET
C          DC(7)=Y(17)-((Y(21)*C(10)-Y(25)*C(9))*C(5)+(Y(25)*C(6)-
C     *          Y(21)*C(9))*C(8))/DET
C
          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)
        END IF
        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
        CALL AP28(3,SSDC,1,1,0.01,FFX,ISRF1,
     *                                  NFFN1,IFFN1,FFN1,NFFN2,IFFN2,DC)
        C(1)=SSDC(1)
        C(2)=SSDC(2)
        C(3)=SSDC(3)
C
C       Calculation of the mean arrival time at the reciever:
        IUF=1
        K1=INT((Y(3)-O1)/D1)
        K2=INT((Y(4)-O2)/D2)
        K3=INT((Y(5)-O3)/D3)
C
        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
        IF (FLOAT(N1-1)*(FLOAT(N1-2)*D1-Y(3)).LT.0.) THEN
          IUF=0
        END IF
        IF (FLOAT(N2-1)*(FLOAT(N2-2)*D2-Y(4)).LT.0.) THEN
          IUF=0
        END IF
        IF (FLOAT(N3-1)*(FLOAT(N3-2)*D3-Y(5)).LT.0.) THEN
          IUF=0
        END IF
        IF (IUF.EQ.1) THEN
          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
C         Writing the output data:
          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
C          E11=-H13*V1-H13*V1+H13*H13*(H13*V1+H23*V2)
C          E12=-H13*V2-H23*V1+H13*H23*(H13*V1+H23*V2)
          E22=-2*H23*V2+H23*H23*(H13*V1+H23*V2)
C
C         CA=C11surf, CB=C22surf, CC=C21surf
C          H21=1.
          IF (H21.LT.0.1) THEN
            H21=0.1
          END IF
          C11=C(1)*H21**2
          C22=C(6)/H21**2
C          CC11=C(4)-C(6)*E22/H21**2
C          CC12=-C(6)*E22/H21**2
          C21=C(4)-C(6)*E22/H21**2
C          WRITE(LU7,'(1(G16.6,1X))') CA/Y(1)
C          WRITE(LU8,'(1(G16.6,1X))') CB/Y(1)**2
C          WRITE(LU9,'(1(G16.6,1X))') CC11/Y(1)
C          WRITE(LU10,'(1(G16.6,1X))') CC12/Y(1)
C          WRITE(LU11,'(1(G16.6,1X))') CC22/Y(1)
C
C          IF (Y0.LT.YMIN) THEN
C            Y0=UNDEF
C          END IF
C          IF ((C11/C22).LT.0.) THEN
C            Y0=UNDEF
C          ELSE
            Y0=SQRT(C11/C22)
C          END IF
C          R011=-CC11/CB
C          R012=-CC12/CB
          R0=-C21/C22
          SIGMA=SQRT(Y0/C22)
          WRITE(LU7,'(1(G16.6,1X))') Y0
          WRITE(LU8,'(1(G16.6,1X))') R0
          WRITE(LU9,'(1(G16.6,1X))') SIGMA
C          WRITE(LU10,'(1(G16.6,1X))') R012
C          WRITE(LU11,'(1(G16.6,1X))') H21
        ELSE
          IPUF=IPUF+1
C          WRITE(LU7,'(1(G16.6,1X))') UNDEF
C          WRITE(LU8,'(1(G16.6,1X))') UNDEF
C          WRITE(LU9,'(1(G16.6,1X))') UNDEF
C          WRITE(LU10,'(1(G16.6,1X))') UNDEF
C          WRITE(LU11,'(1(G16.6,1X))') UNDEF
          WRITE(LU7,'(I7,A)') 1,'*'
          WRITE(LU8,'(I7,A)') 1,'*'
          WRITE(LU9,'(I7,A)') 1,'*'
        END IF
C
      GO TO 40
C
  100 CONTINUE
C      DO 120 ISM=1,IRR
C        TIME=0.
C        IF (ISM.LE.11) THEN
C          H21=SIN(10*0.034907/2)
C        END IF
C        IF (ISM.GT.11) THEN
C          H21=SIN(ISM*0.017165)
C        END IF
C        IF (ISM.GE.171) THEN
C          H21=SIN(172*0.034907/2)
C        END IF
C        DO 110 IP=1,INR
C          TIME=TIME+0.1
C          WRITE(LU10,'(1(G16.6,1X))') TIME
C          WRITE(LU11,'(1(G16.6,1X))') H21**2
C  110   CONTINUE
C  120 CONTINUE
C
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
      CLOSE(LU7)
      CLOSE(LU8)
      CLOSE(LU9)
C      CLOSE(LU10)
C      CLOSE(LU11)
      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.forC
C
C=======================================================================
C