      SUBROUTINE RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX,
     1ASTART,ASTEP,AFIN,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCD)
C
C     TWO-POINT RAY TRACING
C
      COMPLEX AMPX1,AMPX2,AMPY1,AMPY2,AMPZ1,AMPZ2,APX,APY,APZ
      DIMENSION JC(50,2),YDD(2),DEP(6)
      DIMENSION TIME(500),XCOOR(500),ZCOOR(500),INDI(500),AMPX1(500),
     1AMPY1(500),AMPZ1(500),AMPX2(500),AMPY2(500),AMPZ2(500),
     2p(500,3),pol(500,3),pol1(500,3),APX(2),APY(2),APZ(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 /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST,
     1xprf,yprf
      COMPLEX PS,CKMAH
      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/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP
      COMMON/DYN/XDYN(3,3),ydyn(3,3)
      COMMON/KM/KMAH,SPROLD,TSTOLD
C
      ITER=0
      iprof=0
      indr1=0
      JNULL=0
      II=0
      DD=0.
C
      AA=ASTART-ASTEP
      BMIN1=BMIN
      REPS1=.05*REPS
      INDEX=0
      INUM=0
      IK1=0
      ICLS=0
      ICR=0
      INDS=ISOUR
C
C     LOOP IN DECLINATION, FROM VALUE ASTART TO AFIN, WITH THE STEP
C     ASTEP
C
    1 AA=AA+ASTEP
      PNEW=AA
      IF(ASTEP.GT.0..AND.AA.GT.AFIN)GO TO 99
      IF(ASTEP.LT.0..AND.AA.LT.AFIN)GO TO 99
    4 IND=INDS
      CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,AA,PAZM,r,XX,YY,Zz,T,
     1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1)
      IF(IND.EQ.14.OR.IND.EQ.20)RETURN
      IF(ITPR.EQ.43)r=Zz
      IF(NDST.EQ.0.AND.LNDST.EQ.0)THEN
        BMIN=BMIN1
        GO TO 1
      END IF
      IF(NDST.EQ.0)GO TO 65
      IF(IND.EQ.ITPR)rAX=r
      IF(IND.EQ.ITPR)PNEW=AA
      IF(IND.EQ.9)XCR=r
      IF(MOUT.GE.2)WRITE(LOU,100)IND,IND1,KMAH,r,T,AA
      IF(INUM.NE.0)GO TO 2
      dOLD=AA
      aold=pazm
      IOLD=IND
      IOLD1=IND1
      rOLD=r
      xold=xx
      yold=yy
      zold=zz
      TOLD=T
      INUM=1
      GO TO 1
C
C     PARAMETERS FOR THE PRECEDING RAY: 
C       DOLD (AA), AOLD (PAZM), ROLD (R), IOLD (IND)
C     PARAMETERS FOR THE NEW RAY: 
C       DNEW (AA), ANEW (PAZM), RNEW (R), INEW (IND)
C
    2 INEW=IND
      INEW1=IND1
      dNEW=AA
      anew=pazm
      rNEW=r
      xnew=xx
      ynew=yy
      znew=zz
      TNEW=T
      IF(INEW.EQ.ITPR.AND.IOLD.EQ.ITPR)GO TO 21
      IF(INEW.EQ.ITPR)GO TO 50
      IF(IOLD.EQ.ITPR)GO TO 55
      IF(INEW.EQ.9.AND.IOLD.NE.9.AND.IOLD.NE.8)GO TO 30
      IF(INEW.NE.9.AND.INEW.NE.8.AND.IOLD.EQ.9)GO TO 35
      indr1=ind1
      GO TO 3
   21 IF(IOLD1.NE.INEW1)IK1=2
      IF(IK1.EQ.2)GO TO 55
      indr1=ind1
      GO TO 40
C
C     NO ITERATIONS, TAKE A NEW RAY IN THE LOOP
C
    3 CONTINUE
      IOLD=INEW
      IOLD1=INEW1
      rOLD=rNEW
      xold=xnew
      yold=ynew
      zold=znew
      dOLD=dNEW
      aold=anew
      TOLD=TNEW
      GO TO 1
