C 
C     PROGRAM SEIS88
C     **************
C
C     PROGRAM SEIS88 IS DESIGNED FOR A TWO POINT RAY TRACING AND
C     RAY SYNTHETIC SEISMOGRAM COMPUTATION IN A 2-D LATERALLY
C     INHOMOGENEOUS MEDIA WITH CURVED INTERFACES AND BLOCK STRUCTURES.
C     COMPLETE WAVE FIELD (ALL THREE COMPONENTS) GENERATED BY
C     SIMPLE 3-D POINT SOURCES IS COMPUTED.
C
C     ****************************************************************
C
C
      INTEGER CODE
      DIMENSION MTEXT(17),VEL(6),Y(4)
      COMMON/DIST/DST(100),NDST,IDIST,ASTART,ASTEP,AFIN,REPS,REPS1,RSTEP
      COMMON/COD/CODE(50),KREF,KC
      COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20),
     1BRD(2),III(30,20),NPNT(20),NINT
      COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1
      COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT,
     1NTR,IMET
      COMMON/DEN/RHO1(19),RHO2(19),NRHO
      COMMON/SOUR/ROS,VPS,VSS
      COMMON/SH/KSH,NSH
      common/vsp/xvsp,icod,ivsp
C
C     ***
      mode=0
      call serv(mode,lin,lou,llu1,llu2,lu3)
      if(mode.eq.0)lin=5
      if(mode.eq.0)lou=6
C     ***
C
    1 READ(LIN,110)ldum1,ldum2,MPRINT,MTEXT,ldum,MM
      WRITE(LOU,110)ldum1,ldum2,MPRINT,MTEXT,ldum,MM
C
C
C     SPECIFICATION OF THE MODEL
C
      CALL MODEL(lin,lou,MM,MTEXT)
      IF(MM.LT.0)IND=19
      IF(MM.LT.0)WRITE(LOU,112)IND
      IF(MM.LT.0)GO TO 99
C
C     SPECIFICATION OF THE SYNTHETIC SEISMOGRAMS
C
    2 READ(LIN,100)ICONT,MEP,MOUT,MDIM,METHOD,IBP,IBS,iup,ius,IDP,IDS,
     1IREAD,MREG,ITMAX,NLAY,LU1,LU2,LTAPE,KSH,Iloc
      WRITE(LOU,102)ICONT,MEP,MOUT,MDIM,METHOD,IBP,IBS,iup,ius,IDP,IDS,
     1IREAD,MREG,ITMAX,NLAY,LU1,LU2,LTAPE,KSH,Iloc
      if(mode.eq.1)lu1=llu1
      if(mode.eq.1)lu2=llu2
      IF(ICONT.EQ.(-1))GO TO 1
      IF(ICONT.EQ.0)GO TO 99
C
C     ***
C     ***
C
      IF(LU1.NE.0)REWIND LU1
      IF(LU2.NE.0)REWIND LU2
      IF(NRHO.EQ.1)MREG=1
      IF(ITMAX.EQ.0)ITMAX=20
      IREC=1
      IMET=0
      ivsp=0
      itpr=iloc
      IF(ITPR.EQ.0)ITPR=3
      IF(ITPR.GT.100)IREC=ITPR-100
      if(iloc.eq.1)ivsp=1
      if(iloc.eq.1)itpr=22
      if(ivsp.eq.1)mreg=1
      if(ivsp.eq.1)irec=0
      IF(IREC.NE.1)MREG=1
      NLAYER=NINT-NLAY
      NDR=8
      NDER=0
      IF(MDIM.NE.0)NDER=1
      IF(MDIM.NE.0)NDR=11
C
      IF(MEP.GT.0)GO TO 3
      NDST=-MEP
      READ(LIN,104)(DST(I),I=1,NDST),xvsp
      WRITE(LOU,104)(DST(I),I=1,NDST),xvsp
      if(ndst.eq.1)rstep=1.
      if(ndst.eq.1)dst(2)=dst(1)+1.
      if(ndst.eq.1)go to 4
      RSTEP=(DST(NDST)-DST(1))/FLOAT(NDST-1)
      GO TO 4
    3 NDST=MEP
      READ(LIN,104)RMIN,RSTEP,xvsp
      WRITE(LOU,104)RMIN,RSTEP,xvsp
      DO 13 I=1,MEP
   13 DST(I)=RMIN+(I-1)*RSTEP
C
    4 READ(LIN,104)XSOUR,ZSOUR,TSOUR,REPS,reps1,BRD(1),BRD(2)
      WRITE(LOU,104)XSOUR,ZSOUR,TSOUR,REPS,reps1,BRD(1),BRD(2)
      READ(LIN,104)DT,AMIN1,ASTEP1,AMAX1,AMIN2,ASTEP2,AMAX2,AC
      WRITE(LOU,104)DT,AMIN1,ASTEP1,AMAX1,AMIN2,ASTEP2,AMAX2,AC
      IF(DT.LT.0.0000001)DT=1.
      IF(AC.LT.0.0000001)AC=0.0001
      TMAX=10000.
      IF(REPS.LT..00001)REPS=.05
      IF(REPS1.LT..00001)REPS1=.05*REPS
      IND=-1
      NC=NPNT(1)
      BL=X1(1,1)
      BR=X1(NC,1)
      BA=BRD(1)
      BB=BRD(2)
      IF(ABS(BB-BA).GT..00001)GO TO 5
      BRD(1)=BL
      BRD(2)=BR
    5 IF(BA.LT.BL)BRD(1)=BL
      IF(BB.GT.BR)BRD(2)=BR
