subroutine dud_3d
c	s/r to find the values of X and the partial derivatives PXPR,
c	PXPTH and PXPHI to be used by the Jones program, given a gridded
c	array of the parameters required to erect a Dudeney profile.
c	The array is set up by the program IONMOD_DUD  --  see the listing
c	to see what the grid is, and how the array of parameters is set up.
c
c	first call reads in the array P
c
c	written by leo mcnamara 3-jun-87
c-----	last modified 5-nov-87
c
c	s/r DUDNEY is in the file DUD_JONES
c
	real m3000
	logical lcond
	character*8 modx,id
	real*8 x,freq,hh,pxpr,pxpth,pxpph,pxpt,hmax,r,w0,w,ffiray
c
	common / iondat / fof2,fof1,foe,hff2,m3000,hme,yme,hmf2,ymf2,
     a                    h0,xd,dv
	common / parms / beta,ht,ft,at,hf,ad,tau,hd,lcond
	common / xx / modx(2),x,pxpr,pxpth,pxpph,pxpt,hmax
	common r(6) / ww / id(10),w0,w(400)
	dimension p(50,5,20),colat(50),elon(50),xc(4),yc(4),kk(4,2),
     a            xs(4),pr(4),hm(4)
	data ipass / 0 /
c
c---------------------------------
	ENTRY ELECTX
c---------------------------------
c
	pi = 3.141592654
	radeg = 180. / pi
	degrad = pi / 180.
c
c	read in the array P on the first call
	if(ipass.eq.1) go to 50
	open (unit=3,file='ionmod_dud.dat',status='old')
c	skip first 5 lines
	do 20 i = 1,5
20	read(3,300) trash
300	format(a1)
	read(3,301) nop
301	format(i3)
	type*,' # of profiles ',nop
	read(3,301) iray
	firay = float(iray)
	ffiray = dble (firay)
	if(dabs(ffiray-dabs(w(1))).le.1.d-5) go to 1112
	type*,' INCONSISTENT MODEL GRID AND W(1) ',iray,w(1)
	stop
1112	read(3,301) ispec
	do 333 i = 1,nop
	do 333 j = 1,5
333	read(3,334) (p(i,j,k),k=1,20)
334	format(5e16.6)
c
c	arrange to use geographic or geomagnetic coords
	ir = 1
	if(iray.eq.1) ir = 3
	jr = 2
	if(iray.eq.1) jr = 4
c
c	the colatitudes and longitudes of the points along the circuit
c	(including the TX and RX)
c	For E-W propagation, row 3 is a line of latitude, NOT the
c	great circle
	jray = 2 * iray
	do 30 i = 1,nop
	colat(i) = p(i,3,jray+1)
	elon(i) = p(i,3,jray+2)
	if(elon(i).lt.0.) elon(i) = elon(i) + 360.
	type*,'i,colat,elon',i, colat(i),elon(i)
30	continue
c
	ipass = 1
	return
50	continue
c
	freq = w(6)
	hh = r(1) - w(2)
c	coords of point - r,theta,phi - work in degrees
	rr = sngl (r(1))
	th = sngl (r(2)) * radeg
	ph = sngl (r(3)) * radeg
	if(ph.lt.0.) ph = ph + 360.
c
c	print out locations when ray hits the ground
	tz = 90. - th 
	pz = ph 
	if(dabs(r(1)-w(2)).le.0.01) type*,' r,theta,phi ',rr,tz,pz
	if(dabs(r(1)-w(2)).le.0.01.and.abs(w(1)).eq.1.d0)
     a            call ggm (1,pz,tz,ph,90.-th)
	if(dabs(r(1)-w(2)).le.0.01.and.abs(w(1)).eq.1.d0)
     a            type*,' r,theta,phi ',rr,tz,pz
c
	nop1 = nop - 1
	if(ispec.lt.0) go to 59
c
	if(ispec.eq.0.or.ispec.eq.180) go to 230
c	E-W or W-E propagation 
	do 200 i = 1,nop1
	jj = i
	if(ph.ge.elon(i).and.ph.le.elon(i+1)) go to 210
	if(ph.ge.elon(i+1).and.ph.le.elon(i)) go to 210
200	continue
	type*,' Phi out of bounds s/n 200', ph,elon(1),elon(nop)
	go to 6666