C
C     REGULAR CASE: IOLD=3, INEW=3
C
   40 rrNEW=rNEW
      xxnew=xnew
      yynew=ynew
      zznew=znew
      rrOLD=rOLD
      xxold=xold
      yyold=yold
      zzold=zold
      ddNEW=dNEW
      aanew=anew
      ddOLD=dOLD
      aaold=aold
      TTNEW=TNEW
      TTOLD=TOLD
      IBNEW=0
      IBOLD=0
   41 IF(rrNEW.GT.rrOLD.AND.DST(2).GT.DST(1))GO TO 46
      IF(rrNEW.LT.rrOLD.AND.DST(2).LT.DST(1))GO TO 46
C
C     REGULAR CASE: RRNEW.LE.RROLD, ITREND=-1 (REVERSE BRANCH)
C
      ITREND=-1
      DO 42 I=1,NDST
      NDD=NDST-I+1
      DD=DST(NDD)
      dr=rrOLD
      IF(IBOLD.EQ.1)dr=dr+REPS
      IF(DD.GE.dr)GO TO 42
      dr=rrNEW
      IF(IBNEW.EQ.1)dr=dr-REPS
      IF(DD.LT.dr.AND.IK1.NE.0)GO TO 6
      IF(DD.LT.dr)GO TO 3
      II=NDD
      GO TO 43
   42 CONTINUE
      IF(IK1.NE.0)GO TO 6
      GO TO 3
C
C     REGULAR CASE: RRNEW.GT.RROLD, ITREND=1 (REGULAR BRANCH)
C
   46 ITREND=1
      DO 44 I=1,NDST
      DD=DST(I)
      dr=rrOLD
      IF(IBOLD.EQ.1)dr=dr-REPS
      IF(DD.LE.dr)GO TO 44
      dr=rrNEW
      IF(IBNEW.EQ.1)dr=dr+REPS
      IF(DD.GT.dr.AND.IK1.NE.0)GO TO 6
      IF(DD.GT.dr)GO TO 3
      II=I
      GO TO 43
   44 CONTINUE
      IF(IK1.NE.0)GO TO 6
      GO TO 3
   43 d1=ddOLD
      a1=aaold
      d2=ddNEW
      a2=aanew
      x1=xxold
      y1=yyold
      z1=zzold
      x2=xxnew
      y2=yynew
      z2=zznew
      T1=TTOLD
      T2=TTNEW
C
C     I T E R A T I O N S
C
   60 continue
C
C     TRANSFORMATION OF COORDINATES OF A RECEIVER FROM CYLINDRICAL
C     TO CARTESIAN COORDINATES
C
      IF(ITPR.NE.43)THEN
        XD=Xprf+DD*COS(PROF(1))
        YD=Yprf+DD*SIN(PROF(1))
        YDD(1)=XD
        YDD(2)=YD
        INTR=1
        IF(ITPR.GT.100)INTR=ITPR-100
        CALL DISC(YDD,DEP)
        ZD=DEP(1)
      END IF
      IF(ITPR.EQ.43)THEN
        XD=XVSP
        YD=YVSP
        ZD=DD
      END IF
      DELX=XD-X1
      DELY=YD-Y1
      DELZ=ZD-Z1
      dr1=sqrt(delx*delx+dely*dely+delz*delz)
      DELX=XD-X2
      DELY=YD-Y2
      DELZ=ZD-Z2
      dr2=sqrt(delx*delx+dely*dely+delz*delz)
c
      ITER=0
      ISIGN=1
      IPR1=0
      IPR2=0
      GO TO 61
   68 rAX=r
      PAX=PNEW
   61 ITER=ITER+1
      IF(ITER.GT.ITMAX)GO TO 80
      ISIGN=-ISIGN
C
C     PREPARATION FOR ITERATIVE SOLUTION OF TWO-POINT RAY TRACING
C
      daux=0.5*(d1+d2)
      azaux=0.5*(a1+a2)
      if(method.eq.2.or.iter.eq.1.or.itpr.ne.ind)go to 69