c
c     determination of the position of the source in the model
c
      if(xsour.lt.brd(1).or.xsour.gt.brd(2))go to 205
      intr=1
  201 nl=npnt(intr)-1
      do 202 i=1,nl
      j=i
      if(xsour.lt.x1(i+1,intr))go to 203
  202 continue
  203 a2=a1(j,intr)
      b2=b1(j,intr)
      c2=c1(j,intr)
      d2=d1(j,intr)
      x2=x1(j,intr)
      au=xsour-x2
      zint=((d2*au+c2)*au+b2)*au+a2
      if(abs(zsour).gt..0001)go to 204
      isour=1
      zsour=zint
      go to 206
  204 if(zsour.lt.zint.and.intr.eq.1)go to 205
      if(zsour.lt.zint)go to 206
      if(abs(zsour-zint).lt..000001.and.intr.eq.nint)go to 206
      if(intr.eq.nint)go to 205
      isour=intr
      intr=intr+1
      go to 201
  205 ind=50
      write(6,112)ind
      go to 99
  206 lay=isour
      int1=isour
      ITYPE=-1
      Y(1)=XSOUR
      Y(2)=ZSOUR
      Y(3)=AMIN1
      Y(4)=0.
      CALL VELOC(Y,VEL)
      VSS=VEL(1)
      ITYPE=1
      CALL VELOC(Y,VEL)
      VPS=VEL(1)
      ROS=RHO1(Isour)+RHO2(Isour)*VPS
C
C
C     GENERATE FILE LU2 FOR PLOTTING OF SYNTHETIC SEISMOGRAMS
C
      IF(LU2.EQ.0)GO TO 25
      IF(LTAPE.EQ.0)WRITE(LU2,115)MTEXT
      IF(LTAPE.EQ.1)WRITE(LU2)MTEXT
      IF(LTAPE.EQ.0)WRITE(LU2,100)NDST,KSH,itpr
      IF(LTAPE.EQ.1)WRITE(LU2)NDST,KSH,itpr
      IF(LTAPE.EQ.0)WRITE(LU2,104)XSOUR,ZSOUR,TSOUR,RSTEP,ROS,VPS,VSS
      IF(LTAPE.EQ.1)WRITE(LU2)XSOUR,ZSOUR,TSOUR,RSTEP,ROS,VPS,VSS
      IF(LTAPE.EQ.0)WRITE(LU2,104)(DST(I),I=1,NDST)
      IF(LTAPE.EQ.1)WRITE(LU2)(DST(I),I=1,NDST)
   25 CONTINUE
C
C    LOOP FOR ELEMENTARY WAVES
C
      NAWX=0
      IF(MOUT.EQ.0)WRITE(LOU,105)
   20 CALL GENER(NAWX,IDP,IDS,IBP,IBS,IUP,IUS,IREAD,ISOUR,IREC,NLAYER,
     1NLAY,IFIN,NCODE,lin,lou)
      NAWX=1
      IF(IFIN.EQ.1)GO TO 49
      IF(MOUT.NE.0)WRITE(LOU,107)
      WRITE(LOU,106)NCODE,KC,KREF,(CODE(I),I=1,KREF)
      icod=0
      if(ivsp.eq.0)go to 35
      do 34 i=1,kref
      icod=kref-i+1
      if(icod.eq.1)go to 34
      ic1=code(icod)
      ic2=code(icod-1)
      if((ic1-ic2).eq.0)go to 35
      if((ic1*ic2).lt.0)go to 35
   34 continue
   35 continue
      NSH=0
      DO 301 I=1,KREF
      IF(CODE(I).GT.0)GO TO 302
  301 CONTINUE
      NSH=KSH
  302 CONTINUE
      IF(KSH.EQ.1.AND.NSH.EQ.0)GO TO 20
      IF(MOUT.NE.0)WRITE(LOU,108)
      IF(IFIN.EQ.(-1))GO TO 20
C
C     GENERATE FILE LU1 FOR PLOTTING OF RAY DIAGRAMS,TRAVEL TIMES AND
C     AMPLITUDES
C
      IF(LU1.EQ.0)GO TO 21
      if(ltape.eq.0)WRITE(LU1,100)ICONT,itpr
      if(ltape.eq.1)WRITE(LU1)ICONT,itpr
      if(ltape.eq.0)WRITE(LU1,100)NINT,(NPNT(I),I=1,NINT)
      if(ltape.eq.1)WRITE(LU1)NINT,(NPNT(I),I=1,NINT)
      DO 22 IIJ=1,NINT
      NI=NPNT(IIJ)-1
      if(ltape.eq.0)
     1WRITE(LU1,101)(A1(J,IIJ),B1(J,IIJ),C1(J,IIJ),D1(J,IIJ),X1(J,IIJ),
     2III(J,IIJ),J=1,NI)
      if(ltape.eq.1)
     1WRITE(LU1)(A1(J,IIJ),B1(J,IIJ),C1(J,IIJ),D1(J,IIJ),X1(J,IIJ),
     2III(J,IIJ),J=1,NI)
   22 continue
      if(ltape.eq.0)WRITE(LU1,104)XSOUR,ZSOUR,ROS,VPS,VSS
      if(ltape.eq.1)WRITE(LU1)XSOUR,ZSOUR,ROS,VPS,VSS
   21 CONTINUE
