C 
C     PROGRAM  A N R A Y, VERSION 4.30  (PRAHA, JUNE 2001)
C
C*******************************************************************
C
C     PROGRAM ANRAY IS DESIGNED FOR RAY, TRAVEL TIME AND
C     AMPLITUDE COMPUTATIONS IN 3D GENERAL ANISOTROPIC AND ISOTROPIC
C     LATERALLY VARYING LAYERED MEDIA. THE PROGRAM MAKES POSSIBLE
C     COMPUTATION OF RAYS SPECIFIED BY INITIAL ANGLES AT THE SOURCE,
C     I.E., INITIAL-VALUE RAY TRACING, OR RAYS STARTING FROM THE
C     SOURCE AND TERMINATING ON A VERTICAL OR SURFACE PROFILE, I.E.
C     BOUNDARY-VALUE RAY TRACING. RAY AMPLITUDES CAN BE COMPUTED
C     ALONG RAYS.
C
C*******************************************************************
C
C
      CHARACTER*80 MTEXT,FILEIN,FILEOU,FILE1,FILE2,FILE3,FILE4,FILE5
      CHARACTER*80 FILE6,FILE7
      DIMENSION Y(18)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      COMMON /AUXX/  MMX(20),MMY(20),MMXY(20)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST,
     1xprf,yprf
      COMMON /DENS/ RHO(20)
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,20),KINT(20),HHH(3,3),TMAX,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT,
     1   XINTA
      COMMON /ZERO/ RNULL
      COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP
      COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY
C
C**************************************************
C
      LIN=5
      LOU=6
      LU1=1
      LU2=2
      LU3=3
      LUBRD=7
      LUGRD=8
      LUIND=9
      LURAY=10
      FILEIN='anray.dat'
      FILEOU='anray.out'
      WRITE(*,'(2A)') ' SPECIFY NAMES OF INPUT AND OUTPUT FILES',
     1' LIN, LOU, LU1, LU2, LU3, LUBRD, LUGRD, LUIND, LURAY: '
      READ(*,*) FILEIN,FILEOU,FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7
      IF(FILE1.EQ.' ') LU1=0
      IF(FILE2.EQ.' ') LU2=0
      IF(FILE3.EQ.' ') LU3=0
      IF(FILE4.EQ.' ') LUBRD=0
      IF(FILE5.EQ.' ') LUGRD=0
      IF(FILE6.EQ.' ') LUIND=0
      IF(FILE7.EQ.' ') LURAY=0
      LOUT=LOU
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      IF(LU1.NE.0)OPEN(LU1,FILE=FILE1,FORM='FORMATTED')
      IF(LU2.NE.0)OPEN(LU2,FILE=FILE2,FORM='FORMATTED')
      IF(LU3.NE.0)OPEN(LU3,FILE=FILE3,FORM='FORMATTED')
      IF(LUBRD.NE.0)OPEN(LUBRD,FILE=FILE4,FORM='FORMATTED')
      IF(LUGRD.NE.0)OPEN(LUGRD,FILE=FILE5,FORM='FORMATTED')
      IF(LUIND.NE.0)OPEN(LUIND,FILE=FILE6,FORM='FORMATTED')
      IF(LURAY.NE.0)OPEN(LURAY,FILE=FILE7,FORM='FORMATTED')
C
C**************************************************
C
      WRITE(LOU,777)
 777  FORMAT(///,'***********************'
     1,//,'  PROGRAM   A N R A Y      ',//,
     2'***********************',//)
      NCODE=1
      MTEXT='ANRAY'
      INULL=4
      READ(LIN,*)MTEXT
      WRITE(LOU,115)MTEXT
      READ(LIN,*)INULL,ISURF
      RNULL=10.**(-INULL)
      WRITE(LOU,106)INULL,ISURF
C
C
C     SPECIFICATION OF THE MODEL
C
      CALL MODEL(MTEXT,LIN)
C
C     GENERATE FILE FOR PLOTTING VARIOUS CHARACTERISTIC SURFACES
C
      IF(LU3.NE.0)CALL SURFPL(LIN,LU3)
C
C     GENERATE FILE FOR VRML PLOTTING BOUNDARIES OF THE MODEL
C
      IF(LUBRD.NE.0)CALL BOX(BRD)
C
C     GENERATE FILE FOR PLOTTING RAYS
C
      IF(LURAY.NE.0)WRITE(LURAY,113)
      IF(LURAY.NE.0)WRITE(LURAY,105)
C
C     SPECIFICATION OF SYNTHETIC SEISMOGRAMS
C
    2 ICONT=1
      MEP=0
      MOUT=0
      MDIM=0
      METHOD=0
      MREG=0
      ITMAX=10
      IPOL=0
      IPREC=0
      IRAYPL=0
      IPRINT=0
      IAMP=0
      MTRNS=0
      ICOEF=0
      IRT=0
      ILOC=0
      MCOD=0
      MORI=0
      READ(LIN,*)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,
     1IPOL,IPREC,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI
      WRITE(LOU,102)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,
     1IPOL,IPREC,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI
      IF(ICONT.EQ.0)GO TO 99
C
C
      IF(MEP.NE.0.AND.MDIM.EQ.0)MDIM=1
      IVSP=0
      IF(ILOC.EQ.0)ITPR=3
      IF(ILOC.EQ.1)THEN
        IVSP=1
        ITPR=43
        MREG=1
      END IF
      IF(ILOC.GT.1)THEN
        ITPR=ILOC+100
        MREG=1
      END IF
C
      IF(MEP.EQ.0)THEN
        NDST=0
      END IF
C
      IF(MEP.EQ.1)THEN
        NDST=1
        READ(LIN,*)XREC,YREC
        WRITE(LOU,104)XREC,YREC
        GO TO 4
      END IF
      IF(MEP.LT.0)THEN
        NDST=-MEP
        PROF(1)=0.
        XPRF=0.
        YPRF=0.
        READ(LIN,*)PROF(1),(DST(I),I=1,NDST),XPRF,YPRF
        WRITE(LOU,104)PROF(1),(DST(I),I=1,NDST),XPRF,YPRF
        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)
      END IF