C      
C     PARAXIAL RAY APPROXIMATION
C
      AUX1=XDYN(2,1)*XDYN(3,2)-XDYN(3,1)*XDYN(2,2)
      AUX2=XDYN(3,1)*XDYN(1,2)-XDYN(1,1)*XDYN(3,2)
      AUX3=XDYN(1,1)*XDYN(2,2)-XDYN(2,1)*XDYN(1,2)
      DET=AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3)
      IF(ABS(DET).LT..00001)GO TO 69
      AUX1=DELY*XDYN(3,2)-DELZ*XDYN(2,2)
      AUX2=DELZ*XDYN(1,2)-DELX*XDYN(3,2)
      AUX3=DELX*XDYN(2,2)-DELY*XDYN(1,2)
      aux11=(AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3))/DET
      AUX1=DELZ*XDYN(2,1)-DELY*XDYN(3,1)
      AUX2=DELX*XDYN(3,1)-DELZ*XDYN(1,1)
      AUX3=DELY*XDYN(1,1)-DELX*XDYN(2,1)
      aux22=(AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3))/DET
      if(mori.eq.0)pazm=pazm+aux11
      if(mori.ne.0)pazm=pazm+aux22
      if(mori.eq.0)pnew=pnew+aux22
      if(mori.ne.0)pnew=pnew+aux11
      go to 72
C
C     HALVING OF INTERVAL
C
   69 pnew=0.5*(d1+d2)
      pazm=0.5*(a1+a2)
C
C     INITIAL ANGLES FOR A NEW RAY WERE DETERMINED
C
   72 ind=inds
      NOLD=0
      rOLD=0.
      SPROLD=0.
      IF(MDIM.GE.1)NDER=2
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,PX,PY,PZ,XX,YY,Zz,
     1T,DT,AC)
      NDER=1
      XE=XX-Xprf
      YE=YY-Yprf
      r=SQRT(XE*XE+YE*YE)
      delx=xd-xx
      dely=yd-yy
      delz=zd-zz
      drs=sqrt(delx*delx+dely*dely+delz*delz)
      IF(IND.EQ.20)RETURN
      IF(ITPR.EQ.43)r=Zz
      IF(MOUT.GE.2)WRITE(LOU,101)IND,IND1,ITER,KMAH,DD,r,T,PNEW,PAZM
      if(mout.eq.4)write(lou,120)xd,yd,zd,xx,yy,zz
  120 format(1x,'(x,y,z) receiver',3F15.8/1x,'(xx,yy,zz) ray',3F15.8)
      IF(ICLS.NE.0)GO TO 70
      IF(IND.NE.ITPR)then
        d2=PNEW
        a2=pazm
        x2=xx
        y2=yy
        z2=zz
        dr2=drs
        GO TO 61
      end if
      IF(drs.LT.REPS)GO TO 65
      IF(dr2.LT.dr1)GO TO 63
      d2=PNEW
      a2=pazm
      x2=xx
      y2=yy
      z2=zz
      dr2=drs
      T2=T
      IPR2=1
      GO TO 68
   63 continue
      d1=PNEW
      a1=pazm
      x1=xx
      y1=yy
      z1=zz
      dr1=drs
      T1=T
      IPR1=1
      GO TO 68
   70 ICLS=0
