SUBROUTINE PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PAZM,RANG,
     1XXX,YYY,ZZZ,TTT,DT,AC,ASTART,ASTEP,AFIN,ITMAX,MOUT,NCODE,
     2METHOD,ITPR,indr1)
C
C     3-D INITIAL VALUE RAY TRACING AND RAY TRACING FROM THE SOURCE
C     TO A PRESCRIBED PROFILE PASSING THROUGH THE EPICENTER
C
      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/ XDST(200),NDSTX,XREPS,DST(2),NDST,REPS,LNDST,
     1xprf,yprf
      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 /ZERO/ RNULL
      COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP
      COMMON/DYN/XDYN(3,3),ydyn(3,3)
      COMMON/KM/KMAH,SPROLD,TSTOLD
C
      iwave=0
      itp=code(1,2)
      RANG=0.
      XXX=0.
      ZZZ=0.
      TTT=0.
      PI=3.14159265
      REPS1=.05*REPS
      DD=dst(1)
      xcos=cos(dd)
      xsin=sin(dd)
      dd=0.
      X=0.
      ITER=0
      II=0
      LNDST=0
C
      AA=ASTART-ASTEP
      INDEX=0
      IPH=1
      INUM=0
      ICLS=0
      ISUC=0
      INDS=ISOUR
C
C     LOOP IN AZIMUTH, 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
      IND=INDS
      NDER=1
      IF(MDIM.GE.1)NDER=2
      NOLD=0
      SPROLD=0.
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,AA,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      NDER=1
      IF(IND.EQ.14.OR.IND.EQ.20)RETURN
      x=(yprf-yy)*xcos-(xprf-xx)*xsin
      IF(NDSTX.EQ.0)GO TO 65
      IF(IND.EQ.ITPR)XAX=X
      IF(IND.EQ.ITPR)PNEW=AA
      IF(MOUT.EQ.3)WRITE(LOU,100)IND,IND1,X,XX,YY,ZZ,T,AA,PSI0
      IF(INUM.NE.0)GO TO 2
      AOLD=AA
      IOLD=IND
      IOLD1=IND1
      XOLD=X
      TOLD=T
      INUM=1
      GO TO 1
C
C     PARAMETERS FOR THE PRECEDING RAY: AA=AOLD, X=XOLD, IND=IOLD
C     PARAMETERS FOR THE NEW RAY: AA=ANEW, X=XNEW, IND=INEW
C
    2 INEW=IND
      INEW1=IND1
      ANEW=AA
      XNEW=X
      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
      GOTO 3
   21 IF(IOLD1.NE.INEW1)then
        if(inew1.eq.indr1)go to 50
        if(iold1.eq.indr1)go to 55
      else
        GO TO 40
      end if
C
C     NO ITERATIONS, TAKE A NEW RAY IN THE LOOP
C
    3 CONTINUE
      if(isuc.eq.0)ind=0
      IF(IOLD.NE.INEW)IND=0
      IOLD=INEW
      IOLD1=INEW1
      XOLD=XNEW
      AOLD=ANEW
      TOLD=TNEW
      GO TO 1
C
C     REGULAR CASE: IOLD=3, INEW=3
C
   40 XXNEW=XNEW
      XXOLD=XOLD
      AANEW=ANEW
      AAOLD=AOLD
      TTNEW=TNEW
      TTOLD=TOLD
      IBNEW=0
      IBOLD=0
   41 IF(XXNEW.GT.XXOLD.AND.DST(2).GT.DST(1))GO TO 46
      IF(XXNEW.LT.XXOLD.AND.DST(2).LT.DST(1))GO TO 46
C
C     REGULAR CASE: XXNEW.LE.XXOLD, ITREND=-1 (REVERSE BRANCH)
C
      ITREND=-1
      DX=XXOLD
      IF(IBOLD.EQ.1) DX=DX+REPS
      IF(DD.GE.DX) GO TO 3
      DX=XXNEW
      IF(IBNEW.EQ.1) DX=DX-REPS
      IF(DD.LT.DX) GOTO 3
      II=1
      GO TO 43