C
      IF(MEP.GT.0)THEN
        NDST=MEP
        READ(LIN,*)PROF(1),RMIN,RSTEP,XPRF,YPRF
        WRITE(LOU,104)PROF(1),RMIN,RSTEP,XPRF,YPRF
        DO 13 I=1,MEP
   13   DST(I)=RMIN+(I-1)*RSTEP
        IF(NDST.EQ.1)DST(2)=RMIN+RSTEP
      END IF
      PROF(2)=PROF(1)+1.
      NDSTP=1
C
      IF(IVSP.EQ.1.AND.NDST.NE.0)THEN
        READ(LIN,*)XVSP,YVSP
        WRITE(LOU,104)XVSP,YVSP
      END IF
C
    4 TSOUR=0.
      DT=1.
      AC=0.0001
      REPS=0.05
      PREPS=0.05
      READ(LIN,*)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS
      WRITE(LOU,104)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS
C
      IF(ABS(XPRF).LT..000001.AND.ABS(YPRF).LT..000001)THEN
        XPRF=XSOUR
        YPRF=YSOUR
      END IF
      IF(MEP.EQ.1)THEN
        XE=XREC-XPRF
        YE=YREC-YPRF
        RPRF=SQRT(XE*XE+YE*YE)
        XATAN=ATAN2(YE,XE)
        PROF(1)=XATAN
        RMIN=RPRF
        DST(1)=RMIN
        WRITE(LOU,104)RPRF,XATAN
        RSTEP=100.
        DST(2)=DST(1)+100.
        PROF(2)=PROF(1)+1.
        NDSTP=1
      END IF
C
      IF(IVSP.EQ.1.AND.NDST.NE.0)THEN
        XNRM=XVSP-XSOUR
        YNRM=YVSP-YSOUR
        AUX=SQRT(XNRM*XNRM+YNRM*YNRM)
        XNRM=XNRM/AUX
        YNRM=YNRM/AUX
        PROF(1)=ATAN2(YNRM,XNRM)
        PROF(2)=PROF(1)+1.
        XPRF=XSOUR
        YPRF=YSOUR
      END IF
      IF(MCOD.EQ.0)THEN
        READ(LIN,*)AMIN,ASTEP,AMAX
        WRITE(LOU,104)AMIN,ASTEP,AMAX
        BMIN=PROF(1)-.3
        BMAX=PROF(1)+.4
        BSTEP=.6
        READ(LIN,*)BMIN,BSTEP,BMAX
        WRITE(LOU,104)BMIN,BSTEP,BMAX
      END IF
      IF(MREG.EQ.0.AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)')
     1 ' COEFFICIENTS OF CONVERSION ARE APPLIED'
      IF(MREG.NE.0.AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)')
     1 ' COEFFICIENTS OF CONVERSION ARE *** NOT *** APPLIED'
      TMAX=10000.
      IND=-1
      NDER=1
      CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,AMIN1,BMIN,PX,PY,PZ,XX,YY,ZZ,T,
     1DT,AC)
      Y(1)=XSOUR
      Y(2)=YSOUR
      Y(3)=ZSOUR
      IF(IND.EQ.50)WRITE(LOU,111)IND
      IF(IND.EQ.50)GO TO 99
      LAY=IND
      ISOUR=IND
      ITYPE=3
      CALL PARDIS(Y,0)
      VP=SQRT(A11)
      IF(IRHO.EQ.0)RO=1.7+.2*VP
      IF(IRHO.EQ.1)RO=RHO(IND)
C
C     GENERATE FILE LU2 FOR SYNTHETIC SEISMOGRAM COMPUTATIONS
C
      IF(LU2.NE.0.AND.NDST.NE.0)THEN
        WRITE(LU2,115)MTEXT
        KSH=2
        WRITE(LU2,100)NDST,KSH,ILOC
        WRITE(LU2,104)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,RO
        IF(MEP.NE.1)WRITE(LU2,104)(DST(I),I=1,NDST)
        IF(MEP.EQ.1)WRITE(LU2,104)XREC,YREC,RPRF,XATAN
      END IF
C
C    LOOP FOR ELEMENTARY WAVES
C
   20 READ(LIN,*)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF)
      WRITE(LOU,100)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF)
      IF(KREF.EQ.0)GOTO 2
      IF(MOUT.NE.0)WRITE(LOU,107)
      WRITE(LOU,103)NCODE,KC,KREF,((CODE(I,K),K=1,2),I=1,KREF)
