program ION_QP
c	prog to set up a QP ionospheric model for raytracing with the 
c	Jones program.
c
c	To ensure an interpolation scale which is monotonic, a new set of
c	axes is introduced. This is called for convenience the GEOTRANS coord
c	system. It has its north pole at the transmitter. Lines of constant
c	GEOTRANS longitudes and co-latitudes intersect at the required 
c	gridpoints. The longitudes are at +/- nX degrees either side of the
c	longitude of the Rx in GEOTRANS coords, and the co-latitudes are at
c	every Y degrees. The resolutions X and Y must be small enough to 
c	track small-scale structures such as the mid-latitude trough. Note
c	that interpolation between gridpoints is linear.
c
c	written by leo mcnamara 31-jul-90
c-----	last modified 10-aug-90
c
c	The model is defined by a set of profile parameters given at a
c	number of gridpoints covering the part of the ionosphere through
c	which the rays will propagate. It should cover enough area to allow
c	for off-great-circle propagation, and for overshooting the Rx.
c
c	The values of X and its derivatives required by the Jones prog are
c	derived by :-
c	X	Interpolation between gridpoints, using the values of X
c		obtained for a given radius r at the 4 closest surrounding
c		gridpoints. Linear interpolation.
c	pX/pr	As for X. The value of pX/pr is analytic, and merely requires
c		the substitution of the appropriate values of the parameters.
c	pX/pt	Partial derivative wrt theta (co-latitude). Calculate X for 
c		the given r at given theta, and at theta + DELTA theta, 
c		keeping phi (longitude) constant at the given value. Then
c		pX/pt = change of X / DELTA theta.
c	pX/pp	Partial derivative wrt phi (longitude). As for pX/pt.
c
c*****	LINK with files ips,layers,MSIS83,ilios,spoints,erectnh,
c	                profile,gg_gm
c	ips	    loadgp,f2m3,tcon,gyro 
c	msis83	    MSIS83 S/Rs GTS4 etc
c	layers	    Elayer,solar,foeccir
c		    F1layer,moddip,gdmd_interp,fof1OT
c		    F2layer	
c		    Becker [not used - info only]
c	ilios
c	spoints	    spoints,distbear,latlon
c	erectnh	    erectnh,ym,qphfn,qlhfn,qlfnh,qpfnh,parabola,efqp
c	profile	    [only for looking at the N(h) profile]
c	gg_gm
c
c---	Required files
c	1   xgp.***		gridpoint maps of fof2 & m3000 for month ***
c	2   qpnh.dat		output file from s/r PROFILE
c	3   gyro_gp_data.dat
c	4   data.out		generic output file used for debug & profiles
c	7   ionmod_qp.dat	the output file containing the ionospheric model
c
c---	required input data:-
c	month		month
c	iday		UT day
c	ihour		UT hour
c	minute		UT minute
c	tindex		ionospheric index, T
c	r12		12-month smoothed sunspot number
c	tlat		geographic latitude of TX
c	tlon		geographic E longitude of TX
c	rlat		geographic latitude of RX
c	rlon		geographic E longitude of RX
c	extra_range	model is set up for points downrange of the RX by
c			extra_range % of the actual range
c	ng_lat		number of latitude rows to insert between TX and
c			and last latitude of model
c	ng_lon		number of longitudes to set up either side of zero
c			i.e. either side of the TX-RX great circle
c	step_lon	step in longitude between grid longitudes
c	ipro		> 0 if N(h) profile at TX and RX are required 
c			= 0 otherwise
c
c---	output data  .. in the file IONMOD_QP.DAT
c	5 lines of input parameters for identification purposes
c	gmlat_tx,gmlon_tx   geom lat & lon of TX
c	nlat	    number of latitudes in grid
c	grid_colat  colatitudes of grid (GEOTRANS coords)
c	nlon	    number of longitudes of grid
c	grid_lon    GEOTRANS longitudes of grid
c	p(i,j,k)    i=1,nlat
c	p(i,j,k)    j=1,nlon
c	p(i,j,k)    k = 1,19  where
c	k=1	colatitude of point (degrees)
c	k=2	geographic latitude (degrees)
c	k=3	geomagnetic colatitude 
c	k=4	geomagnetic longitude
c	k=5	foE E layer
c	k=6	yme
c	k=7	hme
c	k=8	foF1 F1 layer
c	k=9	ymf1
c	k=10	hmf1
c	k=11	foF2 F2 layer
c	k=12	ymf2
c	k=13	hmf2
c	k=14	hsf
c	k=15	hsf2
c	k=16	fnmax  transition parabola
c	k=17	ymax
c	k=18	hmax
c	k=19	invqp
c
	dimension p(50,5,19),ap(7),nday(12)
	dimension ht(500),fn(500)
	dimension grid_colat(100),grid_lon(50),param(19)
	character*10 munth(12)
	data munth / 'xgp.jan','xgp.feb','xgp.mar','xgp.apr',
     a  'xgp.may','xgp.jun','xgp.jul','xgp.aug','xgp.sep',
     a  'xgp.oct','xgp.nov','xgp.dec' /
	data nday / 31,28,31,30,31,30,31,31,30,31,30,31 /
