SUBROUTINE OUT(X,Y,DERY,IHLF,NDIM,PRMT)
      DIMENSION Y(18),DERY(18),UN(3),POLD(3),PNEW(3),YO(3),
     1    DEP(6),PRMT(5)
      DIMENSION DN(3,3),DNG(3,2),XK(3),XKG(3),YK(3),YKG(3),
     1TM(2)
      COMMON /AUXI/  IANI(20),INTR,INT1,IOUT,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(4),NINT,
     1    XINTA
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1
      COMMON /RAY2/  DRY(3,400)
      COMMON /ZERO/ RNULL
      COMMON /VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP
      COMMON /DYN/XDYN(3,3),ydyn(3,3)
      common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18)
      COMMON /KM/KMAH,SPROLD,TSTOLD
C
      N=N+1
      if(ind.eq.10)go to 25
      NTR=20
      IF(N.GE.400) GOTO 152
      if(x.ge.tmax)go to 155
C
C  CHECK THE POSITION OF THE RAY WITH RESPECT TO BOUNDARIES OF THE MODEL
C
      NTR=21
      IF(N.EQ.1) GOTO 13
      NTR=22
      IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2).
     1GT.BRD(4)) GOTO 20
C
C  CHECK WETHER THE RAY CROSSED THE BOREHOLE
C
      IF(ICOD.EQ.0)GO TO 4
      IF(ICOD.GT.IREF)GO TO 4
      AUX1=(AY(2,N-1)-XVSP)*XNRM+(AY(3,N-1)-YVSP)*YNRM
      AUX2=(Y(1)-XVSP)*XNRM+(Y(2)-YVSP)*YNRM
      IF(AUX1*AUX2.GT.0.)GO TO 4
C
C  THE RAY CROSSED THE VERTICAL PLANE PERPENDICULAR TO THE LINE
C  SOURCE-BOREHOLE AND CONTAINING THE BOREHOLE
C  DETERMINATION OF THE POINT OF INTERSECTION
C
      TR=X
      XRST=DERY(1)
      YRST=DERY(2)
      AUX1=XVSP-Y(1)+XRST*TR
      AUX2=YVSP-Y(2)+YRST*TR
      X=(XNRM*AUX1+YNRM*AUX2)/(XNRM*XRST+YNRM*YRST)
      K=6
      IF(NDER.GT.1)K=18
      DO 2 I=1,K
    2 Y(I)=Y(I)+DERY(I)*(X-TR)
      CALL FCT(X,Y,DERY)
      IND=43
    4 CONTINUE
C
C  CHECK WETHER THE RAY CROSSED AN INTERFACE
C
      NTR=23
      INTR=LAY+1
      CALL DISC (Y,DEP)
      ZPL=DEP(1)
      INTR=LAY
      CALL DISC (Y,DEP)
      ZPU=DEP(1)
      NTR=24
      IF(ZPL.LE.ZPU) THEN
C
C  TWO NEIGHBOURHOOD INTERFACES CROSS EACH OTHER
C
      WRITE(LOU,'(A)') ' MAUX, LAY, IND,IND1,X,Y, Z,ZU,ZL'
      WRITE(LOU,'(4I5,5F12.5,/)')MAUX,LAY,IND,IND1,Y(1),Y(2),Y(3),ZPU,
     1ZPL
      GOTO 150
      END IF
      NTR=25
      is(3,iref)=lay
      IF(Y(3).LT.ZPU.OR.Y(3).GT.ZPL) GOTO 30
      NTR=27
C
C  THE RAY DID NOT CROSS AN INTERFACE
C
C
      IF(IND.EQ.1.OR.IND.EQ.2)GO TO 148
      IF(IND.EQ.30.OR.IND.EQ.31)GO TO 148
   13 CONTINUE
      AY(1,N)=X
      xold=x
      DO 15 I=1,6
        AY(I+1,N)=Y(I)
        yold(i)=y(i)
        dold(i)=dery(i)
   15 continue
      DRY(1,N)=DERY(1)
      DRY(2,N)=DERY(2)
      DRY(3,N)=DERY(3)
      IF(NDER.GT.1)THEN