C
      IF(MCOD.NE.0)THEN
        READ(LIN,*)AMIN,ASTEP,AMAX
        WRITE(LOU,104)AMIN,ASTEP,AMAX
        BMIN=PROF(1)-.3
        BMAX=PROF(1)+.4
        BSTEP=.6
        READ(LIN,*)BMIN,BSTEP,BMAX
        WRITE(LOU,104)BMIN,BSTEP,BMAX
      END IF
C
C     GENERATE FILE LU1 FOR PLOTTING OF RAY DIAGRAMS,
C     TIME-DISTANCE AND AMPLITUDE-DISTANCE CURVES
C
      IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 21
      WRITE(LU1,100)ICONT,NDST,ILOC
      WRITE(LU1,104)RO
      NPN=2
      APN=0.
      WRITE(LU1,100)NPN,NPN,NPN
      WRITE(LU1,101)APN,APN,APN,APN,APN
      WRITE(LU1,101)APN,APN,APN,APN,APN
      WRITE(LU1,104)Xprf,Yprf,0.0,PROF(1)
      WRITE(LU1,104)(DST(I),I=1,NDST)
   21 CONTINUE
C
C
C     SEARCH FOR THE NUMBER OF THE ELEMENT OF THE RAY, STARTING FROM
C     WHICH THE WAVE DOES UNDERTAKE NEITHER REFLECTION NOR CONVERSION
C
      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,1)
      IC2=CODE(ICOD-1,1)
      IF((IC1-IC2).EQ.0)GO TO 35
      IC1=CODE(ICOD,2)
      IC2=CODE(ICOD-1,2)
      IF((IC1-IC2).NE.0)GO TO 35
   34 CONTINUE
   35 CONTINUE
      IF(MOUT.NE.0)WRITE(LOU,108)
C
C
      CALL RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX,AMIN,ASTEP,
     1AMAX,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCODE)
      IF(IND.EQ.14) WRITE(LOU,111) IND
      NCODE=NCODE+1
      GOTO 20
