	subroutine qpar grad
c	a single quasi-parabolic layer with latitudinal gradients in
c	FC and HMAX, used in JONES raytracing
c	written by leo mcnamara 15-may-87
c-----	last modified 23-mar-88
c
	implicit real*8 (a-h,o-z)
	character*8 modx,id
	common r(6) /ww/ id(10),w0,w(400)
	common /xx/ modx(2),x,pxpr,pxpth,pxpph,pxpt,hmax
	equivalence (r0,w(2)), (f,w(6)), (a,w(300)), (b,w(301)),
     a  (th00,w(302)), (fc0,w(303)), (hm0,w(304)), (ym,w(305)),
     a  (theta,r(2)), (rr,r(1)), (fbeg,w(7)), (fend,w(8)),
     a  (fstep,w(9)), (elbeg,w(15)), (elend,w(16)), (elstep,w(17))

	data modx / 'qpar-', 'grad' /
	data ipass / 0 /
c
c................................
	entry ELECTX
c................................
c
	fn = 0.
	x = 0.
	pxpr = 0.
	pxpth = 0.
	pxpph = 0.
	pxpt = 0.
c
c	convert starting latitude to co-latitude
	pi = 3.141592654
	th0 = pi/2. - th00
c
c	check for point below base of layer
	hm = hm0 * (1.+b*(theta - th0))
	rm = r0 + hm
	rb = rm - ym
	if(rr.lt.rb) return
c
c	check for point above the layer
	if(rr.gt.rm+ym) return
c
c	critical frequency of layer
	fc = fc0 * (1.+a*(theta - th0))
c
c	HMAX to be returned to RAYTRACE
	hmax = hm
c
c	the parameter X = (FN/F)**2
	z = (rm - rr) * rb / ym / rr
	fn2 = fc * fc * (1. - z * z)
	fn = dsqrt (fn2)
	x = fn2 / f / f
cc	if(x.gt.1.0) type*,' rr,fn,f,x',rr,fn,f,x
	if(x.gt.1.0) return
c
c	partial derivative of X wrt r
	pxpfn = 2. * fn / f / f
	pfnpz =  - fc * fc * z / fn
	pzpr = - rb * rm / ym / rr / rr
	pxpr = pxpfn * pfnpz * pzpr
c
c	partial derivative of X wrt theta (latitude)
	pfnpfc = dsqrt (1.-z*z)
	pfcpth = a * fc0
	pfnpz = - fc * z / dsqrt(1.- z*z)
c	make sure we are dealing with the bottom of the layer
	pzphm = (rm - rr + rb) / rr / ym
	phmpth = b * hm0
	pfnpth = pfnpfc * pfcpth + pfnpz * pzphm * phmpth
	pxpth = pxpfn * pfnpth 
c
	h = rr - r0
cc	type*,'h theta hm fc x pxpr pxpth',h,theta,hm,fc,x,pxpr,pxpth
cc	type*,' z,fn,tem ',z,fn,tem
cc	type*,' th00,th0,theta,hm,fc ',th00,th0,theta,hm,fc
c
c	if the gradients = 0, profile defaults to the normal qp - calculate
c	the range analytically - see Croft & Hoogasian Radio Science
c	3(1), 69-74, 1968.  --  do on first pass --
	if(a.eq.0..and.b.eq.0.) go to 20
	return
20	if(ipass.eq.1) return
	ipass = 1
	write(2,99)
99	format(//5x,'Single q-p layer results '//)
	nfreq = 1
	if(fstep.ne.0.) nfreq = (fend-fbeg)/fstep + 1.5
	nbeta = 1
	if(elstep.ne.0.) nbeta = (elend-elbeg)/elstep + 1.5
	do 50 nf = 1,nfreq
	fr = fbeg + (nf-1)*fstep
	do 40 i = 1,nbeta
	beta = elbeg + (i-1) * elstep
	ff = fr / fc
	cosbz = dcos (beta)
	aa = 1. - 1./ff/ff + (rb/ff/ym)**2
	bb = - 2.*rm*(rb/ff/ym)**2
	cc = (rb*rm/ff/ym)**2 - (r0*cosbz)**2
	rc = dsqrt(cc)
	tt = bb*bb - 4.*aa*cc
c	ray penetrates if tt < 0
	range = -999.
	if(tt.lt.0.) go to 40
	cosgam = r0*cosbz/rb
	gamma = dacos(cosgam)
	sing = dsin(gamma)
	ttt = sing + rc/rb + bb/2./rc
	ttt = ttt * ttt * 4. * cc
	tt = dlog (tt/ttt) * r0 * cosbz / 2. / rc
	tt = gamma - beta - tt
	range = 2. * r0 * tt
	beta = beta * 180. / pi
	type*,' freq, elev, range ', fr, beta, range
	write(2,100) f,beta,range
100	format(3f20.2)
40	continue
50	continue
c
	return
	end

