subroutine ray_USU
c
c	s/r to find the values of X and its partial derivatives pxpr,
c	pxpth and pxpph to be used with the Jones program. The parameters
c	of the ionospheric model have previously been set up at a special
c	grid by the program IONMOD_USU. This grid has a north pole at the TX.
c
c	PXPR is found by spline interpolation, PXPTH & PXPPH by 2-D linear
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	it is assumed that raytracing is to be done with the earth's
c	magnetic field included    ---   w(1) = +/- 1.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	written by leo mcnamara 10-jun-91
c-----	  last modified 10-jun-91
c
c	USES s/rs SPLINE & SPLINT_MOD
c
	logical lon_flag
	character*8 modx,id
	real *8 x,freq,hh,pxpr,pxpth,pxpph,pxpt,hmax,r,w0,w
c
	common / xx / modx(2),x,pxpr,pxpth,pxpph,pxpt,hmax
	common r(6) / ww / id(10),w0,w(400)
c
	dimension p(50,5,42),kr(4),kc(4),x_colat(4),y_lon(4),xs(4),
     a  pr(4),hm(4),colat(50),glon(5),alt(37),tem(37),tem2(37),
     a  y2(4,37)
c
	data pi,radeg,degrad / 3.141592654,57.29577951,0.017453292 /
	data ipass / 0 /
	data iclat_last, iclon_last / 0, 0 /
c
c-----------------------------
	entry ELECTX
c-----------------------------
c
	if(w(1).eq.0.) stop ' Must have magnetic field '
c
c	read in the ionospheric model on the first pass
	if(ipass.ne.0) go to 100
	lon_flag = .false.
	open (unit=3,file='ionmod_USU.dat',type='old')
c	read in 5 header lines
	do 1 i = 1,5
1	read(3,333) head
333	format(a4)
c	read in the coords of the TX after it has been moved 1 km down range
	read(3,102) tlat_moved,tlon_moved
c	overwrite any specified TX location
	w(4) = dble(tlat_moved) * degrad
	w(5) = dble(tlon_moved) * degrad
c	read in the geom coords of the TX
	read(3,102) tx_gmlat,tx_gmlon
102	format(5e16.8)
c	the North pole is at the geomag coords of the TX
	pole_lat = tx_gmlat
	pole_lon = tx_gmlon
c	read in the co-latitude values of the latitude circles of the grid
	read(3,101) nlat
101	format(i4)
	read(3,102) (colat(i),i=1,nlat)
	do 1101 i = 1,nlat
1101	print*,' colat ', i,colat(i)
c	read in the longitudes of the grid
	read(3,101) nlon
	read(3,102) (glon(i),i=1,nlon)
c	See if the longitudes straddle zero degrees. If so, set lon_flag 
c	to TRUE and add 360 to the easterly longitudes.
	if(glon(1).gt.340.0.and.glon(nlon).lt.20.) lon_flag = .true.
	do 1102 i = 1,nlon
	if(lon_flag.and.glon(i).lt.20.0) glon(i) = glon(i) + 360.