C
C     END OF LOOP FOR ELEMENTARY WAVES
C
C
  100 FORMAT(26I3)
  101 FORMAT(5E15.5)
  102 FORMAT(1H0,////,2X,26I3)
  103 FORMAT(4X,I4,9X,36I3)
  104 FORMAT(8F10.5)
  105 FORMAT('/')
  106 FORMAT(17I5)
  107 FORMAT(//2X,'INT.CODE',5X,'E X T E R N A L   C O D E')
  108 FORMAT(//)
  111 FORMAT(/2X,'IND=',I5,/)
  113 FORMAT(6H'RAYS')
  115 FORMAT(A)
C
   99 CONTINUE
      IF(LURAY.NE.0)WRITE(LURAY,105)
      IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)ICONT,ICONT
      IF(LU1.NE.0)REWIND LU1
      IF(LU2.NE.0)REWIND LU2
C
      STOP
      END
C
C     *********************************************************
C
      SUBROUTINE AMPL (AMPX,AMPY,AMPZ,UU)
C
C     ROUTINE FOR COMPUTING COMPLEX VECTORIAL RAY AMPLITUDES
C
C     OUTPUT PARAMETERS
C     AMPX(2),AMPY(2),AMPZ(2) - X,Y AND Z COMPONENTS OF COMPLEX
C     VECTORIAL RAY AMPLITUDES IN THE MODEL COORDINATES. FOR P WAVE
C     IN ANY MEDIUM AND FOR S WAVES IN AN ANISOTROPIC MEDIUM, I=1.
C     FOR S WAVE GENERATED IN AN ISOTROPIC MEDIUM, I=1,2. I=1 AND 2
C     CORRESPOND TO S WAVES SPECIFIED AT THE SOURCE BY VECTORS E1
C     E2. VECTORS E1 AND E2 TOGETHER WITH UNIT VECTOR TANGENT TO
C     THE RAY FORM A BASIS OF RAY CENTRED COORDINATE SYSTEM.
C     UU - PRODUCT OF RATIOS OF DENSITIES AND COSINES OF INCIDENCE
C     AND OF REFLECTION/TRANSMISSION AT POINTS WHERE THE RAY CROSSES
c     INTERFACES.
C
C     CALLED FROM: RECEIV
C     ROUTINES CALLED: POLAR,TRANSL,COEF
C
      DIMENSION Y(18),UN(3),POLD(3),PNEW(3)
      COMPLEX  AMPX(2),AMPY(2),AMPZ(2),CR(3),UC(3),STU(6)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,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 /DENS/ RHO(20)
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /RAY2/  DRY(3,2000)
C
      KSS=1
      ISHEAR=0
      ITYPE=CODE(1,2)
      IF(IANI(ISOUR).EQ.0.AND.ITYPE.NE.3)then
        ISHEAR=1
        itype=1
      end if
      ITP=ITYPE
      do 1 i=1,2
      AMPX(i)=CMPLX(0.,0.)
      AMPY(i)=CMPLX(0.,0.)
      AMPZ(i)=CMPLX(0.,0.)
    1 continue
C
 3000 NN=N
      IDD=0
      N2=0
      N1=1
      IRE=IREF
      AV=1.
C
C  SPECIFICATION OF DISPLACEMENT VECTOR AT SOURCE
C  IN RAY CENTERED COORDINATES
C
      DO 5 I=1,3
      CR(I)=(0.,0.)
    5 CONTINUE
      CR(ITP)=(1.,0.)
      IREF1=IREF-1
      IF(IREF1.EQ.0) GOTO 100
C
C  LOOP OVER INTERFACES
C
      DO 10 I=1,IREF1
      IREF=I
      IF(KC.NE.0) ITYPE=CODE(IREF,2)
      N=KINT(IREF)
      IF(N.EQ.0) THEN
        IDD=1
        GO TO 10
      ELSE
        N1=N2+1
        N2=N
        IF(IDD.NE.0) N2=-N2
        IDD=0
C
C  COMPUTATION OF POLARIZATION VECTORS
C  CONSIDERED POLARIZATION VECTOR(S) ARE STORED IN CORRESPONDING
C  COLUMNS OF THE MATRIX HHH. OTHER COLUMNS ARE ZERO.
C
        CALL POLAR(N1,N2,NN,IREF)
      END IF
      DO 20 K=1,6
      Y(K)=AY(K+1,N)
  20  CONTINUE
      IF(IAMP.GT.0)WRITE(LOU,'(a,2i5,6f10.5)')' AMPL:I,N,Y',I,N,
     1(Y(L),L=1,6)
      DO 30 K=1,3
      POLD(K)=Y(K+3)
      PS(K,7,IREF)=Y(K+3)
  30  CONTINUE
      DO 40 K=1,3
      UN(K)=DS(K,IREF)
  40  CONTINUE
      LAY=IS(1,IREF)
      ITRANS=IS(2,IREF)
      ITR1=ITRANS
      IF(UN(3).GT.0.0) GOTO 50
C
C  RAY STRIKING THE INTERFACE FROM ABOVE
C
      IF(ITRANS.EQ.0) THEN
        LAY=LAY+1
        ITRANS=1
        GOTO 70
      END IF
      IF(ITRANS.GT.0) THEN
        LAY=LAY-1
        ITRANS=0
        GOTO 70
      END IF
C
C  RAY STRIKING THE INTERFACE FROM BELOW
C
  50  IF(ITRANS.EQ.0) THEN
        LAY=LAY-1
        ITRANS=1
        GOTO 70
      END IF
      IF(ITRANS.GT.0) THEN
        LAY=LAY+1
        ITRANS=0
        GOTO 70
      END IF
C
C  SLOWNESS VECTORS ON THE SIDE OF THE INTERFACE WHERE GENERATED
C  WAVE PROPAGATES WERE DETERMINED DURING THE CALL OF TRANSL IN THE
C  ROUTINE OUT. HERE REMAINING SLOWNESS VECTORS ON THE OTHER SIDE
C  OF THE INTERFACE ARE DETERMINED
C
C  REDEFINITION OF IREF FOR CALL OF ROUTINE TRANSL
C
  70  IF(LAY.EQ.0) THEN
        DO 71 K=4,6
        DO 71 L=1,3
        PS(L,K,IREF)=CMPLX(0.,0.)
  71    CONTINUE
        GO TO 75
      END IF
      IREF=IREF+1
      CALL TRANSL(Y,POLD,PNEW,UN,ITRANS,0)
      IF(IND.EQ.10)RETURN
      IREF=IREF-1
  75  IF(IAMP.ne.0)then
        WRITE(LOU,'(a)')' REFLECTED/TRANSMITTED SLOWNESS VECTORS'
        WRITE(LOU,'(6f12.6)')((PS(L,K,IREF),L=1,3),k=1,6)
      end if
      AV1=(DS(11,IREF)*ds(10,iref))/(DS(8,IREF)*ds(7,iref))
      AV=AV*AV1
      IF(IAMP.GT.0) THEN
        WRITE(LOU,'(A)') 'ROI,ROG,UNVGI,UNVGG,AV1,AV'
        WRITE(LOU,'(6F10.5)') DS(8,IREF),
     1  DS(11,IREF),ds(7,iref),ds(10,iref),AV1,AV
        WRITE(LOU,'(A,/,6F12.5,/,3(3F12.5/))') ' CR,HHH',
     2  CR,((HHH(J,K),J=1,3),K=1,3)
      END IF
C
C  COMPUTATION OF AMPLITUDE COEFFICIENTS OF REFLECTED/TRANSMITTED WAVES
C
C
C  COMPUTATION OF CARTESIAN COMPONENTS OF INCIDENT DISPLACEMENT VECTOR
C
      DO 87 K=1,3
      STU(K)=CMPLX(0.,0.)
      DO 87 J=1,3
      STU(K)=HHH(J,K)*CR(J)+STU(K)
  87  CONTINUE
      IF(IAMP.GT.0)WRITE(LOU,'(A,6F10.5)') ' STU',(STU(K),K=1,3)
      IF(KC.NE.0)ITYPE=CODE(IREF+1,2)
      CALL COEF(STU,CR,ITR1)
      IF(IND.EQ.11)RETURN
      BCR=SQRT(REAL(CR(1)*CONJG(CR(1))+CR(2)*CONJG(CR(2))
     1             +CR(3)*CONJG(CR(3))))
      IF(BCR.LT.1.E-10) THEN
        DO 88 K=1,3
        UC(K)=(0.,0.)
  88    CONTINUE
        GOTO 130
      END IF
  10  CONTINUE
C
C  END OF LOOP OVER INTERFACES
C
C  TERMINATION POINT
C
  100 N1=N2+1
      N2=NN
      IF(KC.NE.0)ITYPE=CODE(IRE,2)
      CALL POLAR(N1,N2,NN,IRE)
C
C COMPUTATION OF CARTESIAN COMPONENTS OF INCIDENT DISPLACEMENT VECTOR
C
      DO 107 K=1,3
      STU(K)=CMPLX(0.,0.)
      DO 107 J=1,3
      STU(K)=HHH(J,K)*CR(J)+STU(K)
 107  CONTINUE
      IF(IAMP.GT.0)WRITE(LOU,'(A,6F10.5)') ' STU',(STU(K),K=1,3)
      DO 105 K=1,6
      Y(K)=AY(K+1,Nn)
      IF(K.LE.3)GO TO 105
      PS(K-3,7,IRE)=Y(K)
  105 CONTINUE
      n=nn
      iref=ire
      IF(MREG.NE.0) THEN
        UC(1)=STU(1)
        UC(2)=STU(2)
        UC(3)=STU(3)
        GOTO 110
      END IF
      IF(IAMP.GT.0)then
        WRITE(LOU,'(a)')
     1  ' REFLECTED SLOWNESS VECTORS AT TERMINATION POINT'
        WRITE(LOU,'(6f12.6)')((PS(L,K,IRE),L=1,3),K=1,3)
      end if
C
C COMPUTATION OF CONVERSION COEFFICIENTS
      KTR=999
      CALL COEF(STU,UC,KTR)
      IF(IND.EQ.11)RETURN
  110 VPEND=1./SQRT(Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6))
      DO 120 K=1,3
      Y(K)=AY(K+4,1)
  120 CONTINUE
      VP0=1./SQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3))
      RHO0=0.2*SQRT(AY(8,1))+1.7
      IF(IRHO.NE.0) RHO0=RHO(ISOUR)
      RHEND=0.2*SQRT(AY(8,NN))+1.7
      IF(IRHO.NE.0) RHEND=RHO(LAY)
      AV=AV*VP0*RHO0
      AV=AV/(VPEND*RHEND)
      UU=SQRT(ABS(AV))
      IF(IAMP.GT.0)
     1WRITE(LOU,'(A,4F12.6)')'VP0,RH0,VPEND,RHEND',Vp0,RHO0,VpEND,RHEND
 130  CONTINUE
      n=nn
      iref=ire
      AMPX(KSS)=UC(1)
      AMPY(KSS)=UC(2)
      AMPZ(KSS)=UC(3)
      IF(ISHEAR.NE.0.AND.KSS.NE.2) THEN
        KSS=2
        ITP=2
        GOTO 3000
      END IF
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE APPROX(X,Y,YD,KDIM)
C
C     THE ROUTINE PERFORMS THIRD-ORDER INTERPOLATION BETWEEN POINTS
C     YOLD AND YNEW PARAMETERIZED BY AN INDEPENDENT VARIABLE X.
C     DOLD, DNEW ARE THE FIRST DERIVATIVES OF Y WITH RESPECT
C     TO X AT THE POINTS YOLD AND YNEW.
C
      DIMENSION Y(18),YD(18)
      COMMON/APPR/ XOLD,XNEW,YOLD(18),DOLD(18),YNEW(18),DNEW(18)
