program ionmod_IRI
c
c	prog to determine IRI profile parameters at a grid of
c	locations, for use in 3-D raytracing using the Jones program
c
c	The grid is rectangular in GEOGRAPHIC coords for NO magnetic
c	                           field raytracing, (iray=0)
c	and                     in GEOMAGNETIC coords for raytracing
c	                           with magnetic field (iray=1)
c
c	written by Leo McNamara 10-jul-87
c---	last modified 11-aug-87
c
c	LINK with files ips,fof2m3_ips,ilios,spoints,irinhpro,
c	                irinhsubs,xesubs,ggm,solar
c	ips	    loadgp,f2m3,tcon,gyro 
c	fof2m3_ips  fof2m3_ips
c	ilios
c	spoints	    spoints,distbear,latlon
c	irinhpro    irinhpro, iondem (currently a DUMMY file)
c	irinhsubs   all the IRI s/rs required for n(h)
c	ggm
c	solar
c
c	The model of the ionosphere is set up between the first end point
c	and the second end point. These points should be chosen so that
c	(1) the TX lies within the body of the grid
c	(2) all rays have all hops without leaving the grid at the
c	    second end point
c	if (1) is not satisfied, the prog will bomb out - be careful
c	   near the magnetic equator
c	if (2) is not satisfied, an error message is printed by RAYTRACE
c	   but the prog will not bomb
c
c	There are N-2 points along the circuit (N including the TX & RX)
c	These are the points (i,3) i=1,n
c	For the general case (see later for special cases),
c	there is a row of points (i,1) for which lat(i,1) = lat(i+2,3)
c	                                         lon(i,1) = lon(i,3)
c 	         a row of points (i,2) for which lat(i,2) = lat(i+1,3)
c	                                         lon(i,2) = lon(i,3)
c	         a row of points (i,4) for which lat(i,4) = lat(i,3)
c	                                         lon(i,4) = lon(i+1,3)
c	and a fifth row for which lat(i,5) = lat(i,3)
c	                          lon(i,5) = lon(i+2,3)
c	This grid facilitates the calculation of the partial derivatives
c	of X wrt theta (constant r & phi) and phi (constant r & theta) 
c	in ELECTX
c	p(n,i,j) j=1,npr are the profile parameters at each point (n,i)
c
c	There are 4 SPECIAL CASES - when propagation is, or is nearly,
c	N,S,E or W. In these cases, the above grid collapses and rectangular
c	grids lying N - S or E - W must be set up. A special case is invoked
c	if the difference from any of the 4 cardinal directions is less than
c	the parameter "BDIFF". The grid lines will be STPLON or STPLAT apart
c	and follow lines of longitude & latitude (NOT great circles)
c
c---	required input data:-
c	mon	month
c	iday	UT day
c	ih	UT hour
c	min	UT minute
c	tindex	ionospheric index, T
c	r12	12-month smoothed sunspot number
c	glat1	geographic latitude of start point of circuit
c	glon1	geographic E longitude of first point
c	glat2	geographic latitude of end point of circuit
c	glon2	geographic E longitude of end point
c	nip	number of profiles to set up between start & end points
c	iray	=0 for grid in geographic coords
c		=1 for grid in geomagnetic coords
c	ipro	> 0 if N(h) profiles are required (only 3*5 profiles)
c		= 0 otherwise
c
c---	output data  .. in the file IONMOD_IRI.DAT
c	5 lines of headings etc
c	nop	total number of points along circuit
c	iray
c	ispec   special case azimuths - 0,90,180,270; -98 for general case
c	p(i,j,k)    i=1,nip+2 points along the circuit
c	p(i,j,k)    j=1,5 "rows" of data points (parallel)
c	p(i,j,k)    k = 1,30  where
c	k=1	colatitude of point (degrees)
c	k=2	geographic latitude (degrees)
c	k=3	geomagnetic colatitude (degrees)
c	k=4	geomagnetic longitude
c	k=5	hmf2   ---  the IRI parameters
c	k=6	nmf2
c	k=7	hmf1
c	k=8	b0
c	k=9	b1
c	k=10	c1
c	k=11	hz
c	k=12	t
c	k=13	hst
c	k=14	hdx
c	k=15	hme
c	k=16	nme
c	k=17	hmd
c	k=18	nmd
c	k=19	hef
c	k=20	d1
c	k=21	xkk
c	k=22	fp30
c	k=23	fpu
c	k=24	fp1
c	k=25	fp2
c	k=26	night
c	k=27	e(1)
c	k=28	e(2)
c	k=29	e(3)
c	k=30	e(4)
c
	dimension p(50,5,30),flatt(50),flonn(50),nday(12),e(4)
	dimension fn(3,5,275)
	character*10 munth(12)
	character*8 modx
      LOGICAL LCOND,night
      REAL M3000,nmd,nme,nmf2