c
c	grid-point data file on logical unit 1 - see later for OPEN
cc	open (unit=1,file=munth(month),status='old')
	lunmap = 1
c	NOTE that PROFILE writes to unit 2
	open (unit=2,file='qpnh.dat',status='old')
c	the gyrofreq, dip & moddip maps are on unit 3
	open (unit=3,file='gyro_gp_data.dat',form='formatted',
     a        status='old')
	lundip = 3
c	the N(h) profiles are written to DATA.OUT
	open (unit=4,file='data.out',status='old')
c	output is to file IONMOD_QP.DAT
	open (unit=7,file='ionmod_qp.dat',status='old')
c
c	hardwire the year
	lyear = 89
c	radius of Earth - consistent for all QP s/rs
	rz = 6370.
c
c	geographic coords of geomagnetic north pole - consistent with IRI
	gmpole_lat = 78.6
	gmpole_lon = 290.2
c
c	number of profile parameters at each grid point
	nparam = 19
c
10	continue
c
c	initialize the arrays
	ap(1) = 0
	do 1 i = 1,50
	do 1 j = 1,5
	do 1 k = 1,19
1	p(i,j,k) = 0.
c
c	Read in universal time 
	print*,' Enter UT - month, day, hour, minute  - mmddhhmm'
	read(5,500) month, iday, ihour, minute
500	format(5i2)
	if(month.eq.0) stop ' Normal exit'
	write(7,701) month,iday,ihour,minute
701	format(5x,' QP IONOSPHERIC MODEL FOR ',4i2)
c
	open (unit=1,file=munth(month),form='unformatted',status='old')
c
c	Read in sunspot number (R12) and T index
	print*,' Enter sunspot number (R12) and T index'
	read(5,501) r12, tindex
501	format(10f5.0)
	write(7,702) r12,tindex
702	format(5x,' INDICES ', 2f10.1)
c
c	Read in geographic coords of transmitter
	print*,' Enter geographic lat & lon of TX - 2f5.0'
	read(5,501) tlat,tlon
c
	if(tlat+tlon.eq.0.)  stop ' LAT=LON=0 '
	write(7,703) tlat,tlon
703	format(5x,' Geographic coords of TX', 2f10.2)
c
c	Read in geographic coords of receiver
	print*,' Enter geographic lat & lon of RX - 2f5.0'
	read(5,501) rlat,rlon
	write(7,704) rlat,rlon
704	format(5x,' Geographic coords of RX', 2f10.2)
c
c	Read in the % increase in TX-RX range for which we are going to
c	set up the model (to allow for overshoot of the RX)
	print*,' Enter % of range required beyond the RX - f5.0'
	read(5,501) extra_range	
	write(7,705) extra_range
705	format(5x,' Range increased by a % of ',f6.1)
c
c	Read in the number of gridpoints to put along the circuit between
c	the TX and the endpoint
	print*,' Enter number of gridpoints to insert between the'
	print*,' TX and the end of the path - equally spaced in '
	print*,' GEOTRANS latitude'
	read(5,500) ng_lat
c
c	Read in the number of GEOTRANS longitudes, and the separation between
c	them, corresponding to the gridpoints
	print*,' Enter number of longitudes either side of the great'
	print*,' circle to use to define the grid'
	read(5,500) ng_lon
	print*,' Enter separation between longitudes - f5.0'
	read(5,501) step_lon
c
c	ccir formula relating r12 and flux12 (supp 252) - used by MSIS
	flux12 = 63.7 + 0.728 * r12 + 0.00089 * r12 * r12
c	use default flux values if not supplied (set to 0.)
	if(f107.eq.0.) f107 = flux12
	if(f107a.eq.0.) f107a = flux12
c
c	move the TX a little down range (1 km) to ensure that it lies
c	within the grid
	call distbear(tlat,tlon,rlat,rlon,range,bearing)
	call latlon(tlat,tlon,1.,bearing,tlat_mod,tlon_mod)
c
c	geomagnetic coords of TX
	call gg_gm(+1,gmpole_lat,gmpole_lon,tlat,tlon,gmlat_tx,gmlon_tx)
c
	uthour = float(ihour) + float(minute)/60.
c	day number of year ( for MSIS )
	idayno = iday
	if(month.eq.1) go to 12
	mon1 = month - 1
	do 11 i = 1,mon1
11	idayno = idayno + nday(i)
	if(mod(lyear,4).eq.0.and.month.ge.3) idayno = idayno + 1
12	continue
c
	ndm = nday(month)
	if(mod(lyear,4).eq.0.and.month.eq.2) ndm = 29
	lmon = month
c
c	coords of RX in GEOTRANS coords
	call gg_gm (+1,tlat,tlon,rlat,rlon,rmlat,rmlon)
	print*,' GEOTRANS coords of RX are ', rmlat,rmlon
c
c	range from TX to RX
	call distbear (0.,0.,rmlat,rmlon,range,bearing)