C
C     DETERMINATION OF KMAH INDEX
C
        IF(N.GT.1)THEN
          SPR1=Y(8)*Y(12)-Y(9)*Y(11)
          SPR2=Y(9)*Y(10)-Y(7)*Y(12)
          SPR3=Y(7)*Y(11)-Y(8)*Y(10)
          SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3)
          IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1
          SPROLD=SPR
          TSTOLD=X
        END IF
        DO 14 I=1,12
          yold(I+6)=Y(I+6)
          dold(i+6)=dery(i+6)
   14   CONTINUE
      END IF
      NTR=99
      IF(IND.EQ.43)GO TO 25
C
C  WRITING IN AY FIELD
C
      CALL PARDIS (Y,1)
      RETURN
C
C  THE RAY CROSSED ONE OF THE VERTICAL BOUNDARIES OF THE MODEL
C
   20 IF(Y(1).LT.BRD(1)) IND=1
      IF(Y(1).GT.BRD(2)) IND=2
      IF(Y(2).LT.BRD(3)) IND=3
      IF(Y(2).GT.BRD(4)) IND=4
      IF(IND.EQ.3.OR.IND.EQ.4) GOTO 23
      AUX=(BRD(IND)-AY(2,N-1))/(Y(1)-AY(2,N-1))
      TR=X
      X=AY(1,N-1)+AUX*(X-AY(1,N-1))
      K=6
      IF(NDER.GT.1)K=18
      DO 21 I=1,K
   21 Y(I)=Y(I)+DERY(I)*(X-TR)
      Y(1)=BRD(IND)
      AY(1,N)=X
      DO 3 I=1,6
   3  AY(I+1,N)=Y(I)
      GO TO 4
  23  AUX=(BRD(IND)-AY(3,N-1))/(Y(2)-AY(3,N-1))
      TR=X
      X=AY(1,N-1)+AUX*(X-AY(1,N-1))
      K=6
      IF(NDER.GT.1)K=18
      DO 16 I=1,K
   16 Y(I)=Y(I)+DERY(I)*(X-TR)
      Y(2)=BRD(IND)
      AY(1,N)=X
      DO 17 I=1,6
  17  AY(I+1,N)=Y(I)
      IF(IND.EQ.3) IND=30
      IF(IND.EQ.4) IND=31
      GO TO 4
C
C  THE RAY CROSSED AN INTERFACE
C
   30 CONTINUE
      xo=xold
      do 35 i=1,3
        yo(i)=yold(i)
   35 continue
      IF(IREF.LE.1.OR.KC.GT.0) GOTO 39
      IF(N-KINT(IREF-1).EQ.2) THEN
      IND=9
      IND1=LAY
      GO TO 25
      END IF
C
C  WHICH INTERFACE WAS CROSSED
C
   39 INTR=LAY
      IF(Y(3).GE.ZPL) INTR=LAY+1
C
C  DETERMINATION OF THE POINT OF INTERSECTION WITH THE INTERFACE
C
      CALL ROOT(XO,YO,X,Y,XINT)
      NTR=10
      IND=0
      IAC=IAC+1
      AUX=XINT-XOLD
      IF(IRT.EQ.2)
     1WRITE(LOU,'(A,/,2F10.5,/)') 'XINT,XINTA',XINT,XINTA
      IF(IAC.EQ.1)GO TO 37
      IF(IAC.EQ.30)GO TO 153
      IF(ABS(XINT-XINTA).LT.RNULL)GO TO 38
      IF(ABS(AUX).LT.RNULL)GO TO 38
   37 XINTA=XINT
      PRMT(3)=AUX
      N=N-1
      PRMT(5)=1.1
      IF(IRT.EQ.2)
     1WRITE(LOU,'(A,/,4F10.5,I5,/)') ' XOLD,X,XINT,DT,IAC',
     2XOLD,X,XINT,PRMT(3),IAC
      RETURN