c	IRI labelled common
      COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR
     &  /BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,XKK,FP30,FP3U,FP1,FP2
     &  /BLOCK5/ZX,TNX,ATN,CTN/BLOCK6/NIGHT,E/BLOTE/AHH(7),ATE1,
     &  STTE(6),DTE(5)/CONST/UMR
     &  /BLOCK8/HS,TNHS,XSM,MM,MXSM/BLO10/BETA,ETA,DELTA,ZETA
	data munth / 'xgp.jan','xgp.feb','xgp.mar','xgp.apr',
     a  'xgp.may','xgp.jun','xgp.jul','xgp.aug','xgpsin.sep',
     a  'xgp.oct','xgp.nov','xgp.dec' /
	data nday / 31,28,31,30,31,30,31,31,30,31,30,31 /
c
c	critical difference in bearing from cardinal direction
	bdiff = 8.0
c	lat and lon spacing of grid for special cases
	stplon = 2.0
	stplat = 2.0
c
	dv = 0.1
c
c	size of array dimension
	ndim = 50
c	set the year in concrete
	lyear = 87
c
c	output is to file IONMOD_IRI.DAT
	open (unit=1,file='ionmod_iri.dat',status='old')
c	grid-point data file on logical unit 2 - see later for OPEN
cc	open (unit=2,file=munth(mon),status='old')
	lunmap = 2
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
c	initialize the arrays
	do 1 i = 1,ndim
	do 1 j = 1,5
	do 1 k = 1,30
1	p(i,j,k) = 0.
	do 2 i = 1,3
	do 2 j = 1,5
	do 2 k = 1,275
2	fn(i,j,k) = 0.
c
	write(6,600)