C
C
      ASTART=AMIN1
      ASTEP=ASTEP1
      AFIN=AMAX1
      IF(KC.LE.0)ASTART=AMIN2
      IF(KC.LE.0)ASTEP=ASTEP2
      IF(KC.LE.0)AFIN=AMAX2
      CALL TRAMP(XSOUR,ZSOUR,TSOUR,DT,AC,TMAX,ITMAX,MOUT,
     1MDIM,NCODE,lou,LU1,LU2,METHOD,ISOUR,LTAPE,ITPR)
      IF(IND.EQ.20.OR.IND.EQ.21)WRITE(LOU,111)IND,LAY
      IF(IND.EQ.20.OR.IND.EQ.21)GO TO 99
      IF(IND.EQ.14)WRITE(LOU,112)IND
      GO TO 20
C
C     END OF LOOP FOR ELEMENTARY WAVES
C
C
   49 IF(LU1.EQ.0)GO TO 50
      JJ=0
      if(ltape.eq.0)WRITE(LU1,100)JJ,JJ,JJ
      if(ltape.eq.1)WRITE(LU1)JJ,JJ,JJ
      REWIND LU1
   50 IF(LU2.EQ.0)GO TO 2
      REWIND LU2
C
C     ***
C     ***
C
      GO TO 2
