      SUBROUTINE INPUT(TITLE,FREQ,M,N,IPRT,IOUT,BSOUT,NSDEP,NDEP,
     1                 NATEN,NDPTH,IEIG,IMAT,RANGE,DEPTH,NREG,MAMP,
     2                 ZS,ZSD,RMIN,RMAX,RINC,ZR,NDINC,MD,NRX,
     3                 RANG,DPTH,IRLIN,IREIG,MX,  
     4                 DEPW,SVPW,DBPWLW,CGRADW,CINTW,RHOW,GRHOW,
     5                 DEPB,SVPB,DBPWLB,CGRADB,CINTB,RHOB,GRHOB,
     6                 NEWSVP,LUSVP)
C
C     THIS SUBROUTINE READS THE INPUT, GENERATES THE INDEX OF REFRACTION
C     TABLES, STORES THE INDEX OF REFRACTION TABLES ON DISK (LUSVP), 
C     INTERPOLATES OR EXTENDS THE BATHYMETRY, COMPUTES AN ESTIMATE 
C     OF THE SUB BOTTOM REFLECTED FIELD (IN THE TWO ISOSPEED LAYERS CASE
C     ONLY) AND GENERATES THE SOURCE AND RECEIVER DEPTHS.
C       
C     IMAT IS THE FLAG WHICH CONTROLS HOW THE COMPLEX EIGENVALUES
C     ARE FOUND AND HOW THE DEPTH FUNCTIONS ARE CONSTRUCTED.
C     IF NSVP(1) .LT. 2 (THERE ARE TWO ISOSPEED LAYERS) THEN IMAT=0.
C     THE EIGENVALUES ARE FOUND IN SUBROUTINE SEARCH AND THE ANALYTIC
C     FORM OF THE DEPTH FUNCTIONS IS USED. IF NSVP(1) .GE. 2 A SOUND
C     SPEED PROFILE IN THE WATER AND A SOUND SPEED AND ATTUNUATION 
C     PROFILE IN THE BOTTOM ARE INPUT (BOTH LAYERS MAY BE NONHOMO-
C     GENEOUS) THEN IMAT=1. THE EIGENVALUES ARE FOUND BY THE GALERKIN
C     MATRIX METHOD IN SUBROUTINE GALRKN. THE DEPTH FUNCTIONS ARE 
C     CONSTRUCTED FROM THE MATRIX EIGENVECTORS.     
C     
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER BSOUT
      DIMENSION RANGE(NRX),DEPTH(NRX)
      DIMENSION ZSD(MD),ZR(MD)
      DIMENSION DEPW(MX),SVPW(MX),DBPWLW(MX),RHOW(MX),GRHOW(MX)
      DIMENSION DEPB(MX),SVPB(MX),DBPWLB(MX),RHOB(MX),GRHOB(MX)
      DIMENSION RANG(NRX),DPTH(NRX),IRLIN(NRX),IREIG(NRX)
C
      COMPLEX*16 CGRADW(MX),CINTW(MX)
      COMPLEX*16 CGRADB(MX),CINTB(MX)
      COMPLEX*16 CI,ONE,FKH,FKVS,FKVW
      COMPLEX*16 CXP
      COMPLEX*16 CDSRT
      COMPLEX*16 CDEM
      COMPLEX*16 CKB,PREFL,CPG
C
      CHARACTER*10 TITLE(8)
      CHARACTER*3 NEWSVP(NRX) 
C
      COMMON /BLKEVN/ HB,CW1,CB,FKW,FKB,ROHW,ROHB,ATEN
      COMMON/FLAGS/IGEOM,NEWSBC
C
      PI=2.*ACOS(0.0)
      CI=DCMPLX(0.0,1.0)
      ONE=DCMPLX(1.0,0.0)
      NREG=0
C
C     READ INPUT 
C
      WRITE(9,50)
   50 FORMAT(1H1)
      READ(5,100) (TITLE(I),I=1,8)
  100 FORMAT(8A10)
      WRITE(6,110) (TITLE(I),I=1,8)
      WRITE(9,110) (TITLE(I),I=1,8)
  110 FORMAT(1X,8A10/)
C
      READ(5,*) HB,FREQ,ROHW,ROHB,DBPWL
      WRITE(6,120) HB,FREQ,ROHW,ROHB,DBPWL
      WRITE(9,120) HB,FREQ,ROHW,ROHB,DBPWL
  120 FORMAT(1X,'HB=',F10.2,2X,'FREQ=',F10.2,2X,'ROHW=',F10.2,2X,
     1       'ROHB=',F10.2,2X,/,1X,'DBPWL=',E14.6,/)
C
      WRITE(9,125)
  125 FORMAT(/)
      READ(5,*) CB,CW1
      WRITE(6,130) CB,CW1
      WRITE(9,130) CB,CW1
  130 FORMAT(1X,'CB=',4PE16.7,2X,'CW=',4PE16.7,2X,/)
      WRITE(9,125)
C
C     FIND THE WAVE NUMBER IN THE BOTTOM
C
      FKB=2.0*PI*FREQ/CB
      WRITE(6,140) FKB   
      WRITE(9,140) FKB   
  140 FORMAT(1X,/,1X,'THE WAVE NUMBER IN THE BOTTOM= ',1P,
     1       E13.5,1X,'1/M'/)
C
C     CONVERT THE ATTENUATION IN THE BOTTOM (DBPWL) TO AN 
C     IMAGINARY PART OF THE SQUARE OF THE INDEX OF REFRACTION 
C     (ATEN).
C
      CALL ESPI(DBPWL,FREQ,CB,FKB,ATEN)
C
C
      WRITE(6,150) ATEN   
      WRITE(9,150) ATEN   
  150 FORMAT(1X,/,'THE IMAGINARY PART OF THE SQUARE OF THE ',
     1       'COMPLEX INDEX OF REFRACTION IN AN',/,1X,
     2       'ISOSPEED BOTTOM BASED ON DBPWL= ',1P,E13.5,///)
C
C
      READ(5,*) M,NATEN,IEIG,NDPTH,IGEOM,NEWSBC
      IF(IGEOM .EQ. 0)   THEN
      WRITE(6,155)
      WRITE(9,155)
      ELSE
      WRITE(6,160)
      WRITE(9,160)
      IGEOM=1
      NEWSBC=1
      END IF
  155 FORMAT(1X,//,' ********  CYLINDRICAL GEOMETRY PROBLEM *******',/)   
  160 FORMAT(1X,//,' ********  PLANE GEOMETRY PROBLEM ******* ',/)   
C
      IF(NEWSBC .EQ. 0)   THEN
      WRITE(6,165)
      WRITE(9,165)
      ELSE
      WRITE(6,170)
      WRITE(9,170)
      NEWSBC=1
      END IF
  165 FORMAT(1X,/,'  ********* ',/,
     & '  OLD INITIAL SOURCE BOUNDARY CONDITIONS PROBLEM',
     & /,'  ********* ',/)
  170 FORMAT(1X,/,' ********',
     & '  NEW INITIAL SOURCE BOUNDARY CONDITIONS PROBLEM',
     & '  ******** ',/)
C
      WRITE(9,125)
      WRITE(6,180) M,NATEN,IEIG,NDPTH
      WRITE(9,180) M,NATEN,IEIG,NDPTH
  180 FORMAT(1X,2HM=,I5,3X,6HNATEN=,I5,3X,5HIEIG=,I5,3X,
     1       6HNDPTH=,I5,/)
C
C     READ SOURCE SPECIFICATIONS
C
      WRITE (6,125)
      READ(5,*) NSDEP,ZSMIN,ZSINC
      WRITE(6,190) NSDEP,ZSMIN,ZSINC
      WRITE(9,190) NSDEP,ZSMIN,ZSINC
  190 FORMAT(1X,'NSDEP= ',I5,5X,'ZSMIN=',F10.2,5X,'ZSINC= ',F10.2)
C
      IF(NSDEP .GT. MD) NSDEP=MD
      DO 200 J=1,NSDEP
      ZSD(J)=ZSMIN+(J-1)*ZSINC
  200 CONTINUE
C
C     FIND THE MIDPOINT OF THE SOURCE ARRAY
C
      ZS=(ZSD(1)+ZSD(NSDEP))/2.0D0
C
C     READ RECEIVER SPECIFICATIONS
C
      WRITE (6,125)
      READ(5,*) NDEP,ZMIN,ZINC
      WRITE(6,210) NDEP,ZMIN,ZINC
      WRITE(9,210) NDEP,ZMIN,ZINC
  210 FORMAT(1X,'NDEP= ',I5,5X,'ZMIN=',F10.2,5X,'ZINC= ',F10.2)
C
      IF(NDEP .GT. MD) NDEP=MD
      NDINC= 1 
C      NDINC=FLOAT(NDEP)/20.+.950001
      DO 240 J=1,NDEP
      ZR(J)=ZMIN+(J-1)*ZINC
  240 CONTINUE
C
      WRITE (6,125)
      READ(5,*) RMIN,RMAX,RINC
      WRITE(6,250) RMIN,RMAX,RINC
      WRITE(9,250) RMIN,RMAX,RINC
  250 FORMAT(1X,'RMIN= ',F10.3,5X,'RMAX= ',F10.3,5X,'RINC= ',F10.3)
C
      READ(5,*) N,IPRT,IOUT,MAMP
C
C     SET IOUT AND BSOUT
C
      IF(IOUT .GT. 0) THEN
C
C     ONE-WAY SOLUTION
C
      IOUT=1
      BSOUT=1
      ELSE IF(IOUT .LT. 0) THEN
C
C     SINGLE SCATTER APPROXIMATION
C
      IOUT=-1
      BSOUT=1
      ELSE
C
C     TWO-WAY SOLUTION
C
      IOUT=0
      BSOUT=0
      END IF
C
      WRITE (6,125)
      WRITE(6,270) N,IPRT,IOUT,MAMP
      WRITE(9,270) N,IPRT,IOUT,MAMP
  270 FORMAT(1X,'N=',I5,3X,'IPRT=',I5,3X,'IOUT=',I5,3X,
     1       'MAMP=',I5,/)
C
       IF(N .GT. NRX) THEN
       NRXX=NRX
       WRITE(9,280) N,NRXX
  280  FORMAT(1X,//,'  TOO MANY INPUT REGIONS  ',
     & /,' NUMBER OF REGIONS  : ',I5,
     & /,' MAX ALLOWED NUMBER : ',I5,
     & /,' EXECUTION IS TERMINATED ')
       STOP
       END IF 
C     WRITE(9,290) N
  290 FORMAT(1X,/,' ******',/,
     & ' THE BATHYMETRY AND SOUND VELOCITY PROFILES',
     & ' FOR THE ',I5,' INPUT REGION(S) ARE:',/,' ******',/)
C
      IMAT=0
      NRTOT=1
      IRLTOT=0
      IRETOT=0
C
      DO 500 J=1,N
C
      WRITE(6,310) J
      WRITE(9,310) J
  310 FORMAT(1X,///,'   REGION NO. ',I4,/)
C
      READ(5,*) RANG(J),IRLIN(J),DPTH(J),NPW,NPB,IREIG(J)
C
      IF(((NPW .LT. 2) .AND. (NPB .GE. 2)) .OR.
     1   ((NPW .GE. 2) .AND. (NPB .LT. 2))) THEN
      WRITE(6,*) 'SPECIFY AT LEAST TWO PROFILE POINTS IN',
     1            ' BOTH THE WATER AND BOTTOM' 
      WRITE(9,*) 'SPECIFY AT LEAST TWO PROFILE POINTS IN',
     1            ' BOTH THE WATER AND BOTTOM' 
      STOP 'INPUT BOTH WATER AND BOTTOM PROFILE'
      END IF
C
      IF( (NPW .LT. 2) .AND. (IMAT .GT. 0) ) THEN
      WRITE(6,*) 'DO NOT MIX SEARCH AND GALERKIN REGIONS' 
      WRITE(9,*) 'DO NOT MIX SEARCH AND GALERKIN REGIONS' 
      STOP 'INPUT PROFILE FOR EVERY REGION'
      END IF
C
      IRLTOT=IRLTOT+IRLIN(J)
      IRETOT=IRETOT+IREIG(J)
      IF( (IRLTOT .NE. 0) .AND. (IRETOT .NE. 0) ) THEN
      WRITE(9,*) ' IREIG OPTION NOT VALID WITH IRLIN .NE. ZERO',
     1           ' REVISE INPUT RUN STREAM '
      STOP
      END IF
C
      IF(J.LT.N)   THEN
      NRTOT=NRTOT+MAX(1,IRLIN(J))
      IF(NRTOT+N-J.GT.NRX)   THEN
      WRITE(9,*) ' TOO MANY REGIONS. REVISE INPUT RUN STREAM '
      STOP
      END IF
      END IF
      IF(NPW .LT. 2)   NPW=0
      IF(NPB .LT. 2)   NPB=0
      WRITE(LUSVP,*) NPW,NPB
      IF(NPW .GE. 2)   THEN
       IF(NPW .GT. MX) THEN
       MXP=MX
       WRITE(9,320) J,NPW,MXP
  320  FORMAT(1X,//,'  TOO MANY POINTS IN PROFILE NO. ',I5,
     & /,' NUMBER OF PROFILE POINTS : ',I5,
     & /,' MAX ALLOWED NUMBER : ',I5,
     & /,' EXECUTION IS TERMINATED ')
       STOP
       END IF 
C
C
C      ********** SET IMAT=1 **********
C
C
       IMAT=1
C
       WRITE(9,330) 
  330  FORMAT(1X,' THE SOUND VELOCITY-DENSITY PROFILE IN THE WATER',
     1        ' IS AS ','FOLLOWS:',//2X,'DEPTH (M) ',3X,
     2        'SOUND SPEED (M/SEC)',4X,
     3        'ATTENUATION (DB/WL)',4X,
     4        'DENSITY (GM/CM**3)',/)
C
       DO 340 K=1,NPW
       READ(5,*) DEPW(K),SVPW(K),DBPWLW(K),RHOW(K)
       WRITE(9,335) DEPW(K),SVPW(K),DBPWLW(K),RHOW(K)
  335  FORMAT(1X,F10.2,7X,F10.4,13X,F10.4,10X,F10.4)
  340  CONTINUE
C
       WRITE(9,125)
       WRITE(9,345) 
  345  FORMAT(1X,' THE SOUND VELOCITY-DENSITY PROFILE IN THE BOTTOM',
     1        ' IS AS ','FOLLOWS:',//2X,'DEPTH (M) ',3X,
     2        'SOUND SPEED (M/SEC)',4X,
     3        'ATTENUATION (DB/WL)',4X,
     4        'DENSITY (GM/CM**3)',/)
       DO 360 K=1,NPB
       READ(5,*) DEPB(K),SVPB(K),DBPWLB(K),RHOB(K)
       WRITE(9,335) DEPB(K),SVPB(K),DBPWLB(K),RHOB(K)
  360  CONTINUE
C
      END IF
C
      IF(IMAT .EQ. 0)   THEN
C      FIND THE WAVE NUMBER IN THE WATER
C
       FKW=2.0*PI*FREQ/CW1
       WRITE(6,365) FKW   
       WRITE(9,365) FKW   
  365  FORMAT(1X,/31H THE WAVE NUMBER IN THE WATER= ,1P,
     &        E13.5,4H 1/M/)
      ELSE
C
C      GENERATE COEFFICIENTS IN PIECEWISE LINEAR REPRESENTATION
C      OF THE SQUARE OF THE INDEX OF REFRACTION IN THE WATER AND
C      BOTTOM
C
       CALL CFACIN(FREQ,DEPW,SVPW,DBPWLW,CGRADW,CINTW,NPW,MX)
       CALL CFACIN(FREQ,DEPB,SVPB,DBPWLB,CGRADB,CINTB,NPB,MX)
C
C     FIND THE AVERAGE OF THE SQUARE OF THE INDEX OF REFRACTION
C     IN THE WATER
C
      H=DEPW(NPW)
      CALL AVGIN(DEPW,CGRADW,CINTW,NPW,MX,H,AVG,LHW)
C
C     COMPUTE THE AVERAGE SOUND SPEED AND CORRESPONDING WAVE NUMBER  
C
      CW=CB/SQRT(AVG)
      FKW=2.0*PI*FREQ/CW
C
      WRITE(9,125)
      WRITE(6,370) AVG,CW
      WRITE(9,370) AVG,CW
C
  370 FORMAT(1X,/' THE AVERAGE (INDEX OF REF.)**2 = ',E12.5,
     1       '  OR ',4PE12.5 ,' M/SEC')
C
      WRITE(6,375) FKW
      WRITE(9,375) FKW
  375 FORMAT(1X,' THE AVERAGE WAVE NUMBER IN THE WATER= ',1P,
     1        E13.5,' 1/M')
C
      WRITE(9,380)
      WRITE(6,380)
  380 FORMAT(1X,'THE PIECEWISE LINEAR (INDEX OF REF.)**2 ',
     1       'IN THE WATER IS AS FOLLOWS:',//2X,'DEPTH (M)',6X,
     2       '(INDEX OF REF.)**2',10X,'GRADIENT'/)
C
C     PRINT INDEX OF REFRACTION TABLE IN THE WATER,
C     DEFINE LOGARITHMIC GRADIENTS OF THE WATER DENSITY
C     AND WRITE RESULTS ON LUSVP FOR USE IN THE MAIN
C     PROGRAM COUPLE
C
      NPWM1=NPW-1
      DO 400 K=1,NPWM1
      GRHOW(K)=DLOG(RHOW(K+1)/RHOW(K))/(DEPW(K+1)-DEPW(K))
      WRITE(6,390) DEPW(K),CINTW(K)+CGRADW(K)*DEPW(K),CGRADW(K)
      WRITE(9,390) DEPW(K),CINTW(K)+CGRADW(K)*DEPW(K),CGRADW(K)
      WRITE(LUSVP,*) DEPW(K),CINTW(K),CGRADW(K),RHOW(K),GRHOW(K)
  390 FORMAT(1X,F10.2,2X,F10.4,E12.4,2X,2E12.4)
  400 CONTINUE
      WRITE(6,390) DEPW(NPW),CINTW(NPWM1)+CGRADW(NPWM1)*DEPW(NPW),
     1             CGRADW(NPW)
      WRITE(9,390) DEPW(NPW),CINTW(NPWM1)+CGRADW(NPWM1)*DEPW(NPW),
     1             CGRADW(NPW)
      WRITE(LUSVP,*) DEPW(NPW),CINTW(NPW),CGRADW(NPW),RHOW(NPW),
     1               GRHOW(NPW)
C
      WRITE(9,125)
      WRITE(6,125)
      WRITE(9,405)
      WRITE(6,405)
  405 FORMAT(1X,'THE PIECEWISE LINEAR (INDEX OF REF.)**2 ',
     1       'IN THE BOTTOM IS AS FOLLOWS:',//2X,'DEPTH (M)',6X,
     2       '(INDEX OF REF.)**2',10X,'GRADIENT'/)
      NPBM1=NPB-1
C
C     PRINT INDEX OF REFRACTION TABLE IN THE BOTTOM,
C     DEFINE LOGARITHMIC GRADIENTS OF THE BOTTOM DENSITY
C     AND WRITE RESULTS ON LUSVP FOR USE IN THE MAIN
C     PROGRAM COUPLE
C
      DO 420 K=1,NPBM1
      GRHOB(K)=DLOG(RHOB(K+1)/RHOB(K))/(DEPB(K+1)-DEPB(K))
      WRITE(6,390) DEPB(K),CINTB(K)+CGRADB(K)*DEPB(K),CGRADB(K)
      WRITE(9,390) DEPB(K),CINTB(K)+CGRADB(K)*DEPB(K),CGRADB(K)
      WRITE(LUSVP,*) DEPB(K),CINTB(K),CGRADB(K),RHOB(K),GRHOB(K)
  420 CONTINUE
      GRHOB(NPB)=GRHOB(NPBM1)
      WRITE(6,390) DEPB(NPB),CINTB(NPBM1)+CGRADB(NPBM1)*DEPB(NPB)
      WRITE(9,390) DEPB(NPB),CINTB(NPBM1)+CGRADB(NPBM1)*DEPB(NPB)
      WRITE(LUSVP,*) DEPB(NPB),CINTB(NPB),CGRADB(NPB),RHOB(NPB),
     1               GRHOB(NPB)
C
      END IF
C
      IF(J.EQ.1)   THEN
      IF(RANG(1).LT.1.0E-6)   THEN
      RANG(1)=1.0E-6
      WRITE(9,430)
  430 FORMAT(1X,/,' WARNING : RANGE TO END OF FIRST REGION HAS BEEN',
     & ' INCREASED TO 0.001 m ')
      END IF
      ELSE
       IF(RANG(J).LT.RANG(J-1))   THEN
       WRITE(9,*) ' ERROR : REVISE END OF REGION ',J
       STOP
      END IF
      END IF
C
      IF(DPTH(J).EQ.0.0)   THEN
      DPTH(J)=1.0E-6
      WRITE(9,*) ' WARNING : '
      WRITE(9,*)' WATER DEPTH IS INCREASED TO 1.E-6 M IN SECTOR ',J
      END IF
C
      WRITE(6,450) RANG(J),DPTH(J),IRLIN(J),IREIG(J)
      WRITE(9,450) RANG(J),DPTH(J),IRLIN(J),IREIG(J)
  450 FORMAT(1X,/,' ****** ',/,
     & ' RANGE, DEPTH OF CURRENT REGION, NUMBER OF',
     & ' SUBDIVISIONS FOR NEXT REGION AND SEQUENCE NUMBER',
     & ' FOR STORAGE OR RETRIEVAL OF EIGENVALUES: ',/,
     & 1X,F10.2,3X,F10.2,5X,I5,5X,I5,/,' ****** ',/)
C     
C     
  500 CONTINUE
C
      ENDFILE LUSVP
C
      WRITE(9,125)
C
C     GENERATE REPRESENTATION OF A SLOPING BOTTOM
C
      CALL BATHY(N,RANG,DPTH,IRLIN,RANGE,DEPTH,NREG,NRX,
     &            IMAT,NEWSVP)
C
C     EXTEND LAST REGION IF NECESSARY
C
      IF(RMAX .GT. RANGE(NREG)) RANGE(NREG)=RMAX
C
C     FIND THE COMPLEX WAVE NUMBER IN THE BOTTOM IN THE TWO
C     ISOSPEED LAYERS CASE
C
      IF(IMAT .EQ. 0) CKB=FKB*CDSRT(ONE+ATEN*CI)
C
C     PRINT RANGE,DEPTH AND THE STEEPEST DESCENT ESTIMATE OF THE
C     REFLECTED FIELD RESULTING FROM THE PRESSURE RELEASE SURFACE AT HB
C
      WRITE(6,510)
      WRITE(9,510)
  510 FORMAT(1X,'THE RANGE DEPTH PAIRS AND A STEEPEST DESCENT',
     1         ' ESTIMATE OF THE SUB-BOTTOM REFLECTED FIELD',
     2         ' (TWO ISOSPEED LAYERS CASE ONLY) ARE AS FOLLOWS:'/)
C
      NREG=MAX0(N,NREG)
      DO 600 J=1,NREG
C
      IF(IMAT .EQ. 0) THEN
C
C     COMPUTE THE SUB-BOTTOM REFLECTED FIELD IN THE TWO ISOSPEED
C     LAYER CASE. FIRST FIND COVERSION DIVISOR: CYLINDRICAL TO 
C     PLANE GEOMETRY
C
      IF( (IGEOM .EQ. 1) .AND. (J .EQ. 1) ) THEN
      CPG=CDSRT(CI*PI*CKB)*DCMPLX(BESSJ0(FKW),BESSY0(FKW))
      END IF
C
      H=DEPTH(J)
      RNG=1000.*RANGE(J)
      R=SQRT(RNG**2+4.*(HB-H)**2)
      FKH=CKB*RNG/R
      FKVS=CKB*2.*(HB-H)/R
      FKVW=CDSRT(FKW**2-FKH**2)
      IF(DIMAG(-CKB*R) .LT. -88.) GO TO 520
      CXP=CDEXP(CI*CKB*R)
      CDEM=FKVW*ROHB*CDCOS(FKVW*H)-CI*FKVS*ROHW*CDSIN(FKVW*H)
      PREFL=CXP/R
      PREFL=PREFL*4.0*ROHB*ROHW*(FKVS*CDSIN(FKVW*H)/CDEM)**2
C
C     ADJUST FOR PLANE GEOMETRY
C
      IF(IGEOM .EQ. 1) PREFL=SQRT(2.0*R)*PREFL/CPG
C
      IF(CDABS(PREFL)**2 .EQ. 0.0) GO TO 520
      TLREFL=-20.*LOG10(CDABS(PREFL))
C
      GO TO 530
  520 CONTINUE
      TLREFL=999.99D0
  530 CONTINUE
      ELSE
C
C     THE ESTIMATE OF THE SUBOTTOM REFLECTED FIELD IS ONLY
C     VALID FOR TWO ISOSPEED LAYERS. IT IS SET TO ZERO
C     IN ALL OTHER CASES.
C
      TLREFL=0.0D0
      END IF
C
      WRITE(6,550) J,RANGE(J),DEPTH(J),TLREFL
      WRITE(9,550) J,RANGE(J),DEPTH(J),TLREFL
  550 FORMAT(1X,I5,5X,F10.3,5X,F10.2,5X,F10.2)
C
  600 CONTINUE
C
C
      RETURN
      END