C
      A=(X-XNEW)/(XNEW-XOLD)
      AUX=A+1.
      A1=(2.*A+3.)*A*A
      A2=1.-A1
      B1=AUX*A*(X-XNEW)
      B2=AUX*A*(X-XOLD)
      AD1=6.*A*AUX/(XNEW-XOLD)
      AD2=-AD1
      BD1=A*(3.*A+2.)
      BD2=AUX*(3.*A+1.)
      DO 1 I=1,KDIM
      Y(I)=A1*YOLD(I)+A2*YNEW(I)+B1*DOLD(I)+B2*DNEW(I)
      YD(I)=AD1*YOLD(I)+AD2*YNEW(I)+BD1*DOLD(I)+BD2*DNEW(I)
    1 CONTINUE
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE BIAP(MX1,MX,MY1,MY,MXY1)
C
      DIMENSION X(200),FX(200),V(1000)
      COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000)
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT,
     1   XINTA
      EQUIVALENCE(Z(1),V(1))
C
C     ROUTINE DETERMINING THE COEFFICIENTS
C     OF BICUBIC SPLINE INTERPOLATION
C
      DO 1 J=1,MX
      L=MX1+J-1
    1 X(J)=SX(L)
      DO 3 I=1,MY
      DO 2 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    2 FX(J)=V(K)
      CALL SPLIN(X,FX,1,MX)
      DO 3 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    3 A20(K)=FX(J)