C
C     SUCCESSFUL ITERATIONS
C
   65 INDEX=INDEX+1
      IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 800
      TIME(INDEX)=T
      XCOOR(INDEX)=r
      ZCOOR(INDEX)=Zz
      INDI(INDEX)=II
  800 CONTINUE
      rAX=r
      IF(MOUT.GE.1)WRITE(LOU,113)DD,r,XX,YY,Zz,T,PNEW,PAZM,
     1IND,IND1,ITER,II,INDEX
      IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)N,II
      IF(LU1.NE.0.AND.NDST.NE.0)
     1WRITE(LU1,108)(AY(2,I),AY(3,I),AY(4,I),I=1,N)
      IF(MDIM.EQ.0)GO TO 80
      CALL AMPL(APX,APY,APZ,UU)
      UU=1.
      SPR=1.
      CKMAH=(1.,0.)
      IF(MDIM.EQ.2)THEN
        IF(KMAH.NE.0)THEN
          PH=-1.57079632*KMAH
          CSKMAH=COS(PH)
          SNKMAH=SIN(PH)
          CKMAH=CMPLX(CSKMAH,SNKMAH)
        END IF
        SPR1=XDYN(2,1)*XDYN(3,2)-XDYN(3,1)*XDYN(2,2)
        SPR2=XDYN(3,1)*XDYN(1,2)-XDYN(1,1)*XDYN(3,2)
        SPR3=XDYN(1,1)*XDYN(2,2)-XDYN(2,1)*XDYN(1,2)
        SPRv=XDYN(1,3)*SPR1+XDYN(2,3)*SPR2+XDYN(3,3)*SPR3
        IF(NDST.EQ.0.and.mori.eq.0)CSF=COS(AA)
        IF(NDST.NE.0.and.mori.eq.0)CSF=COS(PNEW)
        if(mori.ne.0)csf=cos(pazm)
        if(mori.eq.0)sprv=-sprv
        VV=ay(5,n)*ay(5,n)+ay(6,n)*ay(6,n)+ay(7,n)*ay(7,n)
        SPR=SPRv*SQRT(VV)
        SPR=SQRT(ABS(SPR/csf))
        IF(MOUT.GE.2)WRITE(LOU,110)XDYN
        IF(MOUT.GE.2)WRITE(LOU,112)yDYN
      END IF
      DO 802 I=1,2
      APX(I)=APX(I)*UU*CKMAH/SPR
      APY(I)=APY(I)*UU*CKMAH/SPR
      APZ(I)=APZ(I)*UU*CKMAH/SPR
  802 CONTINUE
      IF(MOUT.GE.1.AND.MDIM.EQ.2)
     1WRITE(LOU,'(2X,F8.5,6(2X,E11.5)/10X,6(2X,E11.5),F10.5,I5)')
     2UU,(APX(I),APY(I),APZ(I),I=1,2),SPR,KMAH
      TAUX=T
      TAST=0.
      NCC=code(1,2)
      ncod=ncd
      IF(iani(isour).eq.0.and.ncc.eq.1)NCOD=-NCD
      call polar(1,1,1,1)
      IF(LU2.NE.0.AND.NDST.NE.0)then
        WRITE(LU2,116)ncod,II,T,APX(1),APY(1),APZ(1),TAST
        if(ncc.eq.1.and.ncod.lt.0)WRITE(LU2,115)APX(2),APY(2),APZ(2)
        WRITE(LU2,114)ay(5,1),ay(6,1),ay(7,1)
        if(ncc.eq.1)WRITE(LU2,114)(hhh(1,i),i=1,3)
        if(ncc.eq.1.and.ncod.lt.0)WRITE(LU2,114)(hhh(2,i),i=1,3)
        if(ncc.eq.2)WRITE(LU2,114)(hhh(2,i),i=1,3)
        if(ncc.eq.3)WRITE(LU2,114)(hhh(3,i),i=1,3)
      end if
      IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 801
      AMPX1(INDEX)=APX(1)
      AMPY1(INDEX)=APY(1)
      AMPZ1(INDEX)=APZ(1)
      if(ncc.eq.1.and.ncod.lt.0)then
        AMPX2(INDEX)=APX(2)
        AMPY2(INDEX)=APY(2)
        AMPZ2(INDEX)=APZ(2)
      end if
      p(index,1)=ay(5,1)
      p(index,2)=ay(6,1)
      p(index,3)=ay(7,1)
      if(ncc.eq.1)then
        pol(index,1)=hhh(1,1)
        pol(index,2)=hhh(1,2)
        pol(index,3)=hhh(1,3)
      end if
      if(ncc.eq.1.and.ncod.lt.0)then
        pol1(index,1)=hhh(2,1)
        pol1(index,2)=hhh(2,2)
        pol1(index,3)=hhh(2,3)
      end if
      if(ncc.eq.2)then
        pol(index,1)=hhh(2,1)
        pol(index,2)=hhh(2,2)
        pol(index,3)=hhh(2,3)
      end if
      if(ncc.eq.3)then
        pol(index,1)=hhh(3,1)
        pol(index,2)=hhh(3,2)
        pol(index,3)=hhh(3,3)
      end if
  801 CONTINUE
