C	under name : 'profgentilt.for'
	 subroutine profgentilt
c     a	  (lyear,mon,iday,ih,im,tindex,r12,f107a,f107,ap,flat,flon)
c
c	s/r to set up a tilted Dudeney profile for use by the JONES program
c	written by leo mcnamara 23-mar-87
c-----	last modified 10-jul-87
c
c	note that this is all single precision, except for hooks back
c	into the calling program
c........................................................................
c
c---	Requires files ips,layers,dudpro,dudsub,fixpntSP,msis83,ilios,ggm
c	ips	loadgp,f2m3,tcon,gyro (single precision)
c	layers	Elayer,solar,foeccir
c		F1layer,moddip,gdmd_interp,fof1OT
c		F2layer	
c	dudpro
c	dudsub	dimsol,hmax,ymax
c	interp	(attached here)
c	fixpntsp
c	msis83	MSIS83 s/rs GTS4 etc
c	ilios	
c	ggm
c
c---	required input data (from DUD.DAT):-
c	lyear	year (local time !)
c	mon	month
c	iday	UT day
c	ih	UT hour
c	im	UT minute
c	tindex	ionospheric index, T
c	r12	12-month smoothed sunspot number
c	f107a	3-month average 10.7cm flux (for MSIS83)
c	f107	10.7cm flux for previous day
c	ap	AP index - daily or array (See GTS4 listing)
c	ffof2	optional - directly from ionogram
c	m3000	optional - directly from ionogram
c	hff2	optional - directly from ionogram
c	flat	geographic latitude of SKYLOC (negative for south)
c	flon	EAST geographic longitude of SKYLOC
c	hapog	height corresponding to a tilt measurement 
c		= 0 for short & medium-range circuits
c		= apogee for untilted case for long-range circuits
c	tilt    the tilt of the ionosphere [at a height of hapog for LR]
c	ga	azimuth of tilt [set to az of ray if hapog.ne.0]
c
c---	output data
c	x,pxpr etc in labelled common /XX/
c
	real m3000
c-------------------------------Jones is real*8-----------------
	real*8 x,pxpr,r,w0,w,ff,hh,pxpth,pxpph,pxpt,hmax
	real nlon,nlat
	character*8 modx,id
	common r(6) /ww/ id(10),w0,w(400)
	common /xx/ modx(2),x,pxpr,pxpth,pxpph,pxpt,hmax
c------------------ 3-jul-87 tilt additions --------------------
	common /ioncent/ hapog, tilt, cosba, sinba,
     &		         cosblon, sinblon, cosblat, sinblat
c---------------------------------------------------------------
	common / iondat / fof2,fof1,foe,hff2,m3000,hme,yme,hmf2,
     a	ymf2,h0,xd,dv
	character*10 munth(12)
	dimension nday(12),ap(7)
	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 /
	data nday / 31,28,31,30,31,30,31,31,30,31,30,31 /
	data ipass / 0 /
c
c	conversion from degrees to radians = d2r = pi/180
	data d2r, pi / 1.745329252e-2, 3.1415926536 /
	data big / 1.0e+10 /
c
c
c......................
	ENTRY ELECTX
c......................
c
c	bypass determination of profile parameters on second and subsequent
c	passes
ccc	type*,'rad,co-lat,long=',r(1),r(2)/d2r,r(3)/d2r
	if(ipass.eq.1) go to 100
c
c	lundip	logical unit number for grid-point maps of MODDIP etc
	lundip = 9
	open(unit=9,file='gyro_gp_data.dat',status='old')
c	required input is from DUD.DAT
	open (unit=3,file='dud.dat',status='old')