C
      DO 4 I=1,MY
      L=MY1+I-1
    4 X(I)=SY(L)
      DO 6 J=1,MX
      DO 5 I=1,MY
      K=MXY1+(J-1)*MY+I-1
    5 FX(I)=V(K)
      CALL SPLIN(X,FX,1,MY)
      DO 6 I=1,MY
      K=MXY1+(J-1)*MY+I-1
    6 A02(K)=FX(I)
C
      DO 7 J=1,MX
      L=MX1+J-1
    7 X(J)=SX(L)
      DO 9 I=1,MY
      DO 8 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    8 FX(J)=A02(K)
      CALL SPLIN(X,FX,1,MX)
      DO 9 J=1,MX
      K=MXY1+(J-1)*MY+I-1
    9 A22(K)=FX(J)
C
      RETURN
      END
C
C
C     *********************************************************
C
      SUBROUTINE CHRM(Y)
C
C  ROUTINE FOR THE COMPUTATION OF THE ELEMENTS OF THE CHRISTOFFEL
C  MATRIX FOR AN ARBITRARY ANISOTROPIC MEDIUM
C
      DIMENSION Y(18)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /DJK/  D11,D12,D13,D22,D23,D33,DTR
      COMMON /GAM/  C11,C12,C13,C22,C23,C33
C
      P1=Y(4)
      P2=Y(5)
      P3=Y(6)
      P2P3=P2*P3
      P1P2=P1*P2
      P1P3=P1*P3
      P1P1=P1*P1
      P2P2=P2*P2
      P3P3=P3*P3
      C11=P1P1*A11+P2P2*A66+P3P3*A55
     1+2.*(P2P3*A56+P1P3*A15+P1P2*A16)
      C22=P1P1*A66+P2P2*A22+P3P3*A44
     1+2.*(P2P3*A24+P1P3*A46+P1P2*A26)
      C33=P1P1*A55+P2P2*A44+P3P3*A33
     1+2.*(P2P3*A34+P1P3*A35+P1P2*A45)
      C23=P1P1*A56+P2P2*A24+P3P3*A34
     1   +P2P3*A2344+P1P3*A3645+P1P2*A2546
      C13=P1P1*A15+P2P2*A46+P3P3*A35
     1   +P2P3*A3645+P1P3*A1355+P1P2*A1456
      C12=P1P1*A16+P2P2*A26+P3P3*A45
     1   +P2P3*A2546+P1P3*A1456+P1P2*A1266
      C11N=C11-1.
      C22N=C22-1.
      C33N=C33-1.
      C23SQ=C23*C23
      C13SQ=C13*C13
      C12SQ=C12*C12
      D11=C22N*C33N-C23SQ
      D22=C11N*C33N-C13SQ
      D33=C11N*C22N-C12SQ
      D12=C13*C23-C12*C33N
      D13=C12*C23-C13*C22N
      D23=C12*C13-C23*C11N
      DTR=D11+D22+D33
      IF(ABS(DTR).LT.0.00001)THEN
        IND=10
      END IF
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE CHRM1(C,PN,UN)
C
C  ROUTINE FOR THE COMPUTATION OF THE ELEMENTS OF THE CHRISTOFFEL
C  MATRIX FOR AN ARBITRARY ANISOTROPIC MEDIUM
C
      DIMENSION C(3,3),PN(3),UN(3)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
C
      P1=PN(1)
      P2=PN(2)
      P3=PN(3)
      U1=UN(1)
      U2=UN(2)
      U3=UN(3)
      P2U3=P2*U3
      P3U2=P3*U2
      P1U2=P1*U2
      P2U1=P2*U1
      P1U3=P1*U3
      P3U1=P3*U1
      P1U1=P1*U1
      P2U2=P2*U2
      P3U3=P3*U3
      C(1,1)=P1U1*A11+P2U2*A66+P3U3*A55
     1+(P2U3+P3U2)*A56+(P1U3+P3U1)*A15+(P1U2+P2U1)*A16
      C(2,2)=P1U1*A66+P2U2*A22+P3U3*A44
     1+(P2U3+P3U2)*A24+(P1U3+P3U1)*A46+(P1U2+P2U1)*A26
      C(3,3)=P1U1*A55+P2U2*A44+P3U3*A33
     1+(P2U3+P3U2)*A34+(P1U3+P3U1)*A35+(P1U2+P2U1)*A45
      C(2,3)=P1U1*A56+P2U2*A24+P3U3*A34
     1+0.5*((P2U3+P3U2)*A2344+(P1U3+P3U1)*A3645+(P1U2+P2U1)*A2546)
      C(1,3)=P1U1*A15+P2U2*A46+P3U3*A35
     1+0.5*((P2U3+P3U2)*A3645+(P1U3+P3U1)*A1355+(P1U2+P2U1)*A1456)
      C(1,2)=P1U1*A16+P2U2*A26+P3U3*A45
     1+0.5*((P2U3+P3U2)*A2546+(P1U3+P3U1)*A1456+(P1U2+P2U1)*A1266)
      C(2,1)=C(1,2)
      C(3,2)=C(2,3)
      C(3,1)=C(1,3)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE CHRM2(Y,G,i)
C
C     EVALUATES ELEMENTS OF THE CHRISTOFFEL MATRIX
C
      DIMENSION a(21),Y(18),G(3,3)
      COMMON/GAM/G11,G12,G13,G22,G23,G33
      COMMON /APROX1/ e(21,10)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
