program ionmod_wdc
c
c	prog to determine WDC profile parameters at a grid of
c	locations, for use in 3-D raytracing using the Jones program
c
c	This profile consists of 2 or 3 chapman layers for the E, (F1) &
c	F2 layers; a full E-F valley ; and elliptical sections to smooth
c	the corners at the E-F1/F2 and F1-F2 intersections
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 31-jul-87 (from ionmod_dud)
c---	last modified 5-oct-87
c
c	LINK with files ips,layers,MSIS83,ilios,spoints,chap,
c	               wdcpro,ellipse,elled,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	chap        chap
c	wdcpro
c	ellipse
c	elled
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	ltype   F2 layer type =0 for alpha Chapman, 1 for beta Chapman
c
c---	output data  .. in the file IONMOD_WDC.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,25  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	she
c	k=7	hme
c	k=8	foF1 F1 layer
c	k=9	shf1
c	k=10    hmf1
c	k=11    foF2 F2 layer
c	k=12	shf2 [ <0 for beta Chapman layer ]
c	k=13	hmf2
c	k=14	h0 - base of ionosphere
c	k=15	xce   E - F1/F2 ellipse
c	k=16	yce
c	k=17	ae
c	k=18	be
c	k=19	h1
c	k=20	xcf   F1 - F2 ellipse
c	k=21	ycf
c	k=22	af
c	k=23	bf
c	k=24	h2
c	k=25    width
c
	dimension p(50,5,25),ap(7),flatt(50),flonn(50),nday(12)
	dimension fn(3,5,275),shh(50),en(50)
	character*10 munth(12)
	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 = 15.0
c	lat and lon spacing of grid for special cases
	stplon = 2.0
	stplat = 2.0
c
c	width of E-F1 valley (modified later)
	vwidth = 20.
c
c	types of chapman layer (beta may be OVERWRITTEN later)
	alpha = 0.5
	beta = 1.0
c
c	smoothing ellipses are fitted to the F1 or F2 layers at a height
c	HSTEP above the peak of the underlying layer
	hstep = 2.
c
c	base of ionosphere
	base = 80.
c
c	size of array dimension
	ndim = 50
c	set the year in concrete
	lyear = 87