C
C
   80 IF(NDST.EQ.0.AND.LNDST.EQ.1)THEN
        BMIN=BMIN+BSTEP
        GO TO 4
      END IF
      d1=PNEW
      a1=pazm
      dr1=drs
      x1=xx
      y1=yy
      z1=zz
      T1=TAUX
      d2=ddNEW
      a2=aanew
      x2=xxnew
      y2=yynew
      z2=zznew
      T2=TTNEW
      IF(ITREND.EQ.(-1))GO TO 85
      II=II+1
      IF(II.GT.NDST.AND.IK1.NE.0)GO TO 6
      IF(method.eq.1)then
        aa=pnew
        inum=0
      end if
      IF(II.GT.NDST)GO TO 3
      DD=DST(II)
      if(method.eq.1)go to 60
      dr=rrNEW
      IF(IBNEW.EQ.1)dr=dr+REPS
      IF(DD.GT.dr.AND.IK1.NE.0)GO TO 6
      IF(DD.GT.dr)GO TO 3
      GO TO 60
   85 II=II-1
      IF(II.LT.1.AND.IK1.NE.0)GO TO 6
      IF(method.eq.1)then
        aa=pnew
        inum=0
      end if
      IF(II.LT.1)GO TO 3
      DD=DST(II)
      if(method.eq.1)go to 60
      dr=rrNEW
      IF(IBNEW.EQ.1)dr=dr-REPS
      IF(DD.LT.dr.AND.IK1.NE.0)GO TO 6
      IF(DD.LT.dr)GO TO 3
      GO TO 60
C
C
    6 CONTINUE
      IF(IK1.EQ.1)GO TO 7
      IK1=1
      d1=dNEW
      a1=anew
      d2=ddNEW
      a2=aanew
      IOLD1=INEW1
      indr1=iold1
      GO TO 59
    7 IK1=0
      rrOLD=rrNEW
      xxold=xxnew
      yyold=yynew
      zzold=zznew
      ddOLD=ddNEW
      aaold=aanew
      TTOLD=TTNEW
      IBOLD=IBNEW
      rrNEW=rNEW
      xxnew=xnew
      yynew=ynew
      zznew=znew
      ddNEW=dNEW
      aanew=anew
      TTNEW=TNEW
      IBNEW=0
      GO TO 41
C
C     E N D   O F    I T E R A T I O N S
C
C
C      BOUNDARY RAYS. CASE IOLD.NE.3, INEW=3
C
   50 rrNEW=rNEW
      xxnew=xnew
      yynew=ynew
      zznew=znew
      TTNEW=TNEW
      ddNEW=dNEW
      aanew=anew
      IBNEW=0
      d1=dOLD
      d2=dNEW
      a1=dOLD
      a2=dNEW
   54 IRES=0
      ITER=0
   51 PNEW=0.5*(d1+d2)
      ITER=ITER+1
      IND=INDS
      CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T,
     1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1)
      IF(IND.EQ.20)RETURN
      IF(ITPR.EQ.43)r=Zz
      IF(MOUT.GE.2)WRITE(LOU,103)IND,IND1,ITER,r,T,PNEW
      IF(IND.EQ.ITPR.AND.LNDST.EQ.1)GO TO 52
      d1=PNEW
      a1=pazm
      GO TO 53
   52 rrOLD=r
      xxold=xx
      yyold=yy
      zzold=zz
      ddOLD=PNEW
      aaold=pazm
      TTOLD=T
      IBOLD=1
      IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      rAX=r
      ZAX=Zz
      IAX=IND
      IAX1=IND1
      T1=T
      d2=PNEW
      a2=pazm
   53 IF(ITER.LT.ITMAX)GO TO 51
      IF(MOUT.GE.1.AND.IRES.EQ.1)
     1WRITE(LOU,107)rAX,ZAX,T1,d2,IAX,IAX1,IRES
      IF(MOUT.GE.1.AND.IRES.EQ.0)
     1WRITE(LOU,107)r,Zz,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.1)GO TO 41
      GO TO 3