C
      DO 1 J=1,21
      A(J)=E(J,I)
    1 CONTINUE
      P1=Y(4)
      P2=Y(5)
      P3=Y(6)
      P11=P1*P1
      P12=P1*P2
      P13=P1*P3
      P22=P2*P2
      P23=P2*P3
      P33=P3*P3
      G11=A(1)*P11+A(21)*P22+A(19)*P33+
     1    2.*(A(6)*P12+A(5)*P13+A(20)*P23)
      G22=A(21)*P11+A(7)*P22+A(16)*P33+
     1    2.*(A(11)*P12+A(18)*P13+A(9)*P23)
      G33=A(19)*P11+A(16)*P22+A(12)*P33+
     1    2.*(A(17)*P12+A(14)*P13+A(13)*P23)
      G12=A(6)*P11+A(11)*P22+A(17)*P33+
     1    (A(21)+A(2))*P12+(A(20)+A(4))*P13+(A(10)+A(18))*P23
      G13=A(5)*P11+A(18)*P22+A(14)*P33+
     1    (A(20)+A(4))*P12+(A(19)+A(3))*P13+(A(17)+A(15))*P23
      G23=A(20)*P11+A(9)*P22+A(13)*P33+
     1    (A(10)+A(18))*P12+(A(17)+A(15))*P13+(A(16)+A(8))*P23
      G(1,1)=G11
      G(1,2)=G12
      G(1,3)=G13
      G(2,1)=G12
      G(2,2)=G22
      G(2,3)=G23
      G(3,1)=G13
      G(3,2)=G23
      G(3,3)=G33
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE PCHRM(Y,G,L,I)
C
C     EVALUATES FIRST DERIVATIVES OF ELEMENTS OF CHRISTOFFEL MATRIX
C     WITH RESPECT TO THE L-TH COMPONENT OF THE SLOWNESS VECTOR
C
      DIMENSION A(21),Y(18),G(3,3)
      COMMON /APROX1/ E(21,10)
C
      DO 1 J=1,21
      A(J)=E(J,I)
    1 CONTINUE
      P1=Y(4)
      P2=Y(5)
      P3=Y(6)
      IF(L.EQ.1)THEN
        G(1,1)=2.*(A(1)*P1+A(6)*P2+A(5)*P3)
        G(2,2)=2.*(A(21)*P1+A(11)*P2+A(18)*P3)
        G(3,3)=2.*(A(19)*P1+A(17)*P2+A(14)*P3)
        AUX=2.*A(6)*P1+(A(21)+A(2))*P2+(A(20)+A(4))*P3
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(5)*P1+(A(20)+A(4))*P2+(A(19)+A(3))*P3
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(20)*P1+(A(10)+A(18))*P2+(A(17)+A(15))*P3
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.2)THEN
        G(1,1)=2.*(A(6)*P1+A(21)*P2+A(20)*P3)
        G(2,2)=2.*(A(11)*P1+A(7)*P2+A(9)*P3)
        G(3,3)=2.*(A(17)*P1+A(16)*P2+A(13)*P3)
        AUX=2.*A(11)*P2+(A(21)+A(2))*P1+(A(10)+A(18))*P3
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(18)*P2+(A(20)+A(4))*P1+(A(17)+A(15))*P3
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(9)*P2+(A(10)+A(18))*P1+(A(16)+A(8))*P3
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.3)THEN
        G(1,1)=2.*(A(5)*P1+A(20)*P2+A(19)*P3)
        G(2,2)=2.*(A(18)*P1+A(9)*P2+A(16)*P3)
        G(3,3)=2.*(A(14)*P1+A(13)*P2+A(12)*P3)
        AUX=2.*A(17)*P3+(A(20)+A(4))*P1+(A(10)+A(18))*P2
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(14)*P3+(A(19)+A(3))*P1+(A(17)+A(15))*P2
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(13)*P3+(A(17)+A(15))*P1+(A(16)+A(8))*P2
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE PPCHRM(G,L,M,i)
C
C     EVALUATES SECOND DERIVATIVES OF ELEMENTS OF CHRISTOFFEL MATRIX
C     WITH RESPECT TO THE L-TH AND M-TH COMPONENTS OF THE SLOWNESS
C     VECTOR
C
      DIMENSION a(21),G(3,3)
      COMMON /APROX1/ e(21,10)
C
      do 1 j=1,21
      a(j)=e(j,i)
    1 continue
      IF(L.EQ.1.AND.M.EQ.1)THEN
        G(1,1)=2.*A(1)
        G(2,2)=2.*A(21)
        G(3,3)=2.*A(19)
        AUX=2.*A(6)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(5)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(20)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.2.AND.M.EQ.2)THEN
        G(1,1)=2.*A(21)
        G(2,2)=2.*A(7)
        G(3,3)=2.*A(16)
        AUX=2.*A(11)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(18)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(9)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF(L.EQ.3.AND.M.EQ.3)THEN
        G(1,1)=2.*A(19)
        G(2,2)=2.*A(16)
        G(3,3)=2.*A(12)
        AUX=2.*A(17)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=2.*A(14)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=2.*A(13)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF((L.EQ.1.AND.M.EQ.2).OR.(L.EQ.2.AND.M.EQ.1))THEN
        G(1,1)=2.*A(6)
        G(2,2)=2.*A(11)
        G(3,3)=2.*A(17)
        AUX=A(21)+A(2)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=A(20)+A(4)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=A(10)+A(18)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF((L.EQ.1.AND.M.EQ.3).OR.(L.EQ.3.AND.M.EQ.1))THEN
        G(1,1)=2.*A(5)
        G(2,2)=2.*A(18)
        G(3,3)=2.*A(14)
        AUX=A(20)+A(4)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=A(19)+A(3)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=A(17)+A(15)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      IF((L.EQ.2.AND.M.EQ.3).OR.(L.EQ.3.AND.M.EQ.2))THEN
        G(1,1)=2.*A(20)
        G(2,2)=2.*A(9)
        G(3,3)=2.*A(13)
        AUX=A(10)+A(18)
        G(1,2)=AUX
        G(2,1)=AUX
        AUX=A(17)+A(15)
        G(1,3)=AUX
        G(3,1)=AUX
        AUX=A(16)+A(8)
        G(2,3)=AUX
        G(3,2)=AUX
      END IF
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE FACETS(N1,N2,NSRF)
      INTEGER LU,N1,N2,NSRF
C
C Subroutine FACETS writes the index file listing the vertices of each
C tetragon covering the structural interface.  The vertices are assumed
C to be stored in a separate file, with inner loop over N1 points along
C the first horizontal axis, middle loop over N2 points along the second
C horizontal axis and outer loop over the surfaces.  The vertices are
C indexed by positive integers according to their order in the vertex
C file.
C
C Input:
C     LU...   Logical unit number connected to the output file to be
C             written by this subroutine.
C     N1...   Number of points along the first horizontal axis.
C     N2...   Number of points along the second horizontal axis.
C     NSRF... Number of interfaces.
C The input parameters are not altered.
C
C No output.
C
C Output index file with the tetragons:
C For each tetragon, a line containing I1,I2,I3,I4,/
C     I1,I2,I3,I4... Indices of the vertices of the tetragon.
C             The vertices are indexed by positive integers according to
C             their order in the respective vertex file.
C     /...    List of vertices is terminated by a slash.
C
C Date: 1999, October 4
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      CHARACTER*9 FORMAT
      INTEGER I1,I2,ISRF
      COMMON/VRML/LUBRD,LUGRD,LU,LURAY
C
      IF(LU.EQ.0)RETURN
C     Setting output format:
      FORMAT='(4(I0,A))'
      I1=INT(ALOG10(FLOAT(N1*N2*NSRF)+0.5))+1
      FORMAT(5:5)=CHAR(ICHAR('0')+I1)
C
C     Writing the file:
      DO 33 ISRF=0,N1*N2*(NSRF-1),N1*N2
        DO 32 I2=ISRF,ISRF+N1*(N2-2),N1
          DO 31 I1=I2+1,I2+N1-1
            WRITE(LU,FORMAT) I1,' ',I1+1,' ',I1+1+N1,' ',I1+N1,' /'
   31     CONTINUE
   32   CONTINUE
   33 CONTINUE
C
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE BOX(BRD)
C
      DIMENSION BRD(6)
      COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY
C
        WRITE(LUBRD,109)
        WRITE(LUBRD,105)
        I=1
        WRITE(LUBRD,112)I
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(6)
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6)
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5)
        WRITE(LUBRD,105)
        I=2
        WRITE(LUBRD,112)I
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5)
        WRITE(LUBRD,105)
        I=3
        WRITE(LUBRD,112)I
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5)
        WRITE(LUBRD,110)BRD(2),BRD(4),BRD(5)
        WRITE(LUBRD,110)BRD(2),BRD(4),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5)
        WRITE(LUBRD,105)
        I=4
        WRITE(LUBRD,112)I
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(4),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(4),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5)
        WRITE(LUBRD,105)
        I=1
        WRITE(LUBRD,112)I
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5)
        WRITE(LUBRD,110)BRD(2),BRD(4),BRD(5)
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5)
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5)
        WRITE(LUBRD,105)
        I=1
        WRITE(LUBRD,112)I
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6)
        WRITE(LUBRD,110)BRD(1),BRD(4),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(4),BRD(6)
        WRITE(LUBRD,110)BRD(2),BRD(3),BRD(6)
        WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6)
        WRITE(LUBRD,105)
        WRITE(LUBRD,105)
C
  105 FORMAT('/')
  109 FORMAT(25H'BOUNDARIES OF THE MODEL')
  110 FORMAT(3(F10.5,1X),'/')
  112 FORMAT(6H'BOUND,I1,1H',1X,'/')
C
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'a2.for'
C     a2.for
      INCLUDE 'a3.for'
C     a3.for
      INCLUDE 'a4.for'
C     a4.for
      INCLUDE 'a5.for'
C     a5.for
C
C Interpolation method:
C Include just one of the following files 'mod*.for':
C (a) Isosurface interpolation:
      INCLUDE 'modis.for'
C     modis.for
C (b) (Bi-)(tri-)cubic B-spline interpolation:
*      INCLUDE 'modbs.for'
C     modbs.for
C
C=======================================================================
C