c
c	output is to file IONMOD_WDC.DAT
	open (unit=1,file='ionmod_wdc.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,25
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 WDC 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,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
	write(6,113)
113	format(//' Enter 1 for beta Chapman layer, RETURN for alpha')
	read(5,104) ltype
	if(ltype.eq.0) beta = alpha
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	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
	call spoints (flat1,flon1,flat2,flon2,nip,flatt,flonn)
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 WDC 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,foe,yme,hme)
	enmaxe = 1.24e10 * foe * foe
	type*,' E  layer ', foe,yme,hme
c
c	We want the Chapman E layer to pass through the McNamara(1979)
c	value of h(0.5MHz), h5.  We can reconstruct this value from yme, 
c	and then interpolate to find the scale height SHE which gives 
c	fn = 0.5MHz at h5.  For details of the interpolation, see how shF1 
c	is calculated later.
	rbe = 6370. + hme - yme
	rme = 6370. + hme
	en5 = 1.24e10 * 0.5 * 0.5
c	use positive square root to get the lower solution
	rr = rme / (1. + yme/rbe * sqrt(1. - en5/enmaxe) )
	h5 = rr - 6370.
c	guesses at range of she
	is1 = (yme/2.) * 0.75
	is2 = (yme/2.) * 1.5
	iss = 1
	i = 0
	en5 = alog (en5)
	do 705 ish = is1,is2,iss
	i = i + 1
	shh(i) = float(ish)
	call chap (foe,hme,shh(i),alpha,h5,en(i),dndh)
	if(en(i).le.0.) en(i) = 1.e7 * i
	en(i) = alog (en(i))
	if(i.eq.1) go to 705
	if(en(i).gt.en5.and.en(i-1).le.en5) go to 706
705	continue
	stop ' Problems in calculating SHE'
706	nv = i
	ll = i
	if(ll.gt.4) ll = 4
	call ilios (shh,en,en5,nv,ll,she)
	yme = 2. * she
	type*, 'h(0.5) & modified SHE ',h5,she
c
c	F1-layer parameters
	call F1layer (lundip,
     a  r12,lyear,lmon,lday,lh,lm,flat,flon,reflon,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
	type*,' F1 layer ',fof1,ymf1,hmf1
	type*,' F2 layer ', fof2,ymf2,hmf2,hourlt
c
c	modify hmf1 so that the F2 layer has fn = fof1 at hmf1
c	(for pseudo F1 as well)
c	see calculation of shf1 below for more details of method
c	fof1 for pseudo F1 layer
	if(fof1.eq.0.) psfof1 = 1.4 * foe
	if(psfof1.gt.fof2) psfof1 = 0.5 * (foe + fof2)
	fof = fof1
	if(fof1.eq.0.) fof = psfof1
	ih2 = hmf2 - 20.
	ih1 = hme
	ihs = 5.
	i = 0
	enf1 = alog (fof * fof * 1.24e10)
	shf2 = ymf2 / 2.
	do 315 ihh = ih1,ih2,ihs
	i = i + 1
	hh = float(ihh)
	shh(i) = hh
	call chap (fof2,hmf2,shf2,beta,hh,en(i),dndh)
c	protect against zero/negative Ne
	if(en(i).le.0.) en(i) = 1.e7 * i
	en(i) = alog(en(i))
	if(i.eq.1) go to 315
	if(en(i).gt.enf1.and.en(i-1).le.enf1) go to 316
315	continue
	stop 'trouble in calculating hmF1'
316	nv = i
	ll = i
	if(ll.gt.4) ll = 4
	call ilios (shh,en,enf1,nv,ll,hmf1)
c	drop hmf1 by HSTEP to leave room to fit the smoothing ellipse
	hmf = hmf1
	if(fof1.ne.0.) hmf1 = hmf1 - hstep
	type*,' modified value of hmF1 ',hmf1
c
c	make the F1 chapman layer pass through [foe,hme + vwidth], by
c	modifying shf1
	if(fof1.eq.0.) go to 350
c	to do this we set up a table of shf1 vs Ne, and interpolate to
c	find the value of shf1 corresponding to NmaxE
	is1 = 20.
	is2 = 180.
	iss = 5.
	hh = hme + vwidth
	i = 0
c	need to use logs for interpolation
	enme = alog (enmaxe)
	do 305 ish = is1,is2,iss
	i = i + 1
	shh(i) = float(ish)
	call chap (fof1,hmf1,shh(i),alpha,hh,en(i),dndh)
c	protect against zero/negative Ne
	if(en(i).le.0.) en(i) = 1.e7 * i
	en(i) = alog (en(i))
ccc	type*,'i,shh,NmE,en(i) ', i,shh(i),enme,en(i)
c	see if en(i) and en(i-1) straddle nme
	if(i.eq.1) go to 305
	if(en(i).gt.enme.and.en(i-1).le.enme) go to 306
305	continue
c	prog should not reach here
	stop ' Problems in calculating shF1 '
306	nv = i
	ll = i
c	order of interpolating polynomial
	if(ll.gt.4) ll = 4
	call ilios (shh,en,enme,nv,ll,shf1)
	ymf1 = 2. * shf1
	type*,' modified shf1 ',shf1
350	continue
c
c	scale heights of layers (ym is defined as 2 * sh)
	she = yme / 2.
	shf1 = ymf1 / 2.
	shf2 = ymf2 / 2.
c
	if(fof1.ne.0.) go to 360
c
c	when there is no real f1 layer, represent the e-f transition by
c	an ellipse tangential to the e-layer peak, and to the f2 layer at 
c	the pseudo fof1
c	avoid some problems by insisting that hmf2 lies at least 50km
c	above hmf1 -- fof1 is ignored, and ellipse is tangential to the
c	F2 layer at hmff
	hmff = hmf
	if(hmf2-hmf.lt.50.) hmff = hmf2 - 50.
	if(hmf2-hmf.lt.50.) type*,' hmF1 reduced to ',hmff
	call chap (fof2,hmf2,shf2,beta,hmff,ed,dndh)
	enmaxe = 1.24e10 * foe * foe
	call ellipse (hme,enmaxe,hmff,ed,dndh,xce,yce,ae,be)
c	check for no real solution
	if(be.gt.0.) go to 355
	stop ' trouble fitting the E-F2 ellipse'
355	h1 = hmf
c	there is no f1-f2 ellipse
	xcf = 0.
	ycf = 0.
	af = 0.
	bf = 0.
	h2 = 0.
c	there is now no slab valley
	width = 0.
	go to 380
c
360	continue
c	the e-f1 ellipse is tangential to the f1 layer at the
c	height h1, and to the E-F1 slab valley at
c	the height hme + vwidth - hstep
	h1 = hme + vwidth + hstep
c	h1 will be decreased until a solution is found
	itermax = 10
	hm = hme + vwidth - hstep
	enmaxe = 1.24e10 * foe * foe
	iter = 0
365	call chap (fof1,hmf1,shf1,alpha,h1,ed,dndh)
	call ellipse (hm,enmaxe,h1,ed,dndh,xce,yce,ae,be)
c	no solution was found if be < 0
	if(be.gt.0.) go to 370
	iter = iter + 1
	h1 = h1 - hstep * 0.2
	if(iter.gt.itermax) stop ' TROUBLE in calculating h1'
	go to 365
370	continue
c	decrease the valley width to match the bottom of the fitted ellipse
	width = hm - hme
c
c	the f1-f2 ellipse is tangential to the f2 layer at the
c	height h2
	h2 = hmf1 + 2. * hstep
	enmaxf1 = 1.24e10 * fof1 * fof1
	iter = 0
375	call chap (fof2,hmf2,shf2,beta,h2,ed,dndh)
	call ellipse (hmf1,enmaxf1,h2,ed,dndh,xcf,ycf,af,bf)
	if(bf.gt.0.) go to 380
	h2 = h2 - hstep * 0.2
	iter = iter + 1
	if(iter.gt.10) 	stop ' TROUBLE calculating h2'
	go to 375
c
380	continue
c
c	store the layer parameters in the array P
	p(ii,jj,5) = foe
	p(ii,jj,6) = she
	p(ii,jj,7) = hme
	p(ii,jj,8) = fof1
	p(ii,jj,9) = shf1
	p(ii,jj,10) = hmf1
	p(ii,jj,11) = fof2
	p(ii,jj,12) = shf2
	if(ltype.eq.1) p(ii,jj,12) = -shf2
	p(ii,jj,13) = hmf2
	p(ii,jj,14) = base
	p(ii,jj,15) = xce
	p(ii,jj,16) = yce
	p(ii,jj,17) = ae
	p(ii,jj,18) = be
	p(ii,jj,19) = h1
	p(ii,jj,20) = xcf
	p(ii,jj,21) = ycf
	p(ii,jj,22) = af
	p(ii,jj,23) = bf
	p(ii,jj,24) = h2
	p(ii,jj,25) = width
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
	ihs = 80
	ihf = 350
	if(nop.lt.npro) npro = nop
	if(ii.gt.npro) go to 250
	do 180 ih = ihs,ihf
	iih = ih - ihs + 1
	h = ih
	call wdcpro
     a  (base,foe,hme,she,fof1,hmf1,shf1,fof2,hmf2,shf2,width,
     a   xce,yce,ae,be,h1,xcf,ycf,af,bf,h2,h,ed,dndh)
	fni = sqrt(ed/1.24e10)
	if(ii.eq.1.and.jj.eq.1) write(4,8765) h,fni,ed,dndh
8765	format(f5.0,f8.3,2e20.6)
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
ccc	type*,' ihs,ihf,ih,iih ',ihs,ihf,ih,iih
	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,25)
334	format(5e16.8)
c
	type*,'  Length of circuit & bearing ',distsv,bearsv
	type*,'  Special bearing case ',ispec
c
	go to 10
	end