C
  100 FORMAT(26I3)
  101 format(5e15.5,i5)
  102 FORMAT(1H1,////,2X,26I3)
  104 FORMAT(8F10.5)
  105 FORMAT(//2X,'LISTING OF WAVE CODES'//2X,'INT.CODE',5X,'E X T E R N
     1 A L   C O D E')
  106 FORMAT(I10,5X,31I3/18X,30I3)
  107 FORMAT(//2X,'INT.CODE',5X,'E X T E R N A L   C O D E')
  108 FORMAT(//)
  110 FORMAT(3I2,17A4,2X,2I2)
  111 FORMAT(/2X,'IND=',I5,2X,'LAY=',I5/)
  112 format(/2x,'IND=',i5/)
  115 FORMAT(17A4)
C
C     ***
   99 CONTINUE
C     ***
C
      STOP
      END
C
C     ***************************************************************
C
      SUBROUTINE TRAMP(XSOUR,ZSOUR,TSOUR,DT,AC,TMAX,ITMAX,MOUT,
     1MDIM,NCODE,lou,LU1,LU2,METHOD,ISOUR,LTAPE,ITPR)
C
C     TWO-POINT RAY TRACING
C
      INTEGER CODE
      DIMENSION TIME(400),XCOOR(400),AMPX(400),AMPY(400),AMPZ(400),
     1INDI(400),TAS(400),ANG(400),PHAX(400),PHAY(400),PHAZ(400),JC(50),
     2QP(15)
      COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1
      COMMON/DIST/DST(100),NDST,IDIST,ASTART,ASTEP,AFIN,REPS,REPS1,RSTEP
      COMMON/COD/CODE(50),KREF,KC
      COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),XX(30,20),
     1BRD(2),III(30,20),NPNT(20),NINT
      COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19),
     1NQP(19),NQS(19),NABS
      COMMON/SH/KSH,NSH
      common/vsp/xvsp,icod,ivsp
C

      TAST=0.
      PI=3.14159265
C
      AA=ASTART-ASTEP
      INDEX=0
      INUM=0
      IK1=0
      ICLS=0
      ICR=0
      ISUC=0
      INDS=ISOUR
C
C     LOOP IN RAY PARAMETERS AA, FROM ASTART TO AFIN, WITH THE STEP
C     ASTEP
C
    1 AA=AA+ASTEP
      IF(ASTEP.GT.0..AND.AA.GT.AFIN)GO TO 99
      IF(ASTEP.LT.0..AND.AA.LT.AFIN)GO TO 99
      IND=INDS
      CALL RAY2(XSOUR,ZSOUR,TSOUR,AA,X,Z,T,FI,TMAX,DT,AC)
      IF(IND.EQ.14.OR.IND.EQ.20.OR.IND.EQ.21)RETURN
      if(itpr.eq.22)x=z
      IF(IND.EQ.ITPR)XAX=X
      IF(IND.EQ.ITPR)PNEW=AA
      IF(IND.EQ.9)XCR=X
      IF(MOUT.EQ.2)WRITE(LOU,100)IND,IND1,X,T,AA
      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
      GO TO 3
   21 IF(IOLD1.NE.INEW1)IK1=2
      IF(IK1.EQ.2)GO TO 55
      GO TO 40
C
C     NO ITERATIONS, TAKE A NEW RAY IN THE LOOP
C
    3 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
      DO 42 I=1,NDST
      NDD=NDST-I+1
      DD=DST(NDD)
      DX=XXOLD
      IF(IBOLD.EQ.1)DX=DX+REPS
      IF(DD.GE.DX)GO TO 42
      DX=XXNEW
      IF(IBNEW.EQ.1)DX=DX-REPS
      IF(DD.LT.DX.AND.IK1.NE.0)GO TO 6
      IF(DD.LT.DX)GO TO 3
      II=NDD
      GO TO 43
   42 CONTINUE
      IF(IK1.NE.0)GO TO 6
      GO TO 3
C
C     REGULAR CASE: XXNEW.GT.XXOLD, ITREND=1 (REGULAR BRANCH)
C
   46 ITREND=1
      DO 44 I=1,NDST
      DD=DST(I)
      DX=XXOLD
      IF(IBOLD.EQ.1)DX=DX-REPS
      IF(DD.LE.DX)GO TO 44
      DX=XXNEW
      IF(IBNEW.EQ.1)DX=DX+REPS
      IF(DD.GT.DX.AND.IK1.NE.0)GO TO 6
      IF(DD.GT.DX)GO TO 3
      II=I
      GO TO 43
   44 CONTINUE
      IF(IK1.NE.0)GO TO 6
      GO TO 3
   43 IF(METHOD.NE.3)GO TO 45
   49 NPT=NPNT(1)-1
      DO 47 I=1,NPT
      J=I
      IF(DD.LT.XX(I+1,1))GO TO 48
   47 CONTINUE
   48 AU=DD-XX(J,1)
      DDZ=((D1(J,1)*AU+C1(J,1))*AU+B1(J,1))*AU+A1(J,1)
      IF(ISUC.EQ.1)GO TO 60
   45 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
   60 ITER=0
      IF(ABS(X-DD).LT.REPS)GO TO 65
      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
      IF(METHOD.EQ.3.AND.IND.EQ.itpr)GO TO 66
      ISIGN=-ISIGN
      PNEW=0.5*(P1+P2)
      IF(ISIGN.LT.0.OR.METHOD.GT.1)GO TO 69
      AUX=X2-X1
      IF(ABS(AUX).LT..000001)GO TO 69
      AUX=P1+(DD-X1)*(P2-P1)/AUX
      IF((AUX-P1)*(AUX-P2).GT.0.)GO TO 69
      PNEW=AUX
      GO TO 69
   66 if(itpr.eq.22)then
        dx=0.
        dz=dd-x
      end if
      if(itpr.ne.22)then
        DZ=DDZ-Z
        DX=DD-X
      end if
      NNN=N
      IREF1=IREF
   64 N=KINT(IREF1)
      IF(N.LT.0)IREF1=IREF1-1
      IF(N.LT.0)GO TO 64
      V=AY(6,N)
      TX=AY(4,N)*V
      TZ=AY(5,N)*V
      DN=DX*TZ-DZ*TX
      DSS=-(DX*TX+DZ*TZ)
      N=0
      KMAH=0
      CALL JPAR(qp,afact,KMAH)
      N=NNN
      QS=Qp(3)+V*qP(4)*DSS
      IF(ABS(QS).LT..0001)then
        pnew=0.5*(p1+p2)
        go to 69
      end if
      V0=AY(6,1)
      SN=(DN*V0)/QS
      PNEW=PNEW-SN
      IF(p1.lt.p2.and.(pnew.lt.p1.or.pnew.gt.p2))then
        pnew=0.5*(p1+p2)
        go to 69
      end if
      IF(p1.gt.p2.and.(pnew.gt.p1.or.pnew.lt.p2))then
        pnew=0.5*(p1+p2)
        go to 69
      end if
   69 IND=INDS
      CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC)
      IF(IND.EQ.20.OR.IND.EQ.21)RETURN
      if(itpr.eq.22)x=z
      IF(MOUT.EQ.2)WRITE(LOU,101)IND,IND1,ITER,DD,X,T,PNEW
      IF(ICLS.NE.0)GO TO 70
      IF(IND.NE.ITPR)P2=PNEW
      IF(IND.NE.ITPR)GO TO 61
      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
   62 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
      XAX=X
      IF(LU1.EQ.0)GO TO 800
      TIME(INDEX)=T
      XCOOR(INDEX)=X
      INDI(INDEX)=IND
      ANG(INDEX)=PNEW
  800 CONTINUE
      IF(ITPR.NE.22)AUX=AY(4,N)
      IF(ITPR.EQ.22)AUX=AY(5,N)
      TAUX=T
      T=Taux+AUX*(DD-X)
      IF(MOUT.EQ.2)WRITE(LOU,102)IND,IND1,ITER,II,INDEX,DD,X,T,PNEW
      IF(MOUT.EQ.2.AND.MDIM.EQ.0)WRITE(LOU,114)
      if(itpr.eq.22)z=x
      if(itpr.eq.22)x=xvsp
      IF(MOUT.EQ.1.AND.MDIM.EQ.0)WRITE(LOU,113)DD,X,Z,T,PNEW,IND,
     1IND1,ITER,II
      if(itpr.eq.22)x=z
      if(lu1.eq.0)go to 802
      if(ltape.eq.0)WRITE(LU1,108)N,IND
      if(ltape.eq.0)WRITE(LU1,109)(AY(2,I),AY(3,I),I=1,N)
      if(ltape.eq.1)WRITE(LU1)N,IND
      if(ltape.eq.1)WRITE(LU1)(AY(2,I),AY(3,I),I=1,N)
  802 continue
      PH=0.
      IF(MDIM.EQ.0)GO TO 80
      IF(MDIM.EQ.1)SPR=1.
      IF(MDIM.EQ.1)GO TO 10
      NNN=N
      KMAH=0
      N=0
      CALL JPAR(QP,AFACT,KMAH)
      PH=PH-1.57079632*KMAH
      SPR=ABS(QP(3)/AY(6,1))
      N=NNN
   10 CALL AMPL(SPR,AX,AAY,AZ,PHX,PHY,PHZ)
      PHX=PHX+PH
      PHY=PHY+PH
      PHZ=PHZ+PH
      TAST=QP(5)
      IF(MDIM.LT.3)GO TO 12
      AUX=QP(6)/AY(6,1)
      AUX=SQRT(1./ABS(AUX))
      SPR=SPR/AUX
      AX=AX*AUX
      AAY=AAY*AUX
      AZ=AZ*AUX
   12 CONTINUE
      IF(PHX.GE.(-PI))GO TO 13
      PHX=PHX+2.*PI
      GO TO 12
   13 IF(PHX.LE.PI)GO TO 14
      PHX=PHX-2.*PI
      GO TO 13
   14 IF(PHZ.GE.(-PI))GO TO 15
      PHZ=PHZ+2.*PI
      GO TO 14
   15 IF(PHZ.LE.PI)GO TO 16
      PHZ=PHZ-2.*PI
      GO TO 15
   16 CONTINUE
      IF(LU1.EQ.0)GO TO 801
      TAS(INDEX)=TAST
      PHAX(INDEX)=PHX
      PHAY(INDEX)=PHY
      PHAZ(INDEX)=PHZ
      AMPZ(INDEX)=AZ
      AMPY(INDEX)=AAY
      AMPX(INDEX)=AX
  801 CONTINUE
      if(itpr.eq.22)z=x
      if(itpr.eq.22)x=xvsp
      IF(MOUT.EQ.1)WRITE(LOU,105)DD,X,Z,T,TAST,PNEW,SPR,
     1IND,IND1,ITER,II,KMAH,AX,PHX,AAY,PHY,AZ,PHZ
      if(itpr.eq.22)x=z
      IF(MOUT.EQ.2)WRITE(LOU,110)AX,PHX,AAY,PHY,AZ,PHZ,SPR,TAST,KMAH
      NCD=NCODE
      IF(CODE(1).LT.0)NCD=-NCODE
C
      IF(LU2.NE.0.AND.LTAPE.EQ.0)
     1WRITE(LU2,116)NCD,II,T,AX,AAY,AZ,PHX,PHY,PHZ,PNEW,TAST
      IF(LU2.NE.0.AND.LTAPE.EQ.1)
     1WRITE(LU2)NCD,II,T,AX,AAY,AZ,PHX,PHY,PHZ,PNEW,TAST
C
C
C
   80 P1=PNEW
      X1=X
      T1=TAUX
      IF(ITER.GT.ITMAX)P1=AAOLD
      IF(ITER.GT.ITMAX)X1=XXOLD
      IF(ITER.GT.ITMAX)T1=TTOLD
      P2=AANEW
      X2=XXNEW
      T2=TTNEW
      IF(ITREND.EQ.(-1))GO TO 85
      II=II+1
      IF(II.GT.NDST.AND.IK1.NE.0)GO TO 6
      IF(II.GT.NDST)GO TO 3
      DD=DST(II)
      DX=XXNEW
      IF(IBNEW.EQ.1)DX=DX+REPS
      IF(DD.GT.DX.AND.IK1.NE.0)GO TO 6
      IF(DD.GT.DX)GO TO 3
      ISUC=1
      IF(METHOD.EQ.3)GO TO 49
      GO TO 60
   85 II=II-1
      IF(II.LT.1.AND.IK1.NE.0)GO TO 6
      IF(II.LT.1)GO TO 3
      DD=DST(II)
      DX=XXNEW
      IF(IBNEW.EQ.1)DX=DX-REPS
      IF(DD.LT.DX.AND.IK1.NE.0)GO TO 6
      IF(DD.LT.DX)GO TO 3
      ISUC=1
      IF(METHOD.EQ.3)GO TO 49
      GO TO 60
C
C
    6 IF(IK1.EQ.1)GO TO 7
      IK1=1
      P1=ANEW
      P2=AANEW
      IOLD1=INEW1
      GO TO 59
    7 IK1=0
      XXOLD=XXNEW
      AAOLD=AANEW
      TTOLD=TTNEW
      IBOLD=IBNEW
      XXNEW=XNEW
      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.ITPR, INEW=ITPR
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
      CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC)
      IF(IND.EQ.20.OR.IND.EQ.21)RETURN
      if(itpr.eq.22)x=z
      IF(MOUT.EQ.2)WRITE(LOU,103)IND,IND1,ITER,X,T,PNEW
      IF(IND.EQ.ITPR)GO TO 52
      P1=PNEW
      GO TO 53
   52 XXOLD=X
      AAOLD=PNEW
      TTOLD=T
      IBOLD=1
      IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      XAX=X
      ZAX=Z
      IAX=IND
      IAX1=IND1
      T1=T
      P2=PNEW
   53 IF(ITER.LT.ITMAX)GO TO 51
      IF(MOUT.EQ.1.AND.IRES.EQ.1)
     1WRITE(LOU,107)XAX,ZAX,T1,P2,IAX,IAX1,IRES
      if(itpr.eq.22)z=x
      if(itpr.eq.22)x=xvsp
      IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,107)X,Z,T,PNEW,IND,IND1,IRES
      if(itpr.eq.22)x=z
      IF(IRES.EQ.1)GO TO 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
   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
      CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC)
      IF(IND.EQ.20.OR.IND.EQ.21)RETURN
      if(itpr.eq.22)x=z
      IF(MOUT.EQ.2)WRITE(LOU,103)IND,IND1,ITER,X,T,PNEW
      IF(IND.EQ.ITPR.AND.IBOLD.EQ.1)GO TO 57
      IF(IND.EQ.ITPR.AND.IND1.EQ.IOLD1)GO TO 57
      P2=PNEW
      GO TO 58
   57 XXNEW=X
      AANEW=PNEW
      TTNEW=T
      IBNEW=1
      IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      XAX=X
      ZAX=Z
      IAX=IND
      IAX1=IND1
      T2=T
      P1=PNEW
   58 IF(ITER.LT.ITMAX)GO TO 56
      IF(MOUT.EQ.1.AND.IRES.EQ.1)
     1WRITE(LOU,107)XAX,ZAX,T2,P1,IAX,IAX1,IRES
      if(itpr.eq.22)z=x
      if(itpr.eq.22)x=xvsp
      IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,107)X,Z,T,PNEW,IND,IND1,IRES
      if(itpr.eq.22)x=z
      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.ITPR, 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
      CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC)
      IF(IND.EQ.20.OR.IND.EQ.21)RETURN
      if(itpr.eq.22)x=z
      IF(MOUT.EQ.2)WRITE(LOU,104)IND,IND1,ITER,X,T,PNEW
      IF(IND.EQ.9)GO TO 32
      IF(IND.EQ.ITPR)GO TO 33
      P1=PNEW
      GO TO 34
   32 IF(IND1.NE.INEW1)P1=PNEW
      IF(IND1.NE.INEW1)GO TO 34
      P2=PNEW
      IF(ABS(X-XCR).LT..01.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 89
      XCR=X
      GO TO 34
   89 DO 86 K=1,KREF
   86 JC(K)=CODE(K)
      IRF3=IREF+3
      DO 87 K=IRF3,KREF
   87 CODE(K-2)=JC(K)
      KREF1=KREF
      KREF=KREF-2
      ICR=1
      ITER=ITMAX-1
      GO TO 31
   33 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      XAX=X
      ZAX=Z
      IAX=IND
      IAX1=IND1
      T2=T
      P1=PNEW
      PAP=PNEW
   34 IF(ITER.LT.ITMAX)GO TO 31
      IF(MOUT.EQ.1.AND.IRES.EQ.1)
     1WRITE(LOU,111)XAX,ZAX,T2,PAP,IAX,IAX1,IRES
      if(itpr.eq.22)z=x
      if(itpr.eq.22)x=xvsp
      IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,111)X,Z,T,PNEW,IND,IND1,IRES
      if(itpr.eq.22)x=z
      IF(IRES.EQ.0)GO TO 3
      P2=PAP
      XXNEW=XAX
      AANEW=P2
      TTNEW=T2
      IBNEW=1
      P1=AOLD
      IF(ICR.EQ.0)GO TO 54
      ICR=0
      KREF=KREF1
      DO 88 K=1,KREF
   88 CODE(K)=JC(K)
      GO TO 54
