program ionmod_dud
c
c	prog to determine Dudeney 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 29-may-87
c	modified 15-jul-91 to accept full unformatted maps of fof2
c---	last modified 15-jul-91
c
c	LINK with files ips,layers,MSIS83,ilios,spoints,dudsub,
c	                dud_jones,ggm
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	dudsub	    dimsol,[hmax,ymf2 - not used]
c	dud_jones   dudney [only for looking at the N(h) profile]
c	ggm
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	f107a	3-month average 10.7cm flux (for MSIS)
c	f107	10.7cm flux for previous day
c	ap	AP index - daily or array (See GTS4 listing)
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_DUD.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,20  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 (no need for ymf1 & hmf1)
c	k=9	foF2 F2 layer
c	k=10	ymf2
c	k=11	hmf2
c	k=12	h0 base of ionosphere -- DUDENEY parameters
c	k=13	ad
c	k=14	tau
c	k=15	hd 
c	k=16	beta
c	k=17	ht
c	k=18	ft
c	k=19	at
c
	dimension p(50,5,20),ap(7),flatt(50),flonn(50),nday(12)
	dimension fn(3,5,275)
	character*10 munth(12)
	character*8 modx
      LOGICAL LCOND
      REAL M3000
	real*8 X,freq,hh
      COMMON/IONDAT/FOF2, FOF1, FOE, HFF2, M3000, HME, YME,
     *              HMF2, YMF2, H0, XD, DV
      COMMON/PARMS/BETA, HT, FT, At, HF, AD, TAU, HD, LCOND
	common / XX / modx(2),X, plus, other, stuff
	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	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_DUD.DAT
	open (unit=1,file='ionmod_dud.dat',status='old')
c	grid-point data file on logical unit 2 - see later for OPEN
cc	open (unit=2,file=munth(mon),form='unformatted',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,20
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	the virtual height of the base of the F1/F2 layer is not used
	hff2 = 0.
c
	write(6,600)
600	format(//5x,'Prog to set up the profile parameters for',/
     a  '     a 3-D ionospheric model, using the DUDENEY 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),form='unformatted',
     a	  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,105)
105	format(//2x,'Enter 3-month average 10cm flux, previous days',
     a ' flux',/'   and present days daily Ap --  3F5.0')
	read(5,201) f107a,f107,ap(1)
	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	use default flux values if not supplied (set to 0.)
	if(f107.eq.0.) f107 = flux12
	if(f107a.eq.0.) f107a = flux12
c	radius of Earth
	rz = 6370.
	uthour = float(ih) + float(min)/60.
c	day number of year ( for MSIS )
	idayno = iday
	if(mon.eq.1) go to 12
	mon1 = mon - 1
	do 11 i = 1,mon1
11	idayno = idayno + nday(i)
	if(mod(lyear,4).eq.0.and.mon.ge.3) idayno = idayno + 1
12	continue
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 layer and Dudeney 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	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,mon,uthour,hourlt,idayno,flat,flon,
     a  foe,tindex,f107a,f107,ap,fof2,ymf2,hmf2)
c
c**********************************
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
	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) = fof2
	p(ii,jj,10) = ymf2
	p(ii,jj,11) = hmf2
c
c	calculate the Dudeney profile parameters at each location
c	base of layer
	h0 = 80.
c	xd  --  D and E layers touch at xd * foe
	xd = 0.5
	call DIMSOL (ierror)
	if(ierror.eq.0) go to 1200
	write(6,1150) ierror,flat,flon
1150	format(///5x,' *** error in DIMSOL ',i4, 'lat&lon =',2f6.2)
	stop
1200	continue
c
c	store the Dudeney parameters
	p(ii,jj,12) = h0
	p(ii,jj,13) = ad
	p(ii,jj,14) = tau
	p(ii,jj,15) = hd
	p(ii,jj,16) = beta
	p(ii,jj,17) = ht
	p(ii,jj,18) = ft
	p(ii,jj,19) = at
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
c	we use s/r DUDNEY in file DUD_JONES, which was not written for the 
c	present purposes
c	set frequency = 1., to get fn, not X
	freq = 1.d0
	ihs = 80
	ihf = 350
c	we don't want to set up the profile parameters again
	ipass = 1
	do 180 ih = ihs,ihf
	iih = ih - ihs + 1
	h = ih
	hh = dble (h)
	call dudney (ipass,freq,hh)
c	result is returned in labelled common / xx/
	fni = sngl(x)
	fni = sqrt(fni)
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,20)
334	format(5e16.8)
c
	type*,'  Length of circuit & bearing ',distsv,bearsv
	type*,'  Special bearing case ',ispec
c
	go to 10
	end