program QP_ICED
c
c	prog to set up a QP ionospheric model for raytracing with the 
c	Jones program. This is a version of ION_QP, which includes the
c	ICED mods to the high latitude trough and auroral zone.
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 13-aug-90
c	new ICED mods included 5-nov-90
c	modified 28-may-91 to change ymF2 calculation from MSIS83 method
c	to the use of normalized Millstone Hill profiles. The code referring
c	to MSIS is left in place (you never know). The main change is dropping
c	MSIS83 from the LINK command, and replacing LAYERS by LAYERS_MILL.
c
c-----	last modifed 27-sep-91
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 command is @L_ICED
c	$type L_iced.com
c	$pu
c	$link QP_ICED,ips,layers_mill,ilios,spoints,erectnh,profile,gg_gm,-
c	f2_trough,cnvcrd,eqboun,qoval1,iced_elayer,solpos,flag,feccir,-
c	smooth
c	$pu
c
c	ips		    loadgp,f2m3,tcon,gyro 
c	layers_mill	    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       ICED ->         f2_trough,cnvcrd,eqboun,qoval1,
c			iced_elayer,solpos,flag,feccir
c	f2_trough   [modified version of ICED s/r FLAYR]
c	cnvcrd	    etc are ICED s/rs (untouched)
c	smooth	    smooft,realft,four1 (for smoothing the characteristics)
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	gridpoint maps of gyrofreq, dip & modified dip
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	eff_q		effective Q index - for ICED
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,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
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	k=20	iflag (ICED parameter)
c
c	labelled common required by ICED subroutines
	common / f6parm / iyr,imo,ida,time,dayjul,effssn,effq,effkp,
     a                    gglat,gglon,cglat,cglon,tcgm
	common / coords / cg(2,36,89)
	common / solarc / iitime,ttime,asmlon,eeffssn,eeffkp,eeffq,
     a                    iimo,iida,iiyr
c
	dimension p(50,5,20),ap(7),nday(12),ia(15),item(50,5,3)
	dimension ht(500),fn(500),irgp(50),jtem(50,5,3),tem(1024)
	dimension grid_colat(100),grid_lon(50),param(20),ltem(50,5,2)
	dimension xhi(50,5,3)
	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	The TEM files are used for output of relevant parameters to file
c	for later perusal
c	ITEM	1   Corrected geomagnetic local time
c		2   Ratio of trough-modified value of fof2 to non-modified
c		2   Ratio of trough-modified value of hmf2 to non-modified
c	JTEM	1   10 * foE
c		2   10 * foF1
c		3   10 * foF2
c	LTEM	1   10 * foF1
c		2   ICED FLAG
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-----------	  Some assumptions    --------------------------------
c	number of points to smooth the characteristics over (along the
c	circuit) - this is just a guess. We cannot smooth too heavily, or 
c	the trough will disappear. Too little smoothing will yield errors
c	in the raytracing when the ray hits sudden boundaries.
	n_smooth = 5
c	The minimum value that the ICED value of foE is allowed to have
c	This is anybody's guess
	foe_iced_min = 0.654321
c-----------	  Some assumptions    --------------------------------
c
c	hardwire the year - used only for zenith angle calculations (trivial)
	lyear = 91
c	radius of Earth - consistent for all QP s/rs
	rz = 6370.
c	geographic coords of geomagnetic north pole - consistent with IRI
	gmpole_lat = 78.6
	gmpole_lon = 290.2
c	number of profile parameters at each grid point
	nparam = 20
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,20
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
	write(4,701) month,iday,ihour,minute