C
C     CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.ITPR
C
   35 ITER=0
      XCR=0
      P1=AOLD
      P2=ANEW
      IRES=0
   36 PNEW=0.5*(P1+P2)
      ITER=ITER+1
      IND=INDS
      CALL RAY2(XSOUR,ZSOUR,TSOUR,PNEW,X,Z,T,FI,TMAX,DT,AC)
      IF(IND.EQ.20.OR.IND.EQ.21)RETURN
      if(itpr.eq.22)x=z
      IF(MOUT.EQ.2)WRITE(LOU,104)IND,IND1,ITER,X,T,PNEW
      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..01.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 94
      XCR=X
      GO TO 39
   94 DO 91 K=1,KREF
   91 JC(K)=CODE(K)
      IRF3=IREF+3
      DO 92 K=IRF3,KREF
   92 CODE(K-2)=JC(K)
      KREF1=KREF
      KREF=KREF-2
      ICR=1
      ITER=ITMAX-1
      GO TO 36
   38 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX
      IRES=1
      XAX=X
      ZAX=Z
      IAX=IND
      IAX1=IND1
      P2=PNEW
      PAP=PNEW
      T1=T
   39 IF(ITER.LT.ITMAX)GO TO 36
      IF(MOUT.EQ.1.AND.IRES.EQ.1)
     1WRITE(LOU,111)XAX,ZAX,T1,PAP,IAX,IAX1,IRES
      if(itpr.eq.22)z=x
      if(itpr.eq.22)x=xvsp
      IF(MOUT.EQ.1.AND.IRES.EQ.0)WRITE(LOU,111)X,Z,T,PNEW,IND,IND1,IRES
      if(itpr.eq.22)x=z
      IF(IRES.EQ.0)GO TO 3
      P1=PAP
      XXOLD=XAX
      AAOLD=P1
      TTOLD=T1
      IBOLD=1
      P2=ANEW
      IF(ICR.EQ.0)GO TO 59
      ICR=0
      KREF=KREF1
      DO 93 K=1,KREF
   93 CODE(K)=JC(K)
      GO TO 59