600	format(//5x,'Prog to set up the profile parameters for',/
     a  '     a 3-D ionospheric model, using the IRI86 profile'//)
c
c	read in the data
10	write(6,100)
100	format(//2x,'Enter UNIVERSAL time- month,day,hour,min',
     a  ' ... MMDDHHmm',
     a  //'  hit RETURN to stop the program')
	read(5,200) mon,iday,ih,min
200	format(4i2)
	if(mon.eq.0) stop ' ..... all done'
c
	open (unit=lunmap,file=munth(mon),status='old')
c
	write(6,101)
101	format(//2x,'Enter ionospheric index (T) & sunspot number',
     a	' (R12) .. format 2f5.0')
	read(5,201) tindex,r12
201	format(4f5.0)
	write(6,102)
102	format(//2x,'Enter geographic latitude & EAST longitude of',
     a ' ends of circuit   ....  4f5.0')
	read(5,201) glat1,glon1,glat2,glon2
1022	write(6,103)
103	format(//2x,'Enter number of points BETWEEN end points -- i2')
	read(5,104) nip
104	format(i2)
	if(nip.le.ndim-4) go to 110
	write(6,111) ndim-4
111	format(//'  ***** you can have only ',i3,' points'/)
	go to 1022
110	continue
	write(6,666)
666	format(//' Enter 0 for geographic grid; 1 for ',
     a ' geomagnetic grid')
	read(5,104) iray
	write(6,112)
112	format(//' Enter a positive integer if you want to calculate',
     a ' N(h) profiles for a lookie-loo ')
	read(5,104) ipro
c
c	If iray =1, the grid is to be geomagnetic 
	flat1 = glat1
	flon1 = glon1
	flat2 = glat2
	flon2 = glon2
	if(iray.eq.0) go to 5
	call ggm (0,glon1,glat1,flon1,flat1)
	call ggm (0,glon2,glat2,flon2,flat2)
	type*,' Geom coords ',flat1,flon1,flat2,flon2
c	FLAT & FLON are now geomagnetic
5	continue
c
c	check for problems because of bearing
	call distbear (flat1,flon1,flat2,flon2,distsv,bearsv)
	type*,' Distsv & bearsv ',distsv,bearsv
c	if bearing is less than BDIFF from a cardinal direction, use the
c	corresponding special case
	special = -99.
	if(abs(bearsv-0.00).le.bdiff) special = 0.
	if(abs(bearsv-90.0).le.bdiff) special = 90.
	if(abs(bearsv-180.).le.bdiff) special = 180.
	if(abs(bearsv-270.).le.bdiff) special = 270.
c	recalculate (OVERWRITE) the coords of the second end of the circuit
	if(special.gt.-98.0)
     a    call latlon (flat1,flon1,distsv,special,flat2,flon2)
	if(special.gt.-98.0)
     a    type*,' Coords of end OVERWRITTEN to ',flat2,flon2
c
	write(1,298) mon,iday,ih,min,tindex,r12,flat1,flon1,flat2,
     a  flon2,nip
	write(4,298) mon,iday,ih,min,tindex,r12,flat1,flon1,flat2,
     a  flon2,nip
298	format(/5x,'UT mon,day,hour,min = ',4i3,5x,/
     a  5x,'tindex,r12 =',2f5.1,5x,'lat & long =',4f6.2,i5, ' points')
c
c	ccir formula relating r12 and flux12 (supp 252)
	flux12 = 63.7 + 0.728 * r12 + 0.00089 * r12 * r12
c	radius of Earth
	rz = 6370.
	uthour = float(ih) + float(min)/60.
c
	ndm = nday (mon)
	if(mod(lyear,4).eq.0.and.mon.eq.2) ndm = 29
	lmon = mon
c
c	calculate the positions of the NIP points along the circuit
	type*,' flat1,flon1,flat2,flon2 ', flat1,flon1,flat2,flon2
	call spoints (flat1,flon1,flat2,flon2,nip,flatt,flonn)
	type*,' flat1,flon1,flat2,flon2 ', flat1,flon1,flat2,flon2
c
c	total number of points along circuit, including TX and RX
	nop = nip + 2
c
c	store the geogc lat & lon of the points along the circuit
	p(1,3,1) = flat1
	p(nop,3,1) = flat2
	p(1,3,2) = flon1
	p(nop,3,2) = flon2
	do 1000	 i = 1,nip
	p(i+1,3,1) = flatt(i)
1000	p(i+1,3,2) = flonn(i)
c
c	for the last 2 sets, need the coords of 2 points past the rx,
c	with same spacing
	call distbear (flatt(nip),flonn(nip),flat2,flon2,dist,bear)
	call latlon (flat2,flon2,dist,bear,x1,x2)
	call latlon (flat2,flon2,2.*dist,bear,x3,x4)
	p(nop+1,3,1) = x1
	p(nop+1,3,2) = x2
	p(nop+2,3,1) = x3
	p(nop+2,3,2) = x4
c
c	check for special bearing case
	if(special.ge.0.) go to 205
c
c	store the geogc lat & lon of the other 4 rows of points
	do 1100 i = 1,nop
	p(i,1,1) = p(i+2,3,1)
	p(i,1,2) = p(i,3,2)
	p(i,2,1) = p(i+1,3,1)
	p(i,2,2) = p(i,3,2)
	p(i,4,1) = p(i,3,1)
	p(i,4,2) = p(i+1,3,2)
	p(i,5,1) = p(i,3,1)
1100	p(i,5,2) = p(i+2,3,2)
	go to 220
c
205	if(special.eq.90.0.or.special.eq.270.0) go to 206
c	required propagation is almost N-S or S-N
	if(special.eq.0.0.or.special.eq.180.0) go to 2005
	type*,' Special case error at 2005; special = ',special
2005	continue
c	Note longitudes are all the same as those at the TX end
	do 1101 i = 1,nop
	p(i,1,1) = p(i,3,1)
	p(i,1,2) = p(1,3,2) - 2.*stplon
	p(i,2,1) = p(i,3,1)
	p(i,2,2) = p(1,3,2) - stplon
	p(i,4,1) = p(i,3,1)
	p(i,4,2) = p(1,3,2) + stplon
	p(i,5,1) = p(i,3,1)
1101	p(i,5,2) = p(1,3,2) + 2.*stplon
	go to 220
c	required propagation is E-W or W-E
c	NOTE latitudes are referred to those at the TX site (end of circuit)
206	do 1102 i = 1,nop
	p(i,1,1) = p(1,3,1) + 2.*stplat
	p(i,1,2) = p(i,3,2)
	p(i,2,1) = p(1,3,1) + stplat
	p(i,2,2) = p(i,3,2)
c	the grid lines do not follow the great circle, so we cannot use
c	the previously calculated value of p(i,3,1)
	p(i,3,1) = p(1,3,1) + 0.
	p(i,4,1) = p(1,3,1) - stplat
	p(i,4,2) = p(i,3,2)
	p(i,5,1) = p(1,3,1) - 2.*stplat
1102	p(i,5,2) = p(i,3,2)
c
220	continue
c
c	if the locations are geomagnetic, move them into the right spaces
c	in the arrays, and fill in the geographic values
	do 240 i = 1,nop
	do 240 j = 1,5
	if(iray.ne.0) p(i,j,3) = p(i,j,1)
	if(iray.ne.0) p(i,j,4) = p(i,j,2)
240	call ggm (iray,p(i,j,2),p(i,j,1),p(i,j,4),p(i,j,3))
c
c	the distance between rows 1 and 5
	call distbear (p(nop,1,1),p(nop,1,2),p(nop,5,1),p(nop,5,2),
     a                 dist,bear)
	write(1,1112) dist
1112	format(' lateral distance between rows 1 and 5 is',f7.0,' km')
c
	jray = 2 * iray
	do 1111 ll = 1,nop
1111	type*,' lats & lons ',(p(ll,l,jray+1),p(ll,l,jray+2),l=1,5)
	type*, ' lateral distance between rows 1 and 5 is ',dist
c
c	set up the IRI parameters for each of the 5*NOP locations
	do 2001 ii = 1,nop
	do 2001 jj = 1,5
	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	F2-layer parameters - fof2 and m(3000)f2
	call fof2m3_ips
     a       (lunmap,mon,uthour,flat,flon,tindex,fof2,m3000)
c
c**********************************   can't do this for IRI
c	keep a pretend F1 layer to avoid a very wide E-F step when the
c	calculated value of foF1 is zero
ccc	if(fof1.eq.0.) fof1 = 1.4 * foe
c	make sure that fof2 exceeds fof1
ccc	if(fof1.ge.fof2) fof1 = (foe + fof2) / 2.
c**********************************
c
c	calculate the IRI profile parameters at each location
	axhi = -100.0
c	hmaxf2 is a supplied value (optional)
	hmaxf2 = 0.
	call irinhpro
     a  (flat,flon,0,r12,lmon,lday,hourlt,fof2,m3000,axhi,hmaxf2)
c
c	store the IRI parameters
	p(ii,jj,5) = hmf2
	p(ii,jj,6) = nmf2
	p(ii,jj,7) = hmf1
	p(ii,jj,8) = b0
	p(ii,jj,9) = b1
	p(ii,jj,10) = c1
	p(ii,jj,11) = hz
	p(ii,jj,12) = t
	p(ii,jj,13) = hst
	p(ii,jj,14) = hdx
	p(ii,jj,15) = hme
	p(ii,jj,16) = nme
	p(ii,jj,17) = hmd
	p(ii,jj,18) = nmd
	p(ii,jj,19) = hef
	p(ii,jj,20) = d1
	p(ii,jj,21) = xkk
	p(ii,jj,22) = fp30
	p(ii,jj,23) = fp3u
	p(ii,jj,24) = fp1
	p(ii,jj,25) = fp2
	p(ii,jj,26) = 0.
	if (night) p(ii,jj,26) = 1.
	p(ii,jj,27) = e(1)
	p(ii,jj,28) = e(2)
	p(ii,jj,29) = e(3)
	p(ii,jj,30) = e(4)
c
c	if requested, calculate the N(h) profiles to look at
	if(ipro.eq.0) go to 250
c	the plasma freq will be determined for 80(1)350 km, for the first
c	npro sets of 5 points
	npro = 3
	if(nop.lt.npro) npro = nop
	if(ii.gt.npro) go to 250
	ihs = 80
	ihf = 350
	do 180 ih = ihs,ihf
	iih = ih - ihs + 1
	h = ih
	ed = xe(h)
	fni = sqrt(ed/1.24e10)
180	fn(ii,jj,iih) = fni
c
250	continue
c
2001	continue
c
c	write out the profiles
	if(ipro.eq.0) go to 300
	write(4,400)
400	format(//5x,' First 3 * 5 profiles from 80 to 350 km'/)
c	latitudes
	write(4,401) ((p(l,j,1),j=1,5),l=1,npro)
c	longitudes
	write(4,401) ((p(l,j,2),j=1,5),l=1,npro)
	do 2000	 ih = ihs,ihf
	iih = ih - ihs + 1
	write(4,401) ((fn(l,j,iih),j=1,5),l=1,npro)
	if(mod(ih,10).eq.0) write(4,401) float(ih + 1)
2000	continue
401	format(16f5.1)
300	continue
c
c	convert the latitudes to co-latitudes
	do 280 i = 1,nop
	do 280 j = 1,5
	p(i,j,1) = 90. - p(i,j,1)
280	p(i,j,3) = 90. - p(i,j,3)
c
c	write out the P array
	write(1,332) nop
332	format(/i3,' profiles, including TX and RX')
	write(1,335) iray
335	format(i3,'  iray = 0 for geogc grid; =1 for geomc grid')
	ispec = ifix(special+0.5)
	write(1,336) ispec
336	format(i3,' special bearing case')
	do 333 i = 1,nop
	do 333 j = 1,5
333	write(1,334) (p(i,j,k),k=1,30)
334	format(5e16.8)
c
	type*,'  Length of circuit & bearing ',distsv,bearsv
	type*,'  Special bearing case ',ispec
c
	go to 10
	end