210	continue
c	find the latitudes above and below the point (latitudes at all
c	""jj" are in fact the same
	do 211 i = 1,4
	ii = i
	if(th.ge.p(jj,i,ir).and.th.le.p(jj,i+1,ir)) go to 215
	if(th.ge.p(jj,i+1,ir).and.th.le.p(jj,i,ir)) go to 215
211	continue
	type*,' theta out of range s/n 211', th,p(jj,1,ir),p(jj,5,ir)
	go to 6666
215	continue
c
c	the 4 corners of the box surrounding the point P (for interpolation)
c	The latitudes all correspond to those at the end of the circuit, but 
c	we will leave the coding general
	xc(1) = p(jj,ii,ir)
	yc(1) = elon(jj)
	xc(2) = p(jj+1,ii,ir)
	yc(2) = elon(jj+1)
	xc(3) = p(jj+1,ii+1,ir)
	yc(3) = elon(jj+1)
	xc(4) = p(jj,ii+1,ir)
	yc(4) = elon(jj)
c	the points at which the Dudeney parameters are required
	kk(1,1) = jj
	kk(1,2) = ii
	kk(2,1) = jj + 1
	kk(2,2) = ii
	kk(3,1) = jj + 1
	kk(3,2) = ii + 1
	kk(4,1) = jj
	kk(4,2) = ii + 1
ccc	type*,' 4 corners ',(xc(i),yc(i),i=1,4)
ccc	do 9000 i = 1,4
ccc9000	type*,' the array kk ',(kk(i,j),j=1,2)
	go to 87
c
c	N-S or S-N propagation
230	continue
	do 400 i = 1,4
	jj = i
	if(ph.ge.p(nop,i,jr).and.ph.le.p(nop,i+1,jr)) go to 410
	if(ph.ge.p(nop,i+1,jr).and.ph.le.p(nop,i,jr)) go to 410
400	continue
	type*,' Phi out of bounds s/n 400', ph,p(nop,1,jr),p(nop,5,jr)
	go to 6666
410	continue
c	find the latitudes above and below the point
	do 411 i = 1,nop
	ii = i
	if(th.ge.colat(i).and.th.le.colat(i+1)) go to 415
	if(th.ge.colat(i+1).and.th.le.colat(i)) go to 415
411	continue
	type*,' theta out of range s/n 411', th,colat(1),colat(nop)
	go to 6666
415	continue
c
c	the 4 corners of the box surrounding the point P (for interpolation)
c	The longitudes all correspond to those at the end of the circuit, but 
c	we will leave the coding general
	xc(1) = colat(ii+1)
	yc(1) = p(ii+1,jj,jr)
	xc(2) = colat(ii+1)
	yc(2) = p(ii+1,jj+1,jr)
	xc(3) = colat(ii)
	yc(3) = p(ii,jj+1,jr)
	xc(4) = colat(ii)
	yc(4) = p(ii,jj,jr)
c	the points at which the Dudeney parameters are required
	kk(1,1) = ii + 1
	kk(1,2) = jj
	kk(2,1) = ii + 1
	kk(2,2) = jj + 1
	kk(3,1) = ii
	kk(3,2) = jj + 1
	kk(4,1) = ii 
	kk(4,2) = jj
ccc	type*,' 4 corners ',(xc(i),yc(i),i=1,4)
ccc	do 9000 i = 1,4
ccc9000	type*,' the array kk ',(kk(i,j),j=1,2)
	go to 87
c
c	GENERAL CASE  - find coords of box surrounding the given point P
59	do 60 i = 1,nop1
	ii = i
	if(th.ge.colat(i).and.th.le.colat(i+1)) go to 65
	if(th.ge.colat(i+1).and.th.le.colat(i)) go to 65
60	continue
	type*,' theta out of bounds -- theta,colat(1),colat(nop)',
     a          th,colat(1),colat(nop)
	go to 6666
65	do 70 i = 1,nop1
	jj = i
	if(ph.ge.elon(i).and.ph.le.elon(i+1)) go to 75
	if(ph.ge.elon(i+1).and.ph.le.elon(i)) go to 75
70	continue
	type*,' phi out of bounds -- phi,elon(1),elon(nop)',
     a          ph,elon(1),elon(nop)
	go to 6666
75	continue
ccc	type*,' II & JJ  ',ii,jj
c
c	the 4 corners of the box (clockwise, from top left)
	xc(1) = colat(ii+1)
	yc(1) = elon(jj)
	xc(2) = colat(ii+1)
	yc(2) = elon(jj+1)
	xc(3) = colat(ii)
	yc(3) = elon(jj+1)
	xc(4) = colat(ii)
	yc(4) = elon(jj)
ccc	type*,' 4 corners ',(xc(i),yc(i),i=1,4)
c
c	the points along the rows of the grid at which the Dudeney parameters
c	are required.  kk(l,m) is used for P(kk(l,1),kk(l,2),n)  --  3 cases
c	case 1    ii = jj
	if(ii.ne.jj) go to 80
	kk(1,1) = ii
	kk(1,2) = 2
	kk(2,1) = ii+1
	kk(2,2) = 3
	kk(3,1) = ii
	kk(3,2) = 4
	kk(4,1) = ii
	kk(4,2) = 3
	go to 90
c
80	if(ii.ne.jj-1) go to 85
c	case 2   ii = jj - 1
	kk(1,1) = ii+1
	kk(1,2) = 3
	kk(2,1) = ii + 1
	kk(2,2) = 4
	kk(3,1) = ii
	kk(3,2) = 5
	kk(4,1) = ii
	kk(4,2) = 4
	go to 90
c
c	case 3    ii = jj + 1
85	if(ii.ne.jj+1) type*,'ii,jj+1,  theta,phi ',ii,jj+1,th,ph
	if(ii.ne.jj+1) go to 6666
	kk(1,1) = ii-1
	kk(1,2) = 1
	kk(2,1) = ii
	kk(2,2) = 2
	kk(3,1) = ii
	kk(3,2) = 3
	kk(4,1) = ii-1
	kk(4,2) = 2
c
90	continue
ccc	do 900 i = 1,4
ccc900	type*,' the array kk ',(kk(i,j),j=1,2)
c
c	set up the Dudeney parameters for each corner, and find the values
c	of X and PXPR at the height HH
87	do 100 kc = 1,4
	k1 = kk(kc,1)
	k2 = kk(kc,2)
	foe =	p(k1,k2,5)
	yme =	p(k1,k2,6)
	hme=	p(k1,k2,7)
	fof1 =	p(k1,k2,8)
	fof2 =	p(k1,k2,9)
	ymf2=	p(k1,k2,10)
	hmf2=	p(k1,k2,11)
	h0=	p(k1,k2,12)
	ad=	p(k1,k2,13)
	tau=	p(k1,k2,14)
	hd =	p(k1,k2,15)
	beta=	p(k1,k2,16)
	ht =	p(k1,k2,17)
	ft =	p(k1,k2,18)
	at =	p(k1,k2,19)
c	we don't want to set up the parameters again
	jpass = 1
	call dudney (jpass,freq,hh)
c	results are returned in /XX/
	xs(kc) = sngl (x)
	pr(kc) = sngl (pxpr)
	hm(kc) = hmf2
100	continue
ccc	type*,' X',(xs(i),i=1,4)
ccc	type*,' PXPR',(pr(i),i=1,4)
ccc	type*,' hmF2',(hm(i),i=1,4)
c
c	interpolate at the point P to get X, HMAX & PXPR (2-D interpolation)
	xa = xc(1)
	xd = xc(4)
	ya = yc(1)
	yb = yc(2)
	if(xa.eq.xd) type*,'theta,phi ',th,ph
	if(ya.eq.yb) type*,'theta,phi ',th,ph
	if(xa.eq.xd) type*,' xa & xd are equal --- ',xa
	if(ya.eq.yb) type*, 'ya & yb are equal --- ',ya
	if(xa.eq.xd) type*,' xc', (xc(i),i=1,4)
	if(xa.eq.xd) type*,' yc', (yc(i),i=1,4)
	if(ya.eq.yb) type*,' xc ',(xc(i),i=1,4)
	if(ya.eq.yb) type*,' yc ',(yc(i),i=1,4)
	xp = th
	yp = ph
	za = xs(1)
	zb = xs(2)
	zc = xs(3)
	zd = xs(4)
	fac1 = (yp - ya) / (yb - ya)
	fac2 = (xp - xa) / (xd - xa)
	zp = za + fac1*(zb-za) + fac2*(zd-za+fac1*(za-zb+zc-zd))
	x = dble (zp)
c	calculate the product mu * r for Bouger's rule
	aamu = sqrt(1. - zp)
	h = sngl (hh)
	aar = h + 6370.
	bouger = aamu * aar
cccc	type*,' Bouger -- h, mu, mu * r ', h, aamu, bouger
	za = hm(1)
	zb = hm(2)
	zc = hm(3)
	zd = hm(4)
	zp = za + fac1*(zb-za) + fac2*(zd-za+fac1*(za-zb+zc-zd))
	hmax = dble (zp)
	za = pr(1)
	zb = pr(2)
	zc = pr(3)
	zd = pr(4)
	zp = za + fac1*(zb-za) + fac2*(zd-za+fac1*(za-zb+zc-zd))
	pxpr = dble (zp)
c
c	the partial derivatives wrt theta and phi -- the second derivatives
c	are assumed constant over the box
	pp1 = (xs(1) - xs(4)) / (xc(1) - xc(4))
	pp2 = (xs(2) - xs(3)) / (xc(2) - xc(3))
	pp12 = pp1 + (pp2-pp1) * (yc(1) - ph) / (yc(1) - yc(2))
	pxpth = dble (pp12 * radeg)
ccc	type*,' pp1,pp2,pp12 ',pp1,pp2,pp12
	pt1 = (xs(2) - xs(1)) / (yc(2) - yc(1))
	pt2 = (xs(3) - xs(4)) / (yc(3) - yc(4))
	pt12 = pt1 + (pt2-pt1) * (xc(1) - th) / (xc(1) - xc(4))
ccc	type*,' pt1  pt2  pt12 ',pt1,pt2,pt12
	pxpph = dble (pt12 * radeg)
c
	return
6666	call ggm (1,zlon,zlat,ph,90.-th)
	write(2,6667) zlat,zlon,90.-th,ph
	write(6,6667) zlat,zlon,90.-th,ph
6667	format(//3x,'The ray has gone out of the grid provided',
     a	 /'    at GEOGRAPHIC  coords ',2f8.2,
     a   /'    &  GEOMAGNETIC coords ',2f8.2)
	stop '  TRY it again, Sam .....'
	end