C
C
  100 FORMAT(2I3,2F10.5,F15.10)
  101 FORMAT(1X,'ITERATIONS',5X,3I3,3F10.5,F15.10)
  102 FORMAT(/,1X,'SUCCESSFUL ITERATION',3X,5I3,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)
  105 FORMAT(4F10.5,2F10.6,F10.3,5I3/13X,3(E15.5,F10.5))
  106 FORMAT(2I5,5E15.10)
  107 FORMAT(10X,3F10.5,F15.10,3I3,5X,'BOUNDARY RAY')
  108 format(2i5)
  109 format(2e15.5)
  110 FORMAT(5X,3(E15.5,F10.5),2F15.7,I5/)
  111 FORMAT(10X,3F10.5,F15.10,3I3,5X,'CRITICAL RAY')
  113 FORMAT(4F10.5,F15.10,4I3)
  114 FORMAT(/)
  115 FORMAT(2I3,F10.5,2E12.6,4F10.6)
  116 FORMAT(2I3,F10.5,3E12.6,5F10.6)
  117 format(i5,10e15.5)
C
C
   99 N=0
      NAUX=0
      if(lu1.eq.0)return
      if(ltape.eq.0)WRITE(LU1,108)N,NAUX
      if(ltape.eq.1)WRITE(LU1)N,NAUX
      IF(CODE(1).LT.0)INDEX=-INDEX
      IF(ltape.eq.0)WRITE(LU1,108)INDEX
      IF(ltape.eq.1)WRITE(LU1)INDEX
      IF(INDEX.EQ.0)RETURN
      IF(INDEX.LT.0)INDEX=-INDEX
      IF(ltape.eq.0)WRITE(LU1,117)(INDI(I),XCOOR(I),TIME(I),TAS(I),
     1ANG(I),AMPX(I),AMPY(I),AMPZ(I),PHAX(I),PHAY(I),PHAZ(I),I=1,INDEX)
      IF(ltape.eq.1)WRITE(LU1)(INDI(I),XCOOR(I),TIME(I),TAS(I),ANG(I),
     1AMPX(I),AMPY(I),AMPZ(I),PHAX(I),PHAY(I),PHAZ(I),I=1,INDEX)
      RETURN
      END