1102	print*,' glon ', i,glon(i)
c
c	search relies on the COLAT and GLON increasing monotonically - check
c	that this is true. Could have problems if longitude of RX in GEOTRANS
c	coords is around 0 (in spite of fix above! - let's be careful)
	do 104 i = 2,nlat
	if(colat(i).le.colat(i-1)) stop ' Problems with COLAT grid'
104	continue
	do 103 i = 2,nlon
	print*,' glon(i),glon(i-1) ', i,glon(i),glon(i-1)
	if(glon(i).le.glon(i-1)) stop ' Problems with LON grid'
103	continue
c
c	read in the altitudes at which the USU values of FN are given
	read(3,102) (alt(i),i=1,37)
c	
c	read in the parameters of the ionospheric model at each gridpoint
c	there are NPARAM parameters at each grid point
	read(3,101) nparam
	print*,' NPARAM,nlat,nlon ', nparam,nlat,nlon
	do 105 i = 1,nlat
	do 105 j = 1,nlon
105	read(3,102) (p(i,j,l),l = 1,nparam)
c
c	geographic coords of TX
	gglat_tx = p(1,1,1)
	gglon_tx = p(1,1,2)
	print*,'  Geog coords of TX ',gglat_tx,gglon_tx
c
	ipass = 1
c	return on first call, to avoid problems with rest of s/r
	return
c
100	continue
c
	freq = w(6)
	frequency = sngl(freq)
	hh = r(1) - w(2)
	height = sngl(hh)
c	coords of point - r,theta,phi - work in degrees
	rr = sngl(r(1))
	theta = sngl(r(2)) * radeg
	phi = sngl(r(3)) * radeg
	if(phi.lt.0.) phi = phi + 360.
	gmlat = 90. - theta
	gmlon = phi
ccc	print*,' RAY_QP rr, theta,phi ', rr-6370.,theta,phi
c
c	geog coords - we need this, since we cannot go directly from
c	GEOM coords to GEOTRANS coords (too dumb)
	gmpole_lat = 78.6
	gmpole_lon = 290.2
	call gg_gm (-1,gmpole_lat,gmpole_lon,gglat,gglon,90.-theta,phi)
ccc	print*,' geog coords ',rr,gglat,gglon
c
c	return if point is below the ionosphere
	x = 0.
	pxpr = 0.
	pxpth = 0.
	pxpph = 0.
	pxpt = 0.
	if(height.lt.alt(1)) return
c
c	coords of point in GEOTRANS coord frame (North pole at TX)
c	work from geog coords (since we can't work out how to go directly
c	from geom to geotrans coords)
	call gg_gm (1,gglat_tx,gglon_tx,gglat,gglon,xlat,xlon)
c	modify xlon if the longitudes straddle zero
	if(lon_flag.and.xlon.lt.20.0) xlon = xlon + 360.
ccc	write(6,666)  hh,90.-xlat,xlon
666	format(2x,'GEOTRANS location ', 3f10.2)
c
c	find the gridpoints surrounding (xlat,xlon)
c	co-latitude of the point (xlat,xlon)
	xcolat = 90. - xlat
	nlatm1 = nlat - 1
	do 200 i=1,nlatm1
	iclat = i
	if(xcolat.gt.colat(i).and.xcolat.le.colat(i+1)) go to 222
200	continue
	print*,' xcolat,colat(1),colat(nlat) ', 
     a           xcolat,colat(1),colat(nlat)
	stop ' Point outside grid values of co-latitude'
c
222	continue
	nlonm1 = nlon - 1
	do 250 i = 1,nlonm1
	ilon = i
	if(xlon.gt.glon(i).and.xlon.le.glon(i+1)) go to 226
250	continue
	print*,' xlon,glon(1),glon(nlon) ',
     a           xlon,glon(1),glon(nlon)
	stop ' Point outside grid values of longitude'
226	continue
c
ccc	print*,' iclat,ilon ', iclat,ilon
c
c	point is surrounded by the 4 gridpoints defined by colat(iclat) 
c	and colat(iclat+1) and glon(ilon) and glon(ilon+1)
	kr(1) = iclat
	kr(2) = iclat
	kr(3) = iclat + 1
	kr(4) = iclat + 1
	kc(1) = ilon
	kc(2) = ilon + 1
	kc(3) = ilon + 1
	kc(4) = ilon
c
c	Set up arrays of second derivatives if we have entered a new
c	section of the grid
	if(iclat.eq.iclat_last.and.iclon.eq.iclon_last) go to 50
c	We are in a new section of the grid - set up FN array (TEM)
	do 30 j = 1,4
	do 20 i = 1,37
20	tem(i) = p(kr(j),kc(j),i+5)
c	calculate the second derivatives (Natural spline --> 1.01e30)
	call spline (alt,tem,37,1.01e30,1.01e30,tem2)
	do 25 i = 1,37
25	y2(j,i) = tem2(i)
30	continue
	iclat_last = iclat
	iclon_last = iclon
50	continue
c
c	Interpolate to get values of X, pX/pr and hmF2
	do 300 k = 1,4
c	set up model parameters at each grid point k1,k2 (row,col)
	k1 = kr(k)
	k2 = kc(k)
	do 60 i = 1,37
c	plasma frequency
	tem(i) = p(k1,k2,i+5)
c	second derivative
60	tem2(i) = y2(k,i)
	call splint_mod (alt,tem,tem2,37,height,fn,dfndh)
	xs(k) = (fn/frequency)**2
	pr(k) = dfndh * 2 * fn / frequency**2
	hm(k) = p(k1,k2,5)
c
300	continue
c
ccc	print*,' x ', (xs(k),k=1,4)
c	return if all 4 values of X = 0
	if(xs(1)+xs(2)+xs(3)+xs(4).lt.1.e-10) return
c
c	interpolate at the point P to get x, hmax and pxpr (2-D interpolation)
c	the lat and lon of the 4 corners
	x_colat(1) =  colat(iclat)
	x_colat(2) =  colat(iclat)
	x_colat(3) =  colat(iclat + 1)
	x_colat(4) =  colat(iclat + 1)
	y_lon(1) = glon(ilon)
	y_lon(2) = glon(ilon + 1)
	y_lon(3) = glon(ilon + 1)
	y_lon(4) = glon(ilon)
ccc	print*,' XC ',(x_colat(i),i=1,4)
ccc	print*,' YC ',(y_lon(i),i=1,4)
	call interp_2d (xcolat,xlon,x_colat,y_lon,xs,ans_x)
	call interp_2d (xcolat,xlon,x_colat,y_lon,pr,ans_px)
	call interp_2d (xcolat,xlon,x_colat,y_lon,hm,ans_hm)
	x = dble(ans_x)
ccc	print*,' Interpolated X ', x
	pxpr = dble(ans_px)
	hmax = dble(ans_hm)
c
c	calculate the gradients wrt theta and phi, by finding X at a nearby
c	point P1 with (1) theta + DELTA_theta, same phi, ==> pxpth
c	point P2 with (2) theta, phi + DELTA_phi  ==> pxpph
c	The steps DELTA are made small enough that it is REASONABLY safe 
c	to use the same interpolation grid as before 
	delta = 0.2
c	we need the size of the interpolation box in geom coords
	size_gmlat = abs(p(iclat,ilon,3) - p(iclat+1,ilon,3))	
c	all the longitudes are identical at the TX (iclat = 1)
	jclat = iclat
	if(iclat.eq.1) jclat = 2
	size_gmlon = abs(p(jclat,ilon,4) - p(jclat,ilon+1,4))	
	delta_theta = delta * size_gmlat
	delta_phi = delta * size_gmlon
c	find GEOTRANS coords of the two points, from the geom coords,
c	via the geog coords
	gmlat1 = gmlat + delta_theta
	gmlon1 = gmlon 
	gmlat2 = gmlat
	gmlon2 = gmlon + delta_phi
	call gg_gm 
     a  (-1,gmpole_lat,gmpole_lon,gglat1,gglon1,gmlat1,gmlon1)
	call gg_gm 
     a  (-1,gmpole_lat,gmpole_lon,gglat2,gglon2,gmlat2,gmlon2)
	call gg_gm (1,gglat_tx,gglon_tx,gglat1,gglon1,xlat1,xlon1)
	call gg_gm (1,gglat_tx,gglon_tx,gglat2,gglon2,xlat2,xlon2)
c
	xcolat1 = 90. - xlat1
	xcolat2 = 90. - xlat2
ccc	print*,' GEOTRANS 1 & 2 ', xcolat1,xlon1,xcolat2,xlon2
ccc	print*,' x ', (xs(k),k=1,4)
	call interp_2d (xcolat1,xlon1,x_colat,y_lon,xs,ans_x1)
	call interp_2d (xcolat2,xlon2,x_colat,y_lon,xs,ans_x2)
c
	pxpth = dble(ans_x1 - ans_x) / (delta_theta * degrad)
c	the gradient wrt co-latitude = - gradient wrt latitude
	pxpth = - pxpth
	pxpph = dble(ans_x2 - ans_x) / (delta_phi * degrad)
c
ccc	print*, ans_x,ans_x1,ans_x2,pxpth,pxpph
c
c	there is no variation with time
	pxpt =0.
c
	return
	end
c------------------------------------------------------------------
	subroutine interp_2d (xp,yp,xc,yc,z,zp)
	dimension xc(1),yc(1),z(1)
	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)
	za = z(1)
	zb = z(2)
	zc = z(3)
	zd = z(4)
	fac1 = (yp - ya) / (yb - ya)
	fac2 = (xp - xa) / (xd - xa)
	zp = za + fac1*(zb-za) + fac2*(zd-za+fac1*(za-zb+zc-zd))
	return
	end