C
C     BOUNDARY RAYS. CASE IOLD=3, INEW.NE.3
C     OR CASE  IOLD=3, INEW=3  BUT  IOLD1.NE.INEW1
C
   55 rrOLD=rOLD
      xxold=xold
      yyold=yold
      zzold=zold
      ddOLD=dOLD
      aaold=aold
      TTOLD=TOLD
      IBOLD=0
      d1=dOLD
      d2=dNEW
      a1=aOLD
      a2=aNEW
   59 IRES=0
      ITER=0
   56 PNEW=0.5*(d1+d2)
      ITER=ITER+1
      IND=INDS
      CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T,
     1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1)
      IF(IND.EQ.20)RETURN
      IF(ITPR.EQ.43)r=Zz
      IF(MOUT.GE.2)WRITE(LOU,103)IND,IND1,ITER,r,T,PNEW
      IF(IND.EQ.ITPR.AND.IBOLD.EQ.1.AND.LNDST.EQ.1)GO TO 57
      IF(IND.EQ.ITPR.AND.IND1.EQ.IOLD1.AND.LNDST.EQ.1)GO TO 57
      d2=PNEW
      a2=pazm
      GO TO 58
   57 rrNEW=r
      xxnew=xx
      yynew=yy
      zznew=zz
      ddNEW=PNEW
      aanew=pazm
      TTNEW=T
      IBNEW=1
      IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      rAX=r
      ZAX=Zz
      IAX=IND
      IAX1=IND1
      T2=T
      d1=PNEW
      a1=pazm
   58 IF(ITER.LT.ITMAX)GO TO 56
      IF(MOUT.GE.1.AND.IRES.EQ.1)
     1WRITE(LOU,107)rAX,ZAX,T2,d1,IAX,IAX1,IRES
      IF(MOUT.GE.1.AND.IRES.EQ.0)
     1WRITE(LOU,107)r,Zz,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.1.AND.IK1.EQ.1)GO TO 7
      IF(IRES.EQ.1)GO TO 41
      GO TO 3
C
C     CRITICAL RAYS. CASE IOLD.NE.9, IOLD.NE.3, INEW=9
C
   30 ITER=0
      XCR=rNEW
      d1=dOLD
      d2=dNEW
      IRES=0
   31 PNEW=0.5*(d1+d2)
      ITER=ITER+1
      IND=INDS
      CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T,
     1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1)
      IF(IND.EQ.20)RETURN
      IF(ITPR.EQ.43)r=Zz
      IF(MOUT.GE.2)WRITE(LOU,104)IND,IND1,ITER,r,T,PNEW
      IF(IND.EQ.9)GO TO 32
      IF(IND.EQ.ITPR)GO TO 33
      d1=PNEW
      GO TO 34
   32 IF(IND1.NE.INEW1)d1=PNEW
      IF(IND1.NE.INEW1)GO TO 34
      d2=PNEW
      IF(ITER.EQ.ITMAX.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 89
      XCR=r
      GO TO 34
   89 DO 86 K=1,KREF
      DO 86 L=1,2
   86 JC(K,L)=CODE(K,L)
      IRF3=IREF+3
      DO 87 K=IRF3,KREF
      DO 87 L=1,2
   87 CODE(K-2,L)=JC(K,L)
      KREF1=KREF
      KREF=KREF-2
      ICR=1
      ITER=ITMAX-1
      GO TO 31
   33 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      rAX=r
      rrnew=r
      xxnew=xx
      yynew=yy
      zznew=zz
      ZAX=Zz
      IAX=IND
      IAX1=IND1
      T2=T
      d1=PNEW
      PAP=PNEW
   34 IF(ITER.LT.ITMAX)GO TO 31
      IF(MOUT.GE.1.AND.IRES.EQ.1)
     1WRITE(LOU,111)rAX,ZAX,T2,PAP,IAX,IAX1,IRES
      IF(MOUT.GE.1.AND.IRES.EQ.0)
     1WRITE(LOU,111)r,Zz,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.0)GO TO 3
      d2=PAP
      rrNEW=rAX
      ddNEW=d2
      aanew=anew
      TTNEW=T2
      IBNEW=1
      d1=dOLD
      IF(ICR.EQ.0)GO TO 54
      ICR=0
      KREF=KREF1
      DO 88 K=1,KREF
      DO 88 L=1,2
   88 CODE(K,L)=JC(K,L)
      GO TO 54