C
C     **************************************************************
C
      SUBROUTINE GENER(N,IDP,IDS,IBP,IBS,IUP,IUS,IREAD,IS,IR,NS,NL,IFIN,
     1NCODE,lin,lou)
C
C     THE ROUTINE GENER IS DESIGNED TO GENERATE NUMERICAL CODES OF
C     ELEMENTARY SEISMIC BODY WAVES
C
      DIMENSION JCA(50)
      COMMON/COD/JC(50),KCA,KC
C
C
      IFIN=0
      if(is.le.0.or.is.gt.ns)go to 101
      if(ir.lt.0.or.ir.gt.ns)go to 101
      irc=ir
      IF(N.GT.0)GO TO 1
      NCODE=1
      INEW=1
      INDEX=1
      GO TO 51
    1 NCODE=NCODE+1
      IF(INDEX.GE.11)GO TO 80
      IF(INEW.EQ.1)GO TO 50
      JAUX=0
      DO 20 I=1,KCA
      IF(JCA(I).EQ.JAUX)GO TO 21
      IAUX=I
   20 JAUX=JCA(I)
   21 IF(INDEX.LE.6.AND.JAUX.GE.(NS-1))GO TO 50
      IF(INDEX.EQ.4.AND.NL.EQ.0.AND.JAUX.GE.(NS-2))GO TO 50
      IF(INDEX.EQ.6.AND.NL.EQ.0.AND.JAUX.GE.(NS-2))GO TO 50
      IF(INDEX.GT.6.AND.JAUX.EQ.1)GO TO 50
      IAUX1=IAUX+1
      IAUX2=IAUX+2
      DO 22 I=IAUX1,KCA
      II=KCA+IAUX1-I
      JC(II+2)=JC(II)
   22 JCA(II+2)=JCA(II)
      KCA=KCA+2
      IF(INDEX.GE.7)GO TO 31
      JC(IAUX1)=JAUX+1
      JC(IAUX2)=JAUX+1
      JCA(IAUX1)=JC(IAUX1)
      JCA(IAUX2)=JC(IAUX2)
      IF(INDEX.GE.5)JC(IAUX1)=-JC(IAUX1)
      IF(INDEX.EQ.4.OR.INDEX.EQ.5)JC(IAUX2)=-JC(IAUX2)
       RETURN
   31 JC(IAUX1)=JAUX-1
      JC(IAUX2)=JAUX-1
      JCA(IAUX1)=JC(IAUX1)
      JCA(IAUX2)=JC(IAUX2)
      IF(INDEX.GE.9)JC(IAUX1)=-JC(IAUX1)
      IF(INDEX.EQ.8.OR.INDEX.EQ.9)JC(IAUX2)=-JC(IAUX2)
      RETURN