C
C     TIME CORRESPONDING TO THE POINT OF INCIDENCE OF THE RAY ON THE
C     INTERFACE WAS DETERMINED. THE QUANTITIES OF RAY TRACING AND
C     DYNAMIC RAY TRACING AT THE POINT OF INCIDENCE WILL BE DETERMINED
C
   38 CONTINUE
      IND1=INTR
      IF(IRT.GE.1)WRITE(LOU,'(A,I10)') '  IAC=',IAC
      IAC=0
      xnew=x
      x=xint
      kdim=6
      if(nder.gt.1)kdim=18
      do 40 i=1,kdim
        ynew(i)=y(i)
        dnew(i)=dery(i)
   40 continue
      call approx(x,y,dery,kdim)
      CALL DISC(Y,DEP)
      y(3)=DEP(1)
      AY(1,N)=XINT
      DO 63 I=1,6
        AY(I+1,N)=y(I)
   63 CONTINUE
      CALL PARDIS(Y,1)
      CALL FCT(X,Y,DERY)
      DO 64 I=1,3
        DRY(I,N)=DERY(I)
   64 CONTINUE
      IF(NDER.GT.1)THEN
        IF(N.GT.1)THEN
          SPR1=Y(8)*Y(12)-Y(9)*Y(11)
          SPR2=Y(9)*Y(10)-Y(7)*Y(12)
          SPR3=Y(7)*Y(11)-Y(8)*Y(10)
          SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3)
          IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1
          SPROLD=SPR
          TSTOLD=X
        END IF
        DO 62 I=1,3
          XDYN(I,1)=Y(I+6)
          XDYN(I,2)=Y(I+9)
          XDYN(I,3)=DERY(I)
          YDYN(I,1)=Y(I+12)
          YDYN(I,2)=Y(I+15)
          ydyn(i,3)=dery(i+3)
          XK(I)=DERY(I)
          YK(I)=DERY(I+3)
   62   CONTINUE
      END IF
C
C  DETERMINATION OF THE UNIT NORMAL TO THE INTERFACE AT THE POINT OF
C  INCIDENCE
C
      CALL DISC (Y,DEP)
      DFX=DEP(2)
      DFY=DEP(3)
      ROO=1+DFX*DFX+DFY*DFY
      UN3=ROO**(-0.5)
      UN(3)=-UN3
      UN(2)=UN3*DFY
      UN(1)=UN3*DFX
C
C  SCALAR PRODUCT OF GROUP VELOCITY WITH UNIT NORMAL VECTOR
C
      PN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3)
C
C  ORIENTATION OF THE UNIT NORMAL SO THAT IT POINTS TO THE MEDIUM IN
C  WHICH THE ENERGY OF INCIDENT WAVE PROPAGATES
C
      IF(PN.LT.0) GOTO 120
      UN(1)=-UN(1)
      UN(2)=-UN(2)
      UN(3)=-UN(3)
      UN3=-UN3
  120 CONTINUE
      IF(NDER.GT.1)THEN
        XIINC=Y(4)*UN(1)+Y(5)*UN(2)+Y(6)*UN(3)
        UN33=UN3*UN3*UN3
        DFXX=DEP(4)
        DFXY=DEP(5)
        DFYY=DEP(6)
        TM(1)=0.
        TM(2)=0.
        DO 126 J=1,2
          AUX=0.
          DO 125 I=1,3
            TM(J)=TM(J)+UN(I)*XDYN(I,J)
            AUX=AUX+UN(I)*XK(I)
  125     CONTINUE
          TM(J)=-TM(J)/AUX
  126   CONTINUE
C
C     DETERMINATION OF DERIVATIVES OF THE UNIT NORMAL TO INTERFACE
C
        AUX1=1.+DFX*DFX
        AUX2=1.+DFY*DFY
        AUX3=DFX*DFY
        DN(1,1)=(DFXX*AUX2-AUX3*DFXY)*UN33
        DN(1,2)=(DFXY*AUX2-AUX3*DFYY)*UN33
        DN(1,3)=0.
        DN(2,1)=(DFXY*AUX1-AUX3*DFXX)*UN33
        DN(2,2)=(DFYY*AUX1-AUX3*DFXY)*UN33
        DN(2,3)=0.
        DN(3,1)=(DFX*DFXX+DFY*DFXY)*UN33
        DN(3,2)=(DFX*DFXY+DFY*DFYY)*UN33
        DN(3,3)=0.
        DO 129 L=1,3
          DO 128 M=1,2
            DNG(L,M)=0.
            DO 127 K=1,3
              DNG(L,M)=DNG(L,M)+DN(L,K)*(XDYN(K,M)+XK(K)*TM(M))
  127       CONTINUE
  128     CONTINUE
  129   CONTINUE
      END IF
      NTR=20
      IF(KRE.EQ.1) GOTO 24
      NTR=30
      IF(KRE.EQ.0) GOTO 130