701	format(5x,' QP-ICED 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)
c	Read in the effective Q index (for high latitude work with ICED)
	print*,' enter effective Q index '
	read(5,501) eff_q
	write(7,702) r12,tindex,eff_q
	write(4,702) r12,tindex,eff_q
702	format(5x,' INDICES ', 3f10.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
	write(4,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
	write(4,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
	write(4,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
	call distbear(tlat,tlon,rlat,rlon,range,bearing)
	print*,' Enter number of gridpoints to insert between the'
	print*,' TX and the end of the path - equally spaced in '
	print*,' GEOTRANS latitude'
	print*,' '
	print*,' CIRCUIT LENGTH is',range
	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 for foE
	flux12 = 63.7 + 0.728 * r12 + 0.00089 * r12 * r12
c	use default flux values if not supplied (set to 0.) - used by MSIS
	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	OVERWRITE given location
	tlat = tlat_mod
	tlon = 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
	mon =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 - NOTE the fudge required to get the !!!!!!!!!!!!!
c	correct bearing when the TX lat = 90  (use 89.99)        !!!!!!!!!!!!!
c	TX is always at the north pole of GEOTRANS system
	call distbear (89.99,0.,rmlat,rmlon,range,bearing)
	print*,' GEOTRANS range & bearing ', 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 (89.99,0.,range*fac,bearing,end_lat,end_lon)
	print*,' GEOTRANS coords - new end of circuit', 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	the ranges of the gridpoints from the Tx
	do 215 i = 1,nlat
	call distbear
     a  (89.99,0.0,90.0-grid_colat(i),grid_lon(i),rgp,bear)
215	irgp(i) = rgp + 0.5
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
	print*,' --------- ii,jj --------- ', ii,jj
	flat = p(ii,jj,1)
	flon = p(ii,jj,2)
c
c	calculate the local time (from the UT) 
	time = uthour + flon/15.
	ut_time = time
	lh = time
	lm = (time - lh) * 60 + 0.5
ccc	item(ii,jj,1) = 100 * lh + lm
ccc	if(item(ii,jj,1).ge.2400) item(ii,jj,1) = item(ii,jj,1) - 2400
	lday = iday
	if(lh.ge.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.ge.24) lh = lh - 24
	reflon = flon
	hourlt = float(lh) + float(lm)/60.
	if(hourlt.gt.24.0) hourlt = hourlt - 24.
c
c	The ICED model also requires the parameters passed via /F6PARMS/,
c	/SOLARC/ and /COORDS/.
c	We need to set some of these up.
	iyr = lyear
	imo = month
	ida = iday
	dayjul = idayno
	effssn = r12
	gglat = flat
	gglon = flon
	effq = eff_q
c	effective value of Kp
	effkp = effq - 2.0
	if(effq.lt.3.0) effkp = effq / 3.0
	iitime = 100 * ihour + minute
	ttime = ut_time
	eeffssn = r12
	eeffq = eff_q
	eeffkp = effkp
	iimo = month
	iida = iday
	iiyr = lyear
c
c	corrected geomagnetic coords of this gridpoint
	call cnvcrd (gglat,gglon,0.,4,cglat,cglon,rdum,1,ierr)
c
c	corrected geomagnetic local time and the ICED flag value
c	work via the anti-subsolar point
	call solar 
     a  (lyear,lmon,lday,lh,lm,flat,flon,reflon,dec,zen,r,s)
	sslat =  dec
ccc	print*,' Declination ', dec
	sslon = 180. - 15.0 * uthour
	if(sslon.lt.0.) sslon = sslon + 360.
	asslat = -sslat
	asslon = sslon + 180.
	if(asslon.ge.360.) asslon = asslon -360.
	call cnvcrd (asslat,asslon,0.,4,asmlat,asmlon,rdum,1,ierr)
	call flag (gglat,gglon,cglat,cglon,iflag,tcgm)
ccc	print*,' cglat,cglon,asmlon,tcgm ', cglat,cglon,asmlon,tcgm
	item(ii,jj,1) = tcgm * 100
	ltem(ii,jj,2) = iflag
c
c	E-layer parameters for solar E layer - use foE for Dudeney 
c	calculation of hmF2
	call Elayer 
     a  (flux12,r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,
     a   chi,chi_noon,foe_solar,yme,hme)
ccc	print*,
ccc     a  flux12,r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,
ccc     a   chi,chi_noon,foe_solar,yme,hme
	type*,' E layer ', foe_solar,yme,hme
c
c	ICED E-layer parameters
	call iced_elayer(iflag,foe_iced,hme_iced)
c	ymE for this peak, assuming fn = 0.5 at 90 km (pretty good)
	rme = rz + hme_iced
	r = rz + 90.0
	if(foe_iced.le.0.5) print*,' F2_trough - foe_iced ', foe_iced
c	overwrite foe_iced to be foe_iced_min MHz (anybody's guess)
	if(foe_iced.le.foe_iced_min) foe_iced = foe_iced_min
	fr = 0.5 / foe_iced
	fac = sqrt(1.-fr*fr)
	fac = fac / (rme/r - 1.) + 1.
	yme_iced = rme / fac
c
ccccccccccccccccc	foe = foe_iced
ccccccccccccccccc	hme = hme_iced
	print*,' ICED E layer',foe_iced,hme_iced,yme_iced
c
c	F1-layer parameters
	call F1layer (lundip,r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,
     a  chi_max,fof1,ymf1,hmf1)
c
c	save the zenith angles
	xhi(ii,jj,1) = chi
	xhi(ii,jj,2) = chi_noon
	xhi(ii,jj,3) = chi_max
c
c	F2-layer parameters
	call F2layer (lunmap,month,uthour,hourlt,idayno,flat,flon,
     a  foe_solar,tindex,f107a,f107,ap,fof2,ymf2,hmf2)
c
	print*, ' Zenith angles chi, chi_max ', chi,chi_max
	print*,' F1 layer ',fof1,ymf1,hmf1
	print*,' F2 layer ', fof2,ymf2,hmf2,hourlt
c
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c	ICED modifications to foF2 and hmF2 in the high latitude regions
c=======================================================================
c
c	Modify fof2 and hmf2 according to the ICED model.
c	The calling sequence for this s/r requires everything used to calculate
c	fof2 and M(3000)F2 at the points of an interpolation grid.
c
	fof2_before = fof2
	hmf2_before = hmf2
	call f2_trough 
     a  (lunmap,month,uthour,hourlt,idayno,flat,flon,foe_solar,
     a  tindex,f107a,f107,ap,fof2,ymf2,hmf2,iflag)
	item(ii,jj,2) = (fof2 / fof2_before) * 100 + 0.5
	item(ii,jj,3) = (hmf2 / hmf2_before) * 100 + 0.5
ccc	if(fof2.ne.fof2_before.or.hmf2.ne.hmf2_before)
ccc     a  print*,' Trough modified ==>> fof2,ymf2,hmf2 ', fof2,ymf2,hmf2
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c	ICED modifications to foF2 and hmF2 in the high latitude regions
c=======================================================================
c
c
c	store the layer parameters in the array P
	p(ii,jj,5) = foe_iced
	p(ii,jj,6) = yme_iced
	p(ii,jj,7) = hme_iced
c
c
c	solar values - TEMPORARY FIX
C
	p(ii,jj,5) = foe_solar
	p(ii,jj,6) = yme
	p(ii,jj,7) = hme
c
c
c
c
	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
2001	continue
c
c	Now smooth the values of the profile parameters along the circuit
c	do parameters 5 to 13, skipping 9 & 10 (ymf1 & hmf1 - not used May 91)
	do 2005 ll = 5,13
	if(ll.eq.9) go to 2005
	if(ll.eq.10) go to 2005
	do 2004 jj = 1,nlon
	do 2002 ii = 1,nlat
2002	tem(ii) = p(ii,jj,ll)
	IF(jj.eq.1.and.LL.EQ.5) print*, (tem(i),i=1,nlat)
	call smooft (tem,nlat,float(n_smooth))
	if(jj.eq.1.and.ll.eq.5) print*,(tem(i),i=1,nlat)
	do 2003 ii = 1,nlat
c	we need this little check in case fof1 (especially) goes negative
	if(tem(ii).lt.0.) tem(ii) = 0.
2003	p(ii,jj,ll) = tem(ii)
	if(jj.eq.1.and.ll.eq.5) print*,(tem(i),i=1,nlat)
2004	continue
2005	continue
c
c	we are going to print out 10 * foe, 10 * fof2 & hmf2
	do 2006 ii = 1,nlat
	do 2006 jj = 1,nlon
	jtem(ii,jj,1) = 10. * p(ii,jj,5) + 0.5
	jtem(ii,jj,2) = 10. * p(ii,jj,11) + 0.5
	jtem(ii,jj,3) = p(ii,jj,13) + 0.5
c	also fof1
	ltem(ii,jj,1) = 10. * p(ii,jj,8) + 0.5
2006	continue
c
c	calculate the QP profile parameters at each location
	do 2020 ii = 1,nlat
	do 2020 jj = 1,nlon
	fof2 = p(ii,jj,11)
	foe_iced = p(ii,jj,5)
	yme_iced = p(ii,jj,6)
	hme_iced = p(ii,jj,7)
c	need to watch for the relative value of foE and foF1/foF2
	if(foe_iced.lt.fof2 - 0.1) go to 300
c	take the profile to be just an E layer
	hsf = hme_iced
	hsf2 = hme_iced
	fnmax = 0.
	ymax = 0.
	hmax = hme_iced
	invqp = -1	
	p(ii,jj,11) = foe_iced
	p(ii,jj,12) = yme_iced
	p(ii,jj,13) = hme_iced
c	the ray should now penetrate at hmF2 = hmE
	go to 320
c
300	continue
	foe = p(ii,jj,5)
	yme = p(ii,jj,6)
	hme = p(ii,jj,7)
	fof1 = p(ii,jj,8)
	ymf1 = p(ii,jj,9)
	hmf1 = p(ii,jj,10)
	fof2 = p(ii,jj,11)
	ymf2 = p(ii,jj,12)
	hmf2 = p(ii,jj,13)
	hsf = p(ii,jj,14)
	hsf2 = p(ii,jj,15)
	fnmax = p(ii,jj,16)
	ymax = p(ii,jj,17)
	hmax = p(ii,jj,18)
	invqp = p(ii,jj,19)
	chi      = xhi(ii,jj,1)
	chi_noon = xhi(ii,jj,2)
	chi_max  = xhi(ii,jj,3)
ccc	print*,' II & JJ ', ii,jj
c	throw away the smoothed F1 layer if foe > fof1
	if(fof1.le.foe + 0.1) fof1 = 0.
c
	if(foe.gt.fof1.and.fof1.gt.0) print*,' Critical freqs E & F1! ',
     a  ii,jj,
     a  rz,foe,yme,hme,fof1,ymf1,hmf1,fof2,ymf2,hmf2,chi,chi_noon,
     a   chi_max,hsf,hsf2,fnmax,hmax,ymax,invqp
	if(foe.gt.fof2) print*,' Critical freqs E & F2! ',
     a  ii,jj,
     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	calculate the remaining profile parameters hsf etc
	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
320	continue
	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
	p(ii,jj,20) = iflag
c
c	calculate profile at TX, & write to file
	if(ii.eq.1.and.jj.eq.1) go to 1998
	go to 2020
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.and.mod(i,10).eq.0) write(4,400) i,ht(i),fn(i)
	if(fn(i).gt.0.and.mod(i,10).eq.0) write(6,400) i,ht(i),fn(i)
1999	continue
400	format(i4,f6.1,f6.2)
c
2020	continue
c
c	write the parameters of the model to the output file
	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
c	Write useful stuff to the generic output file for perusal
	write(4,9999)
	write(4,750)
750	format(5x,'Co-latitudes of gridpoints')
	write(4,7000) (grid_colat(i),i=1,nlat)
7000	format(10f8.2)
	write(4,751)
751	format(5x,'Longitudes of gridpoints')
	write(4,7000) (grid_lon(i),i=1,nlon)
c
c	write out lat/lon tales of 10 * foe, fof2 & hmf2
c	can do up to 5 longitudes
	write(4,9999)
	write(4,752)
752	format(5x,'10 * foE, 10 * foF2 & hmF2'/)
	mlon = nlon
	if(mlon.gt.5) mlon = 5
	do 560 i =1,nlat
	do 555 k = 1,15
555	ia (k) = 0
	do 556 j = 1,mlon
	ia(j) = jtem(i,j,1)
	ia(j+5) = jtem(i,j,2)
	ia(j+10) = jtem(i,j,3)
556	continue
	write(6,557) irgp(i),(ia(k),k=1,15)
	write(4,557) irgp(i),(ia(k),k=1,15)
557	format(i5,3(5x,5i4))
560	continue
c
c	write out lat/lon tables of TCGM, fof2 ratio & hmf2 ratio
c	can do up to 5 longitudes
	write(4,9999)
9999	format(//)
	write(4,753)
753	format(5x,'Corrected geomag time, foF2 ratio & hmF2 ratio'/)
	mlon = nlon
	if(mlon.gt.5) mlon = 5
	do 1560 i =1,nlat
	do 1555 k = 1,15
1555	ia (k) = 0
	do 1556 j = 1,mlon
	ia(j) = item(i,j,1) 
	ia(j+5) = item(i,j,2) 
	ia(j+10) = item(i,j,3)
1556	continue
	write(6,557) irgp(i),(ia(k),k=1,15)
	write(4,557) irgp(i),(ia(k),k=1,15)
1560	continue
c
c	write out fof1 and FLAG
	write(4,9999)
	write(4,754)
754	format(5x,'Values of 10*foF1 and the FLAG values (Correct?)'/)
	mlon = nlon
	if(mlon.gt.5) mlon = 5
	do 1570 i =1,nlat
	do 1565 k = 1,15
1565	ia (k) = 0
	do 1566 j = 1,mlon
	ia(j) = ltem(i,j,1)
	ia(j+5) = ltem(i,j,2)
1566	continue
	write(6,557) irgp(i),(ia(k),k=1,10)
	write(4,557) irgp(i),(ia(k),k=1,10)
1570	continue
	go to 10
c
	end