c
10	continue
ccc	write(6,101)
101	format(//' enter UT year,month,day,hour,min  --  yymmddhhmm'/)
	read(3,200) lyear,mon,iday,ih,im
	write(2,500) lyear,mon,iday,ih,im
500	format(/2x,' yymmddhhmm ',5i2)
200	format(5i2)
	if(lyear.eq.0) stop ' all done'
c
c	lunmap	logical unit number for IPS grid-point maps of fof2 & M(3000)
	lunmap = 8
cc	open(unit=8,file='xgp.mar',status='old')
	open (unit=lunmap,file=munth(mon),status='old')
c
ccc	write(6,102)
102	format(//' enter Tindex & R12  --  2f5.0'/)
	read(3,201) tindex,r12
	write(2,501) tindex,r12
501	format(2x,' t index & r12 ',2f5.0)
ccc	write(6,106)
106	format(//' enter hFF2, the minimum value of h"F f5.0'/,
     a  '  PRESS RETURN if no value available'/)
	read(3,201) ffof2,m3000,hff2
	write(2,502) ffof2,m3000,hff2
502	format(2x,' ffof2,m3000,hff2 ',2f6.2,f6.0)
ccc	if(hff2.gt.0.) go to 107
ccc	write(6,104)
104	format(//' enter f107a,f107 & Ap  --  3f5.0'/)
	read(3,201) f107a,f107,ap(1)
c	ccir formula relating r12 and flux12 (supp 252)
	flux12 = 63.7 + 0.728 * r12 + 0.00089 * r12 * r12
c	allow for 3-month & daily fluxes not being specified - derive
c	values from effective sunspot number
	if(f107a.eq.0.) f107a = flux12
	if(f107 .eq.0.) f107  = flux12
	write(2,503) f107a,f107,ap(1)
503	format(2x,' f107a,f107,ap(1) ',3f5.1)
ccc107	write(6,105)
105	format(//' enter lat & lon of reflection point - 2f5.0'/)
c	glat & glon give the position of the beacon site at which the
c	tilt measurement is made. This site is taken as the location
c	of the ionospheric reflection point.
	read(3,201)glat,glon
	write(2,504) glat,glon
504	format(2x,' location of beacon & reflection point ',2f5.1)
c	put the reflection point at glat,glon
	flat = glat
	flon = glon	
201	format(4f5.0)
ccc	write(6,1001)
1001    format(//' enter apogee height, tilt and ',
     &  		'azimuth of apogee layer- 3f5.2'/)
	read(3,1002) tilt, ga,hapog
c	check that w(42) is not too small [we have a few accuracy problems]
	if(tilt.ne.0.0) w(42) = dmax1(w(42),0.001d0)
c	for long-range circuits (defined as having hapog.ne.0) the azimuth
c	of the tilt is set to the first azimuth of the ray
	if(hapog.gt.0.) ga = sngl(w(11)) * 180.0 / pi
	write(2,506) hapog,tilt,ga,w(42)
506	format(2x,' hapog,tilt,azim of tilt,W(42) ',3f5.1,f7.4/)
1002	format(3f5.2)
c
c	UT hour
	uthour = float(ih) + float(im)/60.
c	day of year (for MSIS83)
	idayno = iday
	if(mon.eq.1) go to 20
	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
20	continue
c
c	calculate the local time (from UT)
	ndm = nday (mon)
	if(mod(lyear,4).eq.0.and.mon.eq.2) ndm = 29
	lmon = mon
	time = ih + flon/15. + float(im)/60.
	lh = time
	lm = (time - lh) * 60 + 0.5
	lday = iday
	if(lh.ge.24.and.flon.ge.180.) lday = lday - 1
	if(lh.ge.24.and.flon.le.180.) lday = lday + 1
	if(lday.gt.ndm) lmon = mon + 1
	if(lday.gt.ndm) lday = 1
	if(lmon.eq.13) lmon = 1
	if(lh.ge.24) lh = lh - 24
	if(lday.eq.0) lmon = lmon - 1
	ndmm = nday(lmon)
	if(mod(lyear,4).eq.0.and.lmon.eq.2) ndmm = 29
	if(lday.eq.0) lday = ndmm
	if(lmon.eq.0) lmon = 12
	reflon = flon
	type*,' day no. & local time (ddhhmm)',idayno,lday,lh,lm
	hourlt = float(lh) + float(lm)/60.
c
c	E-layer parameters 
	call Elayer 
     a  (flux12,r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,foe,yme,hme)
	type*,' E layer ', foe,yme,hme
c	make the D layer tangential to the E layer at 0.5MHz
	xd = 0.5 / foe
c	we don't need a D layer at night 
	if(foe.lt.1.0) xd = 0.
c
c	F1-layer parameters
	call F1layer ( lundip,
     a  r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,fof1,ymf1,hmf1)
	type*,' F1 layer ',fof1,ymf1,hmf1
c	can't have E-F valley AND F1 layer 
	dv = 0.1
	if(fof1.ne.0.) dv = 0.
c
c	F2-layer parameters
	if(ffof2.eq.0.0.or.m3000.eq.0.0.or.hff2.eq.0.0)
     a	call F2layer
     a  (lunmap,mon,uthour,hourlt,idayno,flat,flon,foe,tindex,f107a,
     a   f107,ap,fof2,ymf2,hmf2)
c	use scaled value of fof2 if provided (non-zero)
	if(ffof2.gt.0.) fof2 = ffof2
c	use scaled values of hff2 and m3000 to get hmf2 & ymf2 in DUDSUB
	if(m3000.gt.0.) hmf2 = 0.
	if(hff2.gt.0.) ymf2 = 0.
	type*,' F2 layer ', fof2,ymf2,hmf2
c
c	base of layer - 80km for day;set by DIMSOL to hme-yme if h0=0
	h0 = 80.
	if(xd.le.0.) h0 = 0.
c
c-------------------------------------------------------------------------------
c       3-jul-87 (david barker)
c	this section calculates the centre of the tilted ionosphere
c	in cartesian co-ords (xion,yion,zion) wrt GEOM co-ords
c       and also converts GEOG to GEOM co-ords in the case where we have
c	a field (the usual case)
c
	ga = d2r*ga
c
	if (w(1).ne.0.0) then
c
	call GGM(0.,glon,glat,blon,blat)
ccc	blon = blon - 360.
ccc	type*,' glon,glat,blon,blat ',glon,glat,blon,blat
	blon = d2r*blon
	blat = d2r*blat
c
	call GGM(1.,nlon,nlat,0.0,90.0)
ccc	type*,'nlat,nlon=',nlat,nlon
	nlon = d2r*nlon
	nlat = d2r*nlat
ccc	type*,'glat,glon=',glat,glon
	glon = d2r*glon
	glat = d2r*glat
c
	del = glon-nlon
	if (del.gt.pi) del = del -2*pi
	if (del.le.-pi)  del = del +2*pi
ccc	type*,'del=',del/d2r
	top = cos(nlat)*cos(glat)*sin(del)
	down = sin(nlat)*sin(glat)+cos(nlat)*cos(glat)*sin(del) 
	down = sin(nlat) -sin(glat)*down
	if (abs(top).lt.(big*abs(down))) then
	    tang = top/down
	else
	    tang = sign(1.0,top)*sign(1.0,down)*big
	endif
ccc	type*,'top,down=',top,down
ccc	type*,'tang=',tang
	g = atan(tang)
	if (del.ge.0.0) then
	    if (g.lt.0.0) g = g +pi
	else
	    if (g.gt.0.0) g = g -pi
	endif
	dela = -g
c
c
	else
c
c
	blat = d2r*glat
	blon = d2r*glon
	dela = 0.0
c
c
	endif
ccc	type*,'dela=',dela/d2r
	ba = ga + dela
ccc	type*,'blat,blon,ba,tilt=',blat,blon,ba/d2r,tilt
c
	tilt = d2r*tilt
	cosba = cos(ba)
	sinba = sin(ba)
	cosblon = cos(blon)
	sinblon = sin(blon)
	cosblat = cos(blat)
	sinblat = sin(blat)
c
c-------------------------------------------------------------------------------
c
100	continue
	ff = w(6)
	call dudtilt(ipass,ff)
c	the values of ymf2 & hmf2 may have been set up in DIMSOL
	if(ipass.eq.0) type*,' ymF2,hmF2 ',ymf2,hmf2
	ipass = 1
	return
	end