C
C  MULTIPLE REFLECTED WAVE
C
      IF((IREF+1).GT.KRE.AND.INTR.EQ.INT1) GOTO 151
C
C  TERMINATION OF THE RAY AT AN INNER INTERFACE OR AT THE FREE SURFACE
C  OR AT THE BOTTOM OF THE MODEL
C
      IF((IREF+1).GT.KRE) IRR=IREF
      NTR=35
      IF((IREF+1).GT.KRE) GOTO 148
C
C  NCD : INDICATOR WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE
C        WITH RESPECT TO THE CODE OF THE WAVE
C
      NCD=CODE(IREF+1,1)-CODE(IREF,1)
C
C  NCD1 : INDICATOR FOR THE TYPE OF THE WAVE IN THE NEW LAYER
C
      NCD1=CODE(IREF+1,2)-CODE(IREF,2)
      NTR=40
      IF(KC.GT.0.AND.IREF.EQ.1.AND.INTR.EQ.LAY) GOTO 151
      NTR=50
      IF(KC.LT.0.AND.IREF.EQ.1.AND.INTR.NE.LAY) GOTO 151
      NTR=60
      IF(IREF.EQ.1) GOTO 170
      NTR=70
C     IF(INTR.EQ.INT1.AND.NCD.NE.0.OR.NCD1.NE.0) GOTO 151
      IF(INTR.EQ.INT1.AND.NCD.NE.0)GOTO 151
      NTR=75
      IF(INTR.EQ.INT1.AND.NCD1.NE.0)GOTO 151
C
C  REFLECTION OR TRANSMISSION OF THE RAY
C
      NTR=80
      IF(INTR.NE.INT1.OR.NCD.NE.0.OR.NCD1.NE.0) GOTO 170
C
C  COMPOUND ELEMENT OF THE RAY
C
      IREFR=1
      KINT(IREF)=0
      IREF=IREF+1
      NCD=CODE(IREF+1,1)-CODE(IREF,1)
      NCD1=CODE(IREF+1,2)-CODE(IREF,2)
C
C  TERMINATION AT AN INNER INTERFACE
C
      IRR=IREF
      NTR=90
      IF((IREF+1).GT.KRE) GOTO 24
      NTR=100
      GOTO 170
C
C  REFRACTED WAVE
C
C  ORDINARY TERMINATION AT UPPER BOUNDARY
C
  130 NTR=110
      IF(INTR.EQ.LAY.AND.LAY.EQ.1) GOTO 148
      NCD=1
      NCD1=0
C
C  SPECIFICATION OF THE LAYER OF THE GENERATED WAVE
C
  170 INT1=INTR
      IRR=IREF
      IREF=IREF+1
      LAY1=LAY
      NTR=120
      IF(NCD.EQ.0) GOTO 200
      NTR=130
      IF(INTR.EQ.LAY) GOTO 190
C
C  TRANSMISSION AT THE LOWER INTERFACE
C
      NTR=140
      IF(NCD.LT.0) GOTO 151
C
C  ORDINARY TERMINATION AT LOWER BOUNDARY
C
      NTR=150
      IF(KRE.LE.1.AND.INTR.EQ.NINT) GOTO 148
      NTR=160
      IF(INTR.EQ.NINT) GOTO 151
      LAY=LAY+1
      NTR=170
      GOTO 200
C
C  TRANSMISSION AT THE UPPER INTERFACE
C
  190 NTR=180
      IF(NCD.GT.0.AND.KRE.GT.1) GOTO 151
      NTR=190
      IF(KRE.LE.1.AND.LAY.EQ.1) GOTO 24
      NTR=200
      IF(LAY.EQ.1) GOTO 151
      LAY=LAY-1