C
C     CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.3.
C
   35 ITER=0
      XCR=rOLD
      d1=dOLD
      d2=dNEW
      IRES=0
   36 PNEW=0.5*(d1+d2)
      ITER=ITER+1
      IND=INDS
      CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T,
     1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1)
      IF(IND.EQ.20)RETURN
      IF(ITPR.EQ.43)r=Zz
      IF(MOUT.GE.2)WRITE(LOU,104)IND,IND1,ITER,r,T,PNEW
      IF(IND.EQ.9)GO TO 37
      IF(IND.EQ.ITPR)GO TO 38
      d2=PNEW
      GO TO 39
   37 IF(IND1.NE.INEW1)d2=PNEW
      IF(IND1.NE.INEW1)GO TO 39
      d1=PNEW
      IF(ITER.EQ.ITMAX.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 94
      XCR=r
      GO TO 39
   94 DO 91 K=1,KREF
      DO 91 L=1,2
   91 JC(K,L)=CODE(K,L)
      IRF3=IREF+3
      DO 92 K=IRF3,KREF
      DO 92 L=1,2
   92 CODE(K-2,L)=JC(K,L)
      KREF1=KREF
      KREF=KREF-2
      ICR=1
      ITER=ITMAX-1
      GO TO 36
   38 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      rAX=r
      rrold=r
      xxold=xx
      yyold=yy
      zzold=zz
      ZAX=Zz
      IAX=IND
      IAX1=IND1
      d2=PNEW
      PAP=PNEW
      T1=T
   39 IF(ITER.LT.ITMAX)GO TO 36
      IF(MOUT.GE.1.AND.IRES.EQ.1)
     1WRITE(LOU,111)rAX,ZAX,T1,PAP,IAX,IAX1,IRES
      IF(MOUT.GE.1.AND.IRES.EQ.0)WRITE(LOU,111)r,Zz,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.0)GO TO 3
      d1=PAP
      rrOLD=rAX
      ddOLD=d1
      aaold=aold
      TTOLD=T1
      IBOLD=1
      d2=dNEW
      IF(ICR.EQ.0)GO TO 59
      ICR=0
      KREF=KREF1
      DO 93 K=1,KREF
      DO 93 L=1,2
   93 CODE(K,L)=JC(K,L)
      GO TO 59
C
C
  100 FORMAT(3I3,2F10.5,F15.10)
  101 FORMAT(1X,'ITERATIONS',5X,4I3,3F10.5,2F15.10)
  102 FORMAT(5X,6I3,3F10.5,F15.10)
  103 FORMAT(2X,'BOUNDARY RAY',5X,3I3,2F10.5,F15.10)
  104 FORMAT(2X,'CRITICAL RAY',5X,3I3,2F10.5,F15.10)
  107 FORMAT(10X,3F10.5,F15.10,3I3,5X,'BOUNDARY RAY')
  108 FORMAT(3E15.5)
  109 FORMAT(I5,9E15.5)
  110 FORMAT(1X,'XDYN',3F10.5)
  111 FORMAT(10X,3F10.5,F15.10,3I3,5X,'CRITICAL RAY')
  112 FORMAT(1X,'YDYN',3F10.5)
  113 FORMAT(7F10.5,F15.10,5I3)
  114 FORMAT(3F10.5)
  115 format(6e12.6)
  116 FORMAT(2I3,F10.5,12E12.6,F10.6)
  117 FORMAT(6E15.5)
C
C
   99 N=0
      NAUX=0
      IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)N,NAUX
      IF(ncc.eq.1.and.ncod.lt.0)INDEX=-INDEX
      IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)INDEX
      IF(INDEX.EQ.0)RETURN
      IF(INDEX.LT.0)INDEX=-INDEX
      IF(LU1.NE.0.AND.NDST.NE.0)then
        do 200 i=1,index
        WRITE(LU1,109)INDI(I),XCOOR(I),ZCOOR(I),TIME(I),
     1  AMPX1(I),AMPY1(I),AMPZ1(I)
        if(ncc.eq.1.and.ncod.lt.0)
     1  write(lu1,117)AMPX2(I),AMPY2(I),AMPZ2(I)
        write(lu1,108)(p(i,k),k=1,3)
        write(lu1,108)(pol(i,k),k=1,3)
        if(ncc.eq.1.and.ncod.lt.0)
     1  write(lu1,108)(pol1(i,k),k=1,3)
  200   continue
      end if
      RETURN
      END