C
C     REGULAR CASE: XXNEW.GT.XXOLD, ITREND=1 (REGULAR BRANCH)
C
   46 ITREND=1
      DX=XXOLD
      IF(IBOLD.EQ.1) DX=DX-REPS
      IF(DD.LE.DX) GO TO 3
      DX=XXNEW
      IF(IBNEW.EQ.1) DX=DX+REPS
      IF(DD.GT.DX) GOTO 3
      II=1
   43 P1=AAOLD
      P2=AANEW
      X1=XXOLD
      X2=XXNEW
      T1=TTOLD
      T2=TTNEW
C
C     I T E R A T I O N S
C
      ITER=0
      ISIGN=1
      IPR1=0
      IPR2=0
      ISUC=0
      GO TO 61
   68 XAX=X
      PAX=PNEW
   61 ITER=ITER+1
      IF(ITER.GT.ITMAX)GO TO 80
      ISIGN=-ISIGN
      AAUX=0.5*(P1+P2)
      IF(METHOD.LE.1.AND.IND.EQ.ITPR.and.iter.gt.1)GO TO 62
      GO TO 69
   62 if(mori.eq.0)AUX=(XDYN(1,1)*xsin-XDYN(2,1)*xcos)
      if(mori.ne.0)AUX=(XDYN(1,2)*xsin-XDYN(2,2)*xcos)
      IF(ABS(AUX).LT..00001)GO TO 69
      AAUX=PNEW+(DD-X)/AUX
   69 PNEW=AAUX
   71 IND=INDS
      SPROLD=0.
      XOLD=0.
      NOLD=0
      IF(MDIM.GE.1)NDER=2
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      if(iwave.eq.1)code(1,2)=itp
      NDER=1
      IF(IND.EQ.20)RETURN
      XE=XX-Xprf
      YE=YY-Yprf
      rPRF=SQRT(XE*XE+YE*YE)
      x=(yprf-yy)*xcos-(xprf-xx)*xsin
      IF(MOUT.EQ.3)WRITE(LOU,101)
     1IND,IND1,ITER,KMAH,DD,X,XX,YY,T,PNEW,PSI0
C
C     TESTING WHETHER THE RAY OF A QS WAVE DOES NOT TERMINATE
C     OUTSIDE THE RANGE IN WHICH PREVIOUS RAYS TERMINATED; IF YES,
C     A RAY OF THE OTHER QS WAVE WITH THE SAME INITIAL PARAMETERS
C     IS CALCULATED
C
      if((x-x1)*(x-x2).gt.0..and.itype.eq.1)then
        iter=iter+1
        if(iter.gt.itmax)go to 80
        code(1,2)=2
        iwave=1
        go to 71
      end if
      if((x-x1)*(x-x2).gt.0..and.itype.eq.2)then
        iter=iter+1
        if(iter.gt.itmax)go to 80
        code(1,2)=1
        iwave=1
        go to 71
      end if
      IF(ICLS.NE.0)GO TO 70
      IF(IND.NE.ITPR)P2=PNEW
      IF(IND.NE.ITPR)GO TO 61
      IF(ABS(X-XAX).LT..000001)GO TO 67
      IF(ABS(X-DD).LT.REPS)GO TO 65
      IF(X1.LT.X2.AND.DD.GT.X)GO TO 63
      IF(X1.GT.X2.AND.DD.LT.X)GO TO 63
      IF(ABS(X-X1).LT..000001)GO TO 67
      P2=PNEW
      X2=X
      T2=T
      IPR2=1
      GO TO 68
   63 IF(ABS(X-X2).LT..000001)GO TO 67
      P1=PNEW
      X1=X
      T1=T
      IPR1=1
      GO TO 68
   67 IF(ABS(PNEW-PAX).GT..000001)ITER=ITMAX
      AX1=X1-DD
      AX2=X2-DD
      IF((IPR1*IPR2).EQ.0)ITER=ITMAX
      X=X1
      PNEW=P1
      IF(ABS(AX1).GT.ABS(AX2))PNEW=P2
      IF(ABS(AX1).GT.ABS(AX2))X=X2
      IF(ITER.EQ.ITMAX)GO TO 61
      ICLS=1
      GO TO 69
   70 ICLS=0
      GO TO 65