C
C  TRANSFORMATION OF THE SLOWNESS VECTOR
C
  200 IF(INTR.EQ.NINT.AND.MDIM.NE.0) GOTO 154
      DO 210 I=1,3
      POLD(I)=Y(3+I)
  210 CONTINUE
      ITRANS=0
C
C  CHECK WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE
C
      IF(KC.EQ.0) GOTO 217
      IF(CODE(IREF,1).EQ.CODE(IREF-1,1)) GOTO 218
  217 ITRANS=1
  218 CALL TRANSL (Y,POLD,PNEW,UN,ITRANS,1)
      IF(IND.EQ.9) GOTO 25
      DO 220 I=1,3
      Y(I+3)=PNEW(I)
  220 CONTINUE
      IF(NDER.GT.1)THEN
        CALL FCT(X,Y,DERY)
        DO 225 I=1,3
          XKG(I)=DERY(I)
          YKG(I)=DERY(I+3)
  225   CONTINUE
        XI=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3)
        DXTN=XKG(1)*UN(1)+XKG(2)*UN(2)+XKG(3)*UN(3)
        DO 228 M=1,2
          AUX=XKG(1)*DNG(1,M)+XKG(2)*DNG(2,M)+XKG(3)*DNG(3,M)
          DO 227 I=1,3
            AUXX=(DNG(I,M)-UN(I)*AUX/DXTN)*(XI-XIINC)
            DO 226 K=1,3
              AUXX=AUXX-UN(I)*(XKG(K)*(YDYN(K,M)+YK(K)*TM(M))-
     1                  YKG(K)*(XDYN(K,M)+XK(K)*TM(M)))/DXTN
  226       CONTINUE
            XDYN(I,M)=XDYN(I,M)+(XK(I)-XKG(I))*TM(M)
            YDYN(I,M)=YDYN(I,M)+(YK(I)-YKG(I))*TM(M)+AUXX
  227     CONTINUE
  228   CONTINUE
        DO 229 I=1,3
          Y(I+6)=XDYN(I,1)
          Y(I+9)=XDYN(I,2)
          Y(I+12)=YDYN(I,1)
          Y(I+15)=YDYN(I,2)
  229   CONTINUE
      END IF
      PRMT(5)=2.
      RETURN
C
C  ORDINARY TERMINATION OF THE RAY
C
 148  IF(IND.EQ.1.OR.IND.EQ.2.OR.IND.EQ.30.OR.IND.EQ.31) GOTO 25
  24  IND=INTR+100
      IF(KRE.LE.1) IRR=IREF
      IF(INTR.EQ.1) IND=3
      IF(INTR.EQ.NINT) IND=4
      KINT(IRR)=N
      IF(IND.NE.3) GOTO 25
C
C  COMPUTATION OF REFLECTED WAVES AT THE EARTH SURFACE FOR COEFFICIENTS
C  OF CONVERSION
C
      DINC=UN(1)*DRY(1,N)+UN(2)*DRY(2,N)+UN(3)*DRY(3,N)
      DINC=DINC/SQRT(DRY(1,N)*DRY(1,N)+DRY(2,N)*DRY(2,N)
     1     +DRY(3,N)*DRY(3,N))
      DINC=ACOS(DINC)
      IF(MREG.EQ.0.AND.MDIM.NE.0) THEN
        POLD(1)=Y(4)
        POLD(2)=Y(5)
        POLD(3)=Y(6)
        IREF=IREF+1
        CALL TRANSL(Y,POLD,PNEW,UN,0,1)
        IREF=IREF-1
      END IF
      GOTO 25
  150 IND=20
      GOTO 25
  151 IND=8
      GOTO 25
  152 IND=13
      GOTO 25
  153 IND=39
      IAC=0
      GOTO 25
  154 IND=15
      GOTO 25
  155 IND=12
   25 PRMT(5)=3.0
      IF(IND.NE.3.AND.IND.NE.43) RETURN
      ITYPE=CODE(IREF,2)
      CALL PARDIS(Y,2)
      RETURN
      END