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