C
   50 INDEX=INDEX+1
      INEW=1
      IF(INDEX.GE.11)GO TO 80
   51 IF(INDEX.EQ.1.AND.IDP.EQ.0)GO TO 50
      IF(INDEX.EQ.2.AND.IDS.EQ.0)GO TO 50
      IF(INDEX.EQ.3.AND.IBP.EQ.0)GO TO 50
      IF(INDEX.EQ.4.AND.IBP.NE.2)GO TO 50
      IF(INDEX.EQ.5.AND.IBS.EQ.0)GO TO 50
      IF(INDEX.EQ.6.AND.IBS.NE.2)GO TO 50
      IF(INDEX.EQ.7.AND.IUP.EQ.0)GO TO 50
      IF(INDEX.EQ.8.AND.IUP.NE.2)GO TO 50
      IF(INDEX.EQ.9.AND.IUS.EQ.0)GO TO 50
      IF(INDEX.EQ.10.AND.IUS.NE.2)GO TO 50
      if(index.ge.7)irc=ir-1
      IMIN=IS
      IMAX=IRc
      IF(IS.GT.IRc)IMIN=IRc
      IF(IS.GT.IRc)IMAX=IS
      ISUM=IMAX-IMIN+2
      IF(INDEX.LE.2)GO TO 70
      INEW=0
      IF(INDEX.GT.6)GO TO 60
      ISM=IMAX-IS+1
      DO 55 I=1,ISM
      JCA(I)=IS+I-1
      JC(I)=JCA(I)
   55 IF(INDEX.GT.4)JC(I)=-JCA(I)
      ISM2=ISM+1
      DO 56 I=ISM2,ISUM
      JCA(I)=IMAX+ISM2-I
      JC(I)=JCA(I)
   56 IF(INDEX.EQ.4.OR.INDEX.EQ.5)JC(I)=-JCA(I)
      KCA=ISUM
      KC=1
      RETURN
   60 ISM=IS-IMIN+1
      DO 65 I=1,ISM
      JCA(I)=IS+1-I
      JC(I)=JCA(I)
   65 IF(INDEX.GT.8)JC(I)=-JCA(I)
      ISM2=ISM+1
      DO 66 I=ISM2,ISUM
      JCA(I)=IMIN+I-ISM2
      JC(I)=JCA(I)
   66 IF(INDEX.EQ.8.OR.INDEX.EQ.9)JC(I)=-JCA(I)
      KCA=ISUM
      KC=-1
      INEW=0
      RETURN
   70 KCA=IMAX-IMIN+1
      KC=1
      IF(IS.GE.IR)KC=-1
      DO 71 I=1,KCA
      JC(I)=IS+I-1
      IF(IS.GT.IR)JC(I)=IS-I+1
      JCA(I)=JC(I)
      IF(INDEX.EQ.2)JC(I)=-JC(I)
   71 CONTINUE
      RETURN
   80 IF(IREAD.EQ.0)GO TO 101
      READ(LIN,100)KC,KCA,(JC(I),I=1,KCA)
      WRITE(LOU,100)KC,KCA,(JC(I),I=1,KCA)
      IF(KCA.EQ.0)GOTO 101
      IF(KCA.GT.50)GO TO 99
      DO 81 I=1,KCA
   81 JCA(I)=IABS(JC(I))
      JAUX=JCA(1)
      IF(KC.EQ.0)RETURN
      IF(ir.gt.0.and.(JCA(KCA).GT.IR.OR.JCA(KCA).LT.(IR-1)))GO TO 99
      IF(KCA.EQ.1)RETURN
      ID=1
      IF(KC.LT.0)ID=-1
      DO 103 I=2,KCA
      JJ=JCA(I)
      IF(JJ.LE.0.OR.JJ.GE.NS)GO TO 99
      IF(JJ.EQ.JAUX)GO TO 105
      IF(ID.EQ.1.AND.JJ.EQ.(JAUX+1))GO TO 103
      IF(ID.EQ.(-1).AND.JJ.EQ.(JAUX-1))GO TO 103
      GO TO 99
  105 ID=-ID
  103 JAUX=JJ
      RETURN
   99 IFIN=-1
      RETURN
  101 IFIN=1
      RETURN
  100 FORMAT(24I3)
      END
C