c
c	we are going to set up the grid to cover an extra X% of range
	fac = 1. + extra_range / 100.
	call latlon (0.,0.,range*fac,bearing,end_lat,end_lon)
c
c	we are going to set up ng_lat gridpoints between the TX and endpoint
	step_lat = (90. - end_lat) / float(ng_lat + 1)
c
c	number of latitudes in the grid
	nlat = ng_lat + 2
c
c	the GEOTRANS co-latitudes at which the gridpoints are set up
	grid_colat(1) = 0.
	do 200 i = 2,nlat
	grid_colat(i) = grid_colat(i-1) + step_lat
	print*,' colat ',i,grid_colat(i)
200	continue
c
c	the GEOTRANS longitudes at which the gridpoints are to be set up
	grid_lon(1) = rmlon - ng_lon * step_lon
	print*, 'first lon ', grid_lon(1)
	nlon = 2 * ng_lon + 1
	do 210 i = 2,nlon
	grid_lon(i) = grid_lon(i-1) + step_lon
	print*,' lon ', i,grid_lon(i)
210	continue
c
c	calculate the geographic & geomagnetic coords of the gridpoints
	do 230 i = 1,nlat
	grid_lat = 90. - grid_colat(i)
	do 220 j = 1,nlon
	call gg_gm (-1,tlat,tlon,p(i,j,1),p(i,j,2),grid_lat,grid_lon(j))
c	geomagnetic coords
	call gg_gm 
     a  (1,gmpole_lat,gmpole_lon,p(i,j,1),p(i,j,2),p(i,j,3),p(i,j,4))
220	continue
230	continue
c
c	set up the layer and QP parameters for each gridpoint
	do 2001 ii = 1,nlat
	do 2001 jj = 1,nlon
	flat = p(ii,jj,1)
	flon = p(ii,jj,2)
c
c	calculate the local time (from the UT) 
	time = uthour + flon/15.
	lh = time
	lm = (time - lh) * 60 + 0.5
	lday = iday
	if(lh.gt.24) lday = lday + 1
	if(lday.gt.ndm) lmon = mon + 1
	if(lday.gt.ndm) lday = 1
	if(lmon.eq.13) lmon = 1
	if(lh.gt.24) lh = lh - 24
	reflon = flon
	hourlt = float(lh) + float(lm)/60.
c
c	E-layer parameters 
	call Elayer 
     a  (flux12,r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,
     a   chi,chi_noon,foe,yme,hme)
	type*,' E layer ', foe,yme,hme
c
c	F1-layer parameters
	call F1layer (lundip,
     a  r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,
     a  chi_max,fof1,ymf1,hmf1)
c
c	F2-layer parameters
	call F2layer (lunmap,month,uthour,hourlt,idayno,flat,flon,
     a  foe,tindex,f107a,f107,ap,fof2,ymf2,hmf2)
c
	type*,' F1 layer ',fof1,ymf1,hmf1
	type*,' F2 layer ', fof2,ymf2,hmf2,hourlt
c
c	store the layer parameters in the array P
	p(ii,jj,5) = foe
	p(ii,jj,6) = yme
	p(ii,jj,7) = hme
	p(ii,jj,8) = fof1
	p(ii,jj,9) = ymf1
	p(ii,jj,10) = hmf1
	p(ii,jj,11) = fof2
	p(ii,jj,12) = ymf2
	p(ii,jj,13) = hmf2
c
c	calculate the QP profile parameters at each location
	call erectnh 
     a  (rz,foe,yme,hme,fof1,ymf1,hmf1,fof2,ymf2,hmf2,chi,chi_noon,
     a   chi_max,hsf,hsf2,fnmax,hmax,ymax,invqp)
c
c	store the QP parameters
	p(ii,jj,14) = hsf
	p(ii,jj,15) = hsf2
	p(ii,jj,16) = fnmax
	p(ii,jj,17) = ymax
	p(ii,jj,18) = hmax
	p(ii,jj,19) = invqp
c
c	calculate profile at TX, & write to file
	if(ii.eq.1.and.jj.eq.1) go to 1998
	go to 2001
1998	call profile
     a (foe,hme,yme,hsf,fof1,ymf1,hmf1,hsf2,fof2,hmf2,ymf2,rz,
     a  invqp,fnmax,hmax,ymax,ht,fn)
	do 1999 i = 1,500
	if(fn(i).gt.0.) write(4,400) i,ht(i),fn(i)
	if(fn(i).gt.0.) write(6,400) i,ht(i),fn(i)
1999	continue
400	format(i4,f6.1,f6.2)
c
2001	continue
c
c	write out the parameters of the model
	write(7,700) tlat_mod,tlon_mod
	write(7,700) gmlat_tx,gmlon_tx
700	format(5e16.8)
	write(7,706) nlat
706	format(i4)
	write(7,700) (grid_colat(i),i=1,nlat)
	write(7,706) nlon
	write(7,700) (grid_lon(i),i=1,nlon)
	write(7,706) nparam
	do 707 i = 1,nlat
	do 707 j = 1,nlon
	write(7,700) (p(i,j,k),k=1,nparam)
707	continue
c
	go to 10
c
	end