C
C     SUCCESSFUL ITERATIONS
C
   65 INDEX=INDEX+1
      isuc=1
      RANG=rPRF
      if(xx.lt.xprf)rang=-rang
      XXX=XX
      YYY=YY
      ZZZ=ZZ
      TTT=T
      PAZM=PNEW
      XAX=X
      IF(MOUT.EQ.3)WRITE(LOU,113)DD,X,XX,YY,ZZ,T,PNEW,PSI0,
     1IND,IND1,ITER,II,INDEX
      GO TO 98
C
   80 NCD=NCODE
      P1=PNEW
      X1=X
      T1=T
      IF(ITER.GT.ITMAX)P1=AAOLD
      IF(ITER.GT.ITMAX)X1=XXOLD
      IF(ITER.GT.ITMAX)T1=TTOLD
      P2=AANEW
      X2=XXNEW
      T2=TTNEW
      GO TO 3
C
C
    6 continue
      P1=ANEW
      P2=AANEW
      IOLD1=INEW1
      GO TO 59
C
C     E N D   O F    I T E R A T I O N S
C
C     BOUNDARY RAYS: CASE IOLD.NE.ITPR, INEW=ITPR
C     OR CASE  IOLD=ITPR, INEW=ITPR  BUT  IOLD1.NE.INEW1
C     (IOLD1.NE.INDR1, INEW1=INDR1)
C
   50 XXNEW=XNEW
      TTNEW=TNEW
      AANEW=ANEW
      IBNEW=0
      P1=AOLD
      P2=ANEW
   54 IRES=0
      ITER=0
   51 PNEW=0.5*(P1+P2)
      ITER=ITER+1
      IND=INDS
      NDER=1
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      IF(IND.EQ.20)RETURN
      x=(yprf-yy)*xcos-(xprf-xx)*xsin
      IF(MOUT.EQ.3)WRITE(LOU,103)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0
      IF(IND.EQ.ITPR.AND.IND1.EQ.Indr1)GO TO 52
      P1=PNEW
      if((x-dd)*(xnew-dd).gt.0.)iter=itmax
      GO TO 53
   52 XXOLD=X
      AAOLD=PNEW
      TTOLD=T
      IBOLD=1
      if((x-dd)*(xnew-dd).lt.0.)iter=itmax
      IF(ABS(X-XAX).LT.REPS1)ITER=ITMAX
      IRES=1
      XAX=X
      T1=T
      P2=PNEW
   53 IF(ITER.LT.ITMAX)GO TO 51
      IF(MOUT.EQ.3)WRITE(LOU,107)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.1) GOTO 41
      GO TO 3
C
C     BOUNDARY RAYS: CASE IOLD=ITPR, INEW.NE.ITPR
C     OR CASE  IOLD=ITPR, INEW=ITPR  BUT  IOLD1.NE.INEW1
C     (IOLD1=INDR1, INEW1.NE.INDR1)
C
   55 XXOLD=XOLD
      AAOLD=AOLD
      TTOLD=TOLD
      IBOLD=0
      P1=AOLD
      P2=ANEW
   59 IRES=0
      ITER=0
   56 PNEW=0.5*(P1+P2)
      ITER=ITER+1
      IND=INDS
      NDER=1
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      IF(IND.EQ.20)RETURN
      x=(yprf-yy)*xcos-(xprf-xx)*xsin
      IF(MOUT.EQ.3)WRITE(LOU,103)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0
      IF(IND.EQ.ITPR.AND.IND1.EQ.Indr1)GO TO 57
      P2=PNEW
      if((x-dd)*(xold-dd).gt.0.)iter=itmax
      GO TO 58
   57 XXNEW=X
      AANEW=PNEW
      TTNEW=T
      IBNEW=1
      if((x-dd)*(xold-dd).lt.0.)iter=itmax
      IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1) ITER=ITMAX
      IRES=1
      XAX=X
      T2=T
      P1=PNEW
   58 IF(ITER.LT.ITMAX)GO TO 56
      IF(MOUT.EQ.3)WRITE(LOU,107)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.1)GOTO 41
      GO TO 3
C
C     CRITICAL RAYS. CASE IOLD.NE.9, IOLD.NE.3, INEW=9
C
   30 ITER=0
      XCR=XNEW
      P1=AOLD
      P2=ANEW
      IRES=0
   31 PNEW=0.5*(P1+P2)
      ITER=ITER+1
      IND=INDS
      NDER=1
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      IF(IND.EQ.20)RETURN
      x=(yprf-yy)*xcos-(xprf-xx)*xsin
      IF(MOUT.EQ.3)WRITE(LOU,104)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0
      IF(IND.EQ.9)GO TO 32
      IF(IND.EQ.ITPR)GO TO 33
      P1=PNEW
      GO TO 34
   32 CONTINUE
C 32  IF(IND1.NE.INEW1)P1=PNEW
C     IF(IND1.NE.INEW1) GOTO 34
      P2=PNEW
      IF(ABS(X-XCR).LT.0.01.AND.KC.NE.0.AND.IRES.EQ.1) GOTO 89
      XCR=X
      GOTO 34
   89 ITER=ITMAX-1
      GO TO 31
   33 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      XAX=X
      T2=T
      P1=PNEW
      PAP=PNEW
   34 IF(ITER.LT.ITMAX)GO TO 31
      IF(MOUT.EQ.3)WRITE(LOU,111)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.0) GOTO 3
      P2=PAP
      XXNEW=XAX
      AANEW=P2
      TTNEW=T2
      IBNEW=1
      P1=AOLD
      GO TO 54
C
C     CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.3.
C
   35 ITER=0
      P1=AOLD
      P2=ANEW
      IRES=0
   36 PNEW=0.5*(P1+P2)
      ITER=ITER+1
      IND=INDS
      NDER=1
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      IF(IND.EQ.20)RETURN
      x=(yprf-yy)*xcos-(xprf-xx)*xsin
      IF(MOUT.EQ.3)WRITE(LOU,104)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0
      IF(IND.EQ.9)GO TO 37
      IF(IND.EQ.ITPR)GO TO 38
      P2=PNEW
      GO TO 39
   37 IF(IND1.NE.IOLD1)P2=PNEW
      IF(IND1.NE.IOLD1) GO TO 39
      P1=PNEW
      IF(ABS(X-XCR).LT.0.01.AND.KC.NE.0.AND.IRES.EQ.1) GO TO 94
      XCR=X
      GO TO 39
   94 ITER=ITMAX-1
      GOTO 36
   38 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1) ITER=ITMAX
      IRES=1
      XAX=X
      P2=PNEW
      PAP=PNEW
      T1=T
   39 IF(ITER.LT.ITMAX)GO TO 36
      IF(MOUT.EQ.3)WRITE(LOU,111)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES
      IF(IRES.EQ.0) GOTO 3
      P1=PAP
      XXOLD=XAX
      AAOLD=P1
      TTOLD=T1
      IBOLD=1
      P2=ANEW
      GO TO 59
C
C
  100 FORMAT('*',2I3,5F10.5,2F15.10)
  101 FORMAT(1X,'*','ITERATIONS',5X,4I3,5F10.5,2F15.10)
  103 FORMAT(2X,'*','BOUNDARY RAY',5X,3I3,4F10.5,2F15.10)
  104 FORMAT(2X,'*','CRITICAL RAY',5X,3I3,4F10.5,2F15.10)
  107 FORMAT(10X,'*',5F10.5,F15.10,3I3,5X,'BOUNDARY RAY')
  111 FORMAT(10X,'*',5F10.5,F15.10,3I3,5X,'CRITICAL RAY')
  113 FORMAT('*',7F10.5,F15.10,5I3)
C
C
   98 CONTINUE
      LNDST=1
   99 CONTINUE
      RETURN
      END