	PROGRAM		WKBZ

	IMPLICIT	REAL*8(A-H, O-Z)

	PARAMETER	( NIU = 1, NPU = 6, NOU = 2,
     &		    	MAXLYR = 100,  MAXLYR1 = 101, MAXSSP = 10, 
     &		    	MAXMODE = 2000,  MAXNRC = 1000, 
     &		    	MAXSOU = 2, MAXREC = 4 ) 

	REAL*4		RANGE00(MAXSSP), DEPTH00(MAXSSP), 
     &			Z00(MAXSSP, MAXLYR1), C00(MAXSSP, MAXLYR1)
	INTEGER		NLAYER00(MAXSSP)

	REAL*8		ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR),
     &			DEPSOU(MAXSOU), DEPREC(MAXREC), RANREC(MAXNRC), 
     &			EIGVL1(MAXMODE), EIGVL2(MAXMODE), UNRMT(MAXREC),
     &			UNRM1(MAXMODE,MAXSOU), UNRM2(MAXMODE,MAXREC)
	COMPLEX*16	YFACTOR, YY(MAXSOU,MAXREC,MAXNRC)
	REAL*4		RTL, YYTL, TLS(MAXSOU, MAXREC)

	CHARACTER*80	TITLE

	COMMON		NLAYER
	COMMON /EIGEN/  PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS,
     &			WFREQ, WFREQ2, WKZ0, WKSUR2 	


	PI   = DATAN( 1.D0 ) * 4.
	PREC = 1.D-8

CCCCCC	PART 1	CCCCCC

	OPEN(UNIT=NIU, FILE='WKBZ.IN', FORM='FORMATTED', STATUS='OLD')
	OPEN(UNIT=NOU, FILE='WKBZ.OUT', FORM='FORMATTED', STATUS='NEW')

	READ(NIU, *)	TITLE
	WRITE(NPU, *)	TITLE

	READ(NIU, *)	FREQ, RMKM, DRKM

	WFREQ  = FREQ * PI * 2.
	WFREQ2 =  WFREQ * WFREQ
	NRCALC = RMKM / DRKM + 1./2.
	IF ( NRCALC .GT. MAXNRC ) THEN
	  WRITE(NPU, 1002) MAXNRC, NRCALC
	  STOP
	ENDIF
	RMKM   = DRKM * NRCALC
	DR     = DRKM * 1000.
	DO 105 I = 1, NRCALC
	RANREC(I) = DR * I
105	CONTINUE
	WRITE(NPU,1000)	FREQ
	WRITE(NPU,1001) RMKM
	READ(NIU, *)	NSOU
	IF ( NSOU .GT. MAXSOU ) THEN
	  WRITE(NPU, 1003) MAXSOU, NSOU
	  STOP
	ENDIF
	READ(NIU, *)	(DEPSOU(I), I=1,NSOU)

	READ(NIU, *)	NREC
	IF ( NREC .GT. MAXREC ) THEN
	  WRITE(NPU, 1004) MAXREC, NREC
	  STOP
	ENDIF
	READ(NIU, *)	(DEPREC(I), I=1,NREC)

1000	FORMAT( 4X, 'Frequency = ', G11.4, 'Hz', / )
1001	FORMAT( 4X, 'MAX. RANGE = ', G11.4, 'km', // )
1002	FORMAT( 4X, 'MAXNRC =', I5, 4X, 'NRCALC = ', I5, / )
1003	FORMAT( 4X, 'MAXSOU =', I5, 4X, 'NSOU = ', I5, / )
1004	FORMAT( 4X, 'MAXREC =', I5, 4X, 'NREC = ', I5, / )

CCCCCC	PART 2	CCCCCC

	NSSP = 0
110	CONTINUE
	READ(NIU, *, END=120)	RANGE1, DPH11
	WRITE(NPU, 1010)	RANGE1, DPH11
1010	FORMAT(8X, 'SSP RANGE =',G11.4, 'km, ', 4X, 'DEPTH =',G11.4, 'm',/)

	NSSP = NSSP + 1
	RANGE00(NSSP) = RANGE1 * 1000.
	DEPTH00(NSSP) = DPH11

	DO 115 I = 1, MAXLYR1
	READ(NIU, *)  ZTMP, CTMP
	IF ( DABS( ZTMP - DPH11 ) .LT. PREC ) THEN
	    Z00(NSSP,I) = DPH11
	    C00(NSSP,I) = CTMP
	    NLAYER00(NSSP) = I - 1
	    GOTO 116
	  ELSE
	    Z00(NSSP,I) = ZTMP
	    C00(NSSP,I) = CTMP
	ENDIF
115	CONTINUE
116	CONTINUE

	IF( RANGE1 .LT. RMKM )  GOTO 110

120	CONTINUE
	CLOSE(NIU)


CCCCCC	PART 3	CCCCCC

1020	FORMAT( /4X, 'ISSP = ', I4, 4X, 'RANGE = ', G11.4, 'km', // ) 

	NMODE = MAXMODE
	ISSP  = 1	

	RANGE1 = 0.
	RANGE2 = 0.
	WRITE(NPU,1020) ISSP, RANGE2/1000.

	DPH11  = DEPTH00(ISSP)
	NLAYER = NLAYER00(ISSP)
	NPLUS1 = NLAYER + 1

	DO 125 I = 1, NPLUS1	    
	ZI( I ) = Z00(ISSP,I)
	CI0(I)  = C00(ISSP,I)
125	CONTINUE

	LZ0AXIS = 1
	DO 129 I = 1, NLAYER
	  I1    = I + 1
	  IF( CI0(I1) .EQ. CI0(I) ) CI0(I1) = CI0(I1) - 0.01
	  AI(I) = ( CI0(I1) - CI0(I) )/( ZI(I1) - ZI(I) )
	  IF ( CI0(I1) .LT. CI0(LZ0AXIS) ) LZ0AXIS = I1
129	CONTINUE
	Z0AXIS  = ZI(LZ0AXIS)			!!!! CHANNEL AXIS 
	CZ0AXIS = CI0(LZ0AXIS)

	WKZ0    = WFREQ / CZ0AXIS
	WKSUR2  = WFREQ2 / ( CI0(1) * CI0(1) )
	IF( LZ0AXIS .GT. 1 )  THEN
	    ULIN = WFREQ/CI0( LZ0AXIS-1 )
	  ELSE
	    ULIN = WFREQ/CI0( CI0(1) + CI0(2) ) * 2.
	ENDIF

	LTHVAL = 0 
130	CONTINUE

	LTHVAL1 = LTHVAL + 1

	CALL	EIGENVAL( LTHVAL, ULIN, ULOUT, XSP, PHADJ,
     &			ZREA, LREA, ZREB, LREB, MMSUR, MMBOT )
	IF ( MMBOT .NE. 0 ) THEN
	    NMODE = LTHVAL
	    GOTO 140
	ENDIF
	EIGVL1(LTHVAL1) = 0.
	EIGVL2(LTHVAL1) = ULOUT
	write(*,*) lthval1, ulout, wfreq/ulout

	DO 135 I = 1, NSOU
	ZDPH = DEPSOU(I)
	CALL	EIGENFUN( LTHVAL, ULOUT, XSP, PHADJ,
     &  			ZREA, LREA, ZREB, LREB, ZDPH, YL )
	UNRM1( LTHVAL1, I) = YL
135	CONTINUE
	DO 139 I = 1, NREC
	ZDPH = DEPREC(I)
	CALL	EIGENFUN( LTHVAL, ULOUT, XSP, PHADJ,
     &				ZREA, LREA, ZREB, LREB, ZDPH, YL )
	UNRM2( LTHVAL1, I) = YL
139	CONTINUE
	
	LTHVAL = LTHVAL1
	IF( LTHVAL .LT. NMODE ) THEN
	  ULIN = ULOUT - PI * 2. /XSP
	  GOTO 130
	ENDIF

CCCCCC	PART 4	CCCCCC

140	CONTINUE
	NEEDBG = 1
	NEEDLF = NRCALC - NEEDBG + 1

150	CONTINUE
	IF ( ISSP .EQ. NSSP ) GOTO 200

CCCCCC	PART 5	CCCCCC

	ISSP = ISSP + 1

	RANGE1 = RANGE2
	RANGE2 = RANGE00(ISSP)
	WRITE(NPU,1020) ISSP, RANGE2/1000.

	DPH11  = DEPTH00(ISSP)
	NLAYER = NLAYER00(ISSP)
	NPLUS1 = NLAYER + 1

	INR    = NEEDBG
	NEEDCA = 0
161	CONTINUE
	IF ( RANREC( INR ) .GT. RANGE2 ) GOTO 162
	INR    = INR + 1
	NEEDCA = NEEDCA + 1
	IF ( INR .LE. NRCALC ) THEN
	    GOTO 161
	  ELSE
	    NEEDLF = 0
	    GOTO 163
	ENDIF
162	NEEDLF = NRCALC - INR + 1
163	CONTINUE

	DO 165 I = 1, NPLUS1	    
	ZI( I ) = Z00(ISSP,I)
	CI0(I)  = C00(ISSP,I)
165	CONTINUE

	LZ0AXIS = 1
	DO 169 I = 1, NLAYER
	  I1    = I + 1
	  IF( CI0(I1) .EQ. CI0(I) ) CI0(I1) = CI0(I1) - 0.01
	  AI(I) = ( CI0(I1) - CI0(I) )/( ZI(I1) - ZI(I) )
	  IF ( CI0(I1) .LT. CI0(LZ0AXIS) ) LZ0AXIS = I1
169	CONTINUE
	Z0AXIS  = ZI(LZ0AXIS)			!!!! CHANNEL AXIS 
	CZ0AXIS = CI0(LZ0AXIS)

	WKZ0 = WFREQ / CZ0AXIS
	IF( LZ0AXIS .GT. 1 )  THEN
	    ULIN = WFREQ/CI0( LZ0AXIS-1 )
	  ELSE
	    ULIN = WFREQ/CI0( CI0(1) + CI0(2) ) * 2.
	ENDIF

	LTHVAL = 0 
170	CONTINUE

	LTHVAL1 = LTHVAL + 1

	CALL	EIGENVAL( LTHVAL, ULIN, ULOUT, XSP, PHADJ,
     &				ZREA, LREA, ZREB, LREB, MMSUR, MMBOT )
	IF ( MMBOT .NE. 0 ) THEN
	    NMODE = LTHVAL
	    GOTO 190
	ENDIF
	write(*,*) lthval1, ulout, wfreq/ulout

	DO 172 I = 1, NREC
	ZDPH = DEPREC(I)
	CALL	EIGENFUN( LTHVAL, ULOUT, XSP, PHADJ,
     &				ZREA, LREA, ZREB, LREB, ZDPH, YL )
	UNRMT(I) = YL
172	CONTINUE

	IF ( NEEDCA .EQ. 0 ) GOTO 180 
	DO 179 INR = 1, NEEDCA
	INRR    = INR + NEEDBG - 1
	RANRR   = RANREC(INRR) - RANGE1
	RANRR1  = RANRR / ( RANGE2 - RANGE1 )
	ULOUTR  = EIGVL2(LTHVAL1) + ( ULOUT - EIGVL2(LTHVAL1) ) * RANRR1
	WNPHA   = EIGVL1(LTHVAL1) + (EIGVL2(LTHVAL1)+ULOUTR) * RANRR / 2.
	YFACTOR = CDEXP( DCMPLX(0.D0, WNPHA) ) / DSQRT( WNPHA )
	DO 175 IR = 1, NREC
	YLR     = UNRM2(LTHVAL1, IR) + 
     &			( UNRMT(IR)-UNRM2(LTHVAL1, IR) )*RANRR1
	DO 174 IS = 1, NSOU
	YRINRR  = UNRM1(LTHVAL1, IS) * YLR 
	YY(IS, IR, INRR)  = YY(IS, IR, INRR) + YFACTOR * YRINRR
174	CONTINUE	
175	CONTINUE
179	CONTINUE

180	CONTINUE		
	EIGVL1(LTHVAL1) = EIGVL1(LTHVAL1) + 
     &			( EIGVL2(LTHVAL1) + ULOUT ) * (RANGE2-RANGE1) / 2.
	EIGVL2(LTHVAL1) = ULOUT
	DO 182 I = 1, NREC
	UNRM2(LTHVAL1, I) = UNRMT(I)
182	CONTINUE

	
	LTHVAL = LTHVAL1
	IF( LTHVAL .LT. NMODE ) THEN
	  ULIN = ULOUT - PI * 2. /XSP
	  GOTO 170
	ENDIF

190	CONTINUE
	IF ( NEEDLF .EQ. 0 ) THEN
	    GOTO 300
	  ELSE
	    NEEDBG = NRCALC - NEEDLF + 1
	    GOTO 150
	ENDIF
	
CCCCCC	PART 6	CCCCCC

200	CONTINUE
	RANGE1 = RANGE2
	DO 210 INRR = NEEDBG, NRCALC
	RANRR   = RANREC( INRR ) - RANGE1
	DO 209 LTHVAL = 1, NMODE
	WNPHA   = EIGVL1(LTHVAL) + EIGVL2(LTHVAL) * RANRR
	YFACTOR = CDEXP( DCMPLX(0.D0, WNPHA) ) / DSQRT( WNPHA )
	DO 205 IR = 1, NREC
	YLR     = UNRM2(LTHVAL, IR)
	DO 204 IS = 1, NSOU
	YRINRR  = UNRM1(LTHVAL, IS) * YLR 
	YY(IS, IR, INRR)  = YY(IS, IR, INRR) + YFACTOR * YRINRR
204	CONTINUE	
205	CONTINUE
209	CONTINUE
210	CONTINUE
	
CCCCCC	PART 7	CCCCCC

300	CONTINUE

	WRITE(NOU,*) TITLE
	WRITE(NOU, 1000) FREQ
	WRITE(NOU, 2000) NSOU
	WRITE(NOU, 2001) (DEPSOU(IS), IS=1, NSOU)
	WRITE(NOU, 2002) NREC
	WRITE(NOU, 2003) (DEPREC(IR), IR=1, NREC)
2000	FORMAT(/4X, 'SOURCE', 6X, I5, 6x, 'DEPTH IN METERS')
2001	FORMAT(/4X, 4G11.4)
2002	FORMAT(/4X, 'RECEIVER', 4X, I5, 6x, 'DEPTH IN METERS')
2003	FORMAT(/4X, 4G11.4)
	WRITE(NOU, 2004)
2004	FORMAT(//X, 'RANGE (km)   &   TRANSMISSION LOSS (dB)', //)

	FACTOR  = 20. / LOG(10.)
	YFACTOR = CDEXP( DCMPLX(0.D0, PI/4.) ) * DSQRT( PI * 2. )
	IF ( (NSOU + NREC) .GT. 2 ) GOTO 320

	DO 310 INR = 1, NRCALC
	RTL = RANREC(INR) / 1000.
	YY( 1, 1, INR) = YFACTOR * YY( 1, 1, INR )
	TLS( 1, 1 )    = FACTOR * LOG( CDABS( YY( 1, 1, INR ) ) )
	YYTL           = CDABS( YY( 1, 1, INR ) )
	IF ( YYTL .LE. 1.D-10 ) THEN
	    TLS( 1, 1 ) = -200.
	  ELSE
	    TLS( 1, 1 ) = FACTOR * LOG( YYTL )
	ENDIF
	WRITE(NOU, *) RTL, TLS(1,1)
310	CONTINUE
	GOTO 350
	
320	CONTINUE

	DO 329 INR = 1, NRCALC
	RTL = RANREC(INR) / 1000.
	DO 329 IS = 1, NSOU
	DO 328 IR = 1, NREC
	YY( IS, IR, INR) = YFACTOR * YY( IS, IR, INR )
	YYTL             = CDABS( YY( IS, IR, INR ) )
	IF ( YYTL .LE. 1.D-10 ) THEN
	    TLS( IS, IR ) = -200.
	  ELSE
	    TLS( IS, IR ) = FACTOR * LOG( YYTL )
	ENDIF
328	CONTINUE
	WRITE(NOU,*) RTL, ( TLS(IS,IR), IR=1, NREC )	
329	CONTINUE

350	CLOSE(NOU)

	STOP
	END
	

	SUBROUTINE	EIGENVAL( LTHVAL, ULIN, ULOUT, XSP, PHADJSUR,
     &				ZREA, LREA, ZREB, LREB, MMSUR, MMBOT)

	IMPLICIT	REAL*8(A-H, O-Z)

	PARAMETER	( MAXLYR = 100,  MAXLYR1 = 101, PREC = 1.D-8 ) 

	REAL*8		ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR)

	COMMON		NLAYER
	COMMON /EIGEN/  PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS,
     &			WFREQ, WFREQ2, WKZ0, WKSUR2 	


	NITER  = 0
	UL0    = ULIN

110	NITER  = NITER + 1
	CALL	SPANPHASE( UL0,XSP,XPH,ZREA,LREA,ZREB,LREB,MMSUR,MMBOT )

	TVARAI0 = WKSUR2 - UL0 * UL0
	IF ( MMSUR .NE. 0 ) THEN
	    OMEGA0 = DABS( TVARAI0 )
	    OMEGA0 = -( CI0(1)*( OMEGA0**(3./2.) ) )/(AI(1) * WKSUR2 * 3.)
	  ELSE
	    OMEGA0 = 0.
	ENDIF
	DGRAD23 = ( DABS( -AI(1) * WKSUR2 * 2. / CI0(1) ) ) ** (2./3.)
	TVARAI0 = - TVARAI0 / DGRAD23
	CALL	AIRYFUNC( TVARAI0, THETA )

	PHADJSUR  = PI / 2. + DATAN( DTAN( ( THETA - OMEGA0 ) * 2. ) )

	UL1  = UL0 - ( ( PI*LTHVAL*2. + PI/2. + PHADJSUR ) - XPH ) / XSP
	IF ( DABS( UL1 - UL0 )/UL1 .GT. PREC ) THEN
	    UL0 = UL1
	    GOTO 110
	ENDIF

	ULOUT = UL1

	RETURN
	END




	SUBROUTINE	EIGENFUN( LTHVAL, UL, XSP, PHADJSUR,
     &					ZREA, LREA, ZREB, LREB, ZD, YZ)

	IMPLICIT	REAL*8(A-H, O-Z)

	PARAMETER	( MAXLYR = 100,  MAXLYR1 = 101, PREC = 1.D-8 ) 

	REAL*8		ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR)

	COMMON		NLAYER
	COMMON /EIGEN/  PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS,
     &			WFREQ, WFREQ2, WKZ0, WKSUR2 	


	PARAMETER	( BPARA = 2.152D0, DPARA = 1.619D0 )


C
	LZD    = 0	
110	LZD    = LZD + 1
	IF ( ZD .GT. ZI( LZD + 1 ) )  GOTO 110
	CZD    = CI0( LZD ) + AI( LZD ) * ( ZD - ZI( LZD ) )
C
	UKZ     = WFREQ / CZD
	SQUKZ   = UKZ * UKZ
	EKSQZ   = DABS( -SQUKZ * AI( LZD ) * 2. / CZD )
	EKS23   = EKSQZ ** ( 2./3. )
	EKS43   = EKSQZ ** ( 4./3. )
 	SQK     = ( UKZ + UL ) * ( UKZ - UL )
	SQK2    = SQK * SQK
C
	IF ( ZD .GE. ZREB ) GOTO 400
	IF ( ZD .GE. ZREA ) GOTO 300

CCCCCC	[ 0, ZREA )	CCCCCC

	CONST  = WFREQ/UL

	XCOSA	= CONST / CZD
	XTANA	= DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA
	XANGL   = DATAN( XTANA )
	XXP0	= XTANA - XANGL

	XPHASE	= 0.
	IF ( LZD .GE. LREA ) GOTO 210
	DO II = LZD+1, LREA
	  XCOSA  = CONST / CI0( II )
	  XTANA  = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA
	  XANGL  = DATAN( XTANA )
	  XXPP	 = XTANA - XANGL
	  XPHASE = XPHASE + ( XXPP - XXP0 ) / AI( II-1 )
	  XXP0   = XXPP
	END DO

210	XPHASE = XPHASE - XXP0 / AI ( LREA )
	XPHASE = WFREQ * XPHASE
C
	FRACT  = XPHASE 
	FRACT  = DEXP( -FRACT ) * DCOS( PHADJSUR / 2. )
	DENOM  = BPARA * EKS43 - DPARA * EKS23 * SQK + SQK2 * 16.
	GOTO 500

CCCCCC	[ ZREA, ZREB )	CCCCCC

300	CONTINUE

	CONST = UL / WFREQ
	
	IF (ZREA .EQ. 0. ) THEN
	    XSINA  = CI0(1) * CONST
	    XSINA  = DSQRT( DABS( 1.D0 - XSINA*XSINA ) )
	  ELSE
	    XSINA  = 0.
	END IF
	XXP0 = DLOG ( ( 1.D0 + XSINA )/( 1.D0 - XSINA ) ) - XSINA * 2.

	XPHASE = 0
	IF ( LREA .GE. LZD ) GOTO 310
	DO II = LREA+1, LZD
	  XSINA  = CI0( II ) * CONST
	  XSINA  = DSQRT( DABS( 1.D0 - XSINA*XSINA ) )
	  XXPP   = DLOG( ( 1.D0 + XSINA )/( 1.D0 - XSINA ) ) - XSINA * 2.
	  XPHASE = XPHASE + ( XXP0-XXPP ) / AI( II-1 )
	  XXP0   = XXPP
	END DO

310	XSINA  = CZD * CONST
	XSINA  = DSQRT( DABS( 1.D0 - XSINA * XSINA ) )
	XXPP   = DLOG( ( 1.D0 + XSINA )/( 1.D0 - XSINA ) ) - XSINA * 2.

	XPHASE = XPHASE + ( XXP0-XXPP ) / AI(LZD)
	XPHASE = WFREQ * XPHASE / 2.
C
	FRACT  = XPHASE + ( PI - PHADJSUR ) / 2.
	FRACT  = DSIN( FRACT )
	DENOM  = BPARA * EKS43 - DPARA * EKS23 * SQK + SQK2

	GOTO 500

CCCCCC [ZREB, DPH11 ]	CCCCCC

400	CONTINUE
	
	CONST = WFREQ/UL

	IF ( ZREB .EQ. DPH11 ) THEN
	    XCOSA  = CONST / CI0( NLAYER + 1 )
	    XTANA  = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA
	    XANGL  = DATAN( XTANA )
	  ELSE
	    XTANA=0.
	    XANGL=0.
	ENDIF
	XXP0  = XTANA - XANGL
	

	XPHASE = 0
	IF ( LREB .GE. LZD ) GOTO 410
	DO II = LREB+1, LZD
	  XCOSA  = CONST / CI0( II )
	  XTANA  = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA
	  XANGL  = DATAN( XTANA )
	  XXPP   = XTANA - XANGL
	  XPHASE = XPHASE + DABS( ( XXPP - XXP0 ) / AI( II-1 ) )
	  XXP0   = XXPP
	END DO

410	XCOSA  = CONST / CZD
	XTANA  = DSQRT( DABS( 1.D0 - XCOSA * XCOSA ) ) / XCOSA
	XANGL  = DATAN( XTANA )
	XXPP   = XTANA - XANGL

	XPHASE = XPHASE + DABS( ( XXPP - XXP0 ) / AI( LZD ) )
	XPHASE = WFREQ * XPHASE
C
	FRACT  = DFLOAT((-1)**LTHVAL) * DEXP( -XPHASE ) / DSQRT( 2.D0 )
	DENOM  = BPARA * EKS43 - DPARA * EKS23 * SQK + SQK2 * 16.

CCCCCC

500	FRACTOR  = DSQRT( UL * 4. / XSP )
	YZ       = FRACTOR * FRACT / DENOM ** ( 1./8. )

	RETURN
	END


	SUBROUTINE	SPANPHASE( UL, XSPAN, XPHASE, ZREA, LREA,
     &					ZREB, LREB, MMSUR, MMBOT )


	IMPLICIT	REAL*8(A-H, O-Z)

	PARAMETER	( MAXLYR = 100,  MAXLYR1 = 101 ) 

	REAL*8		ZI(MAXLYR1), CI0(MAXLYR1), AI(MAXLYR)
	REAL*8		XSINA(MAXLYR1)

	COMMON		NLAYER
	COMMON /EIGEN/  PI,DPH11,ZI,CI0,AI,Z0AXIS,CZ0AXIS,LZ0AXIS,
     &			WFREQ, WFREQ2, WKZ0, WKSUR2 	

C
	MMSUR=0
	MMBOT=0
	
C	CALCULATE THE SOUND SPEED AT THE TURNING-POINT DEPTH

	CTURN  = WFREQ/UL

C	ZREA -- THE RAY TURNING-DEPTH ABOVE THE CHANNEL AXIS
	LREA = LZ0AXIS - 1

110	CONTINUE
	IF ( LREA .EQ. 0 ) THEN
	    MMSUR = 1
	    ZREA  = 0.
	    GOTO 120
	ENDIF	
	IF ( CTURN .LT. CI0( LREA ) ) THEN	
	    GOTO 111
	  ELSE
	    LREA  = LREA - 1	
	    GOTO 110
	ENDIF
111	ZREA  = ZI( LREA ) + ( CTURN - CI0( LREA ) ) / AI( LREA )
120	CONTINUE

C	ZREB -- THE RAY TURNING-DEPTH ABOVE THE CHANNEL AXIS
	LREB = LZ0AXIS

130	CONTINUE
	IF ( LREB .GT. NLAYER ) THEN
	    MMBOT = 1
	    ZREB = DPH
	    GOTO 140
	ENDIF
	IF ( CTURN .LT. CI0( LREB+1 ) ) THEN
	    GOTO 131
	  ELSE
	    LREB  = LREB + 1
	    GOTO 130
	ENDIF
131	ZREB  = ZI( LREB ) + ( CTURN - CI0( LREB ) ) / AI( LREB )
140	CONTINUE


C	SOUND RAY SPAN

	CONST = UL / WFREQ

	XSPAN  = 0.
	XPHASE = 0.
	
	IF ( LREA .EQ. 0) THEN
	    LREA     = 1
	    XXSINA   = CI0(1) * CONST
	    XSINA(1) = DSQRT( DABS( 1.D0 - XXSINA * XXSINA ) )
	  ELSE
	    XSINA(LREA) = 0.
	ENDIF
	IF ( LREB .GT. NLAYER ) THEN
	    LREB            = NLAYER
	    XXSINA          = CI0( NLAYER+1 ) * CONST
	    XSINA(NLAYER+1) = DSQRT( DABS( 1.D0 - XXSINA * XXSINA ) )
	  ELSE
	    XSINA( LREB+1 ) = 0.
	END IF

	XXP0= DLOG( ( 1.D0 + XSINA(LREA) ) / ( 1.D0 - XSINA(LREA) ) )
     1			- XSINA( LREA ) * 2.


	DO II = LREA+1, LREB
	  XXSINA    = CI0( II ) * CONST
	  XSINA(II) = DSQRT( DABS( 1.D0 - XXSINA * XXSINA ) )
	  XXPP      = DLOG( ( 1.D0 + XSINA(II) ) / ( 1.D0 - XSINA(II) ) )
     &			- XSINA( II ) * 2.
	  XSPAN     = XSPAN + ( XSINA( II-1 ) - XSINA( II ) ) / AI( II-1 )
	  XPHASE    = XPHASE + ( XXP0 - XXPP ) / AI( II-1 )
	  XXP0      = XXPP 
	END DO

	XSPAN  = XSPAN + ( XSINA( LREB ) - XSINA( LREB+1 ) ) / AI( LREB )
	XSPAN  = XSPAN * 2. / CONST

	XXPP   = DLOG( ( 1.D0 + XSINA(LREB+1) )/( 1.D0 - XSINA(LREB+1) ) )
     1			- XSINA(LREB+1) * 2.
	XPHASE = XPHASE + ( XXP0 - XXPP ) / AI( LREB )
	XPHASE = WFREQ * XPHASE


	RETURN
	END


	SUBROUTINE	AIRYFUNC( XTIN, THETA )

	IMPLICIT	REAL*8(A-H, O-Z)
	PARAMETER	( PREC = 1.D-6 )

	IF ( XTIN .EQ. 0.D0 ) THEN
	  UAIX=1.0899290710			!!XTIN = 0.
	  VAIX=0.6292708425
	  GOTO 2500
	ENDIF
C
C
	IF ( ( XTIN .LT. 5. ) .AND. ( XTIN .GT. -8. ) ) GOTO 2460

	IF ( XTIN .LT. 0. ) GOTO 2450			
C
C
	XW   = ( XTIN**(3./2.) )*2./3.			!! XTIN > 0.
	EXW  = DEXP( -XW * 2. )
	XW72 = 1.D0 / (XW*72.)

	NAIITER=1
	XXAI=DFLOAT(5)*XW72
	UAIX=DFLOAT(1)+XXAI
	VAIX=DFLOAT(1)-XXAI
	
2401	NAIITER  = NAIITER+1
	NAIITER6 = 6*NAIITER
	XXAI     = XXAI*XW72*(NAIITER6-1)*(NAIITER6-5)/NAIITER
	UAIX     = UAIX+XXAI
	VAIX     = VAIX+XXAI*((-1)**NAIITER)
	IF ( DABS(XXAI/UAIX) .GE. PREC ) GOTO 2401

	THETA  = DATAN( VAIX/UAIX * EXW * 0.5 )

	GOTO 2500
C
C
2450	XTIN  = DABS(XTIN)				!! XTIN < 0.
	XW    = (XTIN**(3./2.))*2./3.
	XWSS  = 1.D0 / ( XTIN**(1./4.) )
	PI4   = DATAN( 1.D0 )
	XWCOS = DCOS( XW+PI4 )
	XWSIN = DSIN( XW+PI4 )
	XW72  = 1.D0 / ( XW*72. )

	NAIITER = 1
	XXUAI   = 1.D0
	XXVAI   = XW72 * 5.
	UAIX    = 1.D0
	VAIX    = XXVAI

2451	NAIITER  = NAIITER+1
	NAIITER6 = 6*NAIITER
	XXUAI    = -XXVAI*XW72*(NAIITER6-1)*(NAIITER6-5)/NAIITER
	NAIITER  = NAIITER+1
	NAIITER6 = 6*NAIITER
	XXVAI    = XXUAI*XW72*(NAIITER6-1)*(NAIITER6-5)/NAIITER
	UAIX     = UAIX + XXUAI
	VAIX     = VAIX + XXVAI

	IF((DABS(XXUAI/UAIX).GE. PREC ).OR.(DABS(XXVAI/VAIX).GE. PREC )) 
     1			GOTO 2451
	XXUAI  = XWCOS*UAIX + XWSIN*VAIX
	VAIX   = XWSIN*UAIX - XWCOS*VAIX
	UAIX   = XXUAI

	THETA  = DATAN( VAIX/UAIX )

	GOTO 2500
C
2460	UAIX0  = 1.0899290710			!!!!
	VAIX0  = 0.6292708425
	UAIX10 = 0.7945704238
	VAIX10 = -0.4587454481
	XTIN3  = XTIN ** 3.

	XXUAI  = 1.D0
	XXVAI  = 1.D0
	UAIX   = 1.D0
	VAIX   = 1.D0

	NAIITER = 0

2461	NAIITER = NAIITER + 1
	NAI3    = NAIITER * 3
	XXUAI   = XXUAI * XTIN3 / ( ( NAI3-1 ) * NAI3 )		
	XXVAI   = XXVAI * XTIN3 / ( NAI3 * ( NAI3+1 ) )
	UAIX    = UAIX + XXUAI
	VAIX    = VAIX + XXVAI

	IF ( DABS( XXVAI/VAIX ).GT. PREC ) GOTO 2461
	XXUAI   = UAIX0 * UAIX + UAIX10 * VAIX * XTIN
	VAIX    = VAIX0 * UAIX + VAIX10 * VAIX * XTIN
	UAIX    = XXUAI

	THETA  = DATAN( VAIX/UAIX )

2500	CONTINUE

	RETURN
	END


