program ION_QP
c prog to set up a QP ionospheric model for raytracing with the
c Jones program.
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 31-jul-90
c----- last modified 10-aug-90
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 with files ips,layers,MSIS83,ilios,spoints,erectnh,
c profile,gg_gm
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 erectnh erectnh,ym,qphfn,qlhfn,qlfnh,qpfnh,parabola,efqp
c profile [only for looking at the N(h) profile]
c gg_gm
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
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 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,19 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
dimension p(50,5,19),ap(7),nday(12)
dimension ht(500),fn(500)
dimension grid_colat(100),grid_lon(50),param(19)
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 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 hardwire the year
lyear = 89
c radius of Earth - consistent for all QP s/rs
rz = 6370.
c
c geographic coords of geomagnetic north pole - consistent with IRI
gmpole_lat = 78.6
gmpole_lon = 290.2
c
c number of profile parameters at each grid point
nparam = 19
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,19
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
701 format(5x,' QP 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)
write(7,702) r12,tindex
702 format(5x,' INDICES ', 2f10.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
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
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
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
print*,' Enter number of gridpoints to insert between the'
print*,' TX and the end of the path - equally spaced in '
print*,' GEOTRANS latitude'
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 by MSIS
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
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
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
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
call distbear (0.,0.,rmlat,rmlon,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 (0.,0.,range*fac,bearing,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 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
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,month,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 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) = ymf1
p(ii,jj,10) = hmf1
p(ii,jj,11) = fof2
p(ii,jj,12) = ymf2
p(ii,jj,13) = hmf2
c
c calculate the QP profile parameters at each location
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
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
c
c calculate profile at TX, & write to file
if(ii.eq.1.and.jj.eq.1) go to 1998
go to 2001
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.) write(4,400) i,ht(i),fn(i)
if(fn(i).gt.0.) write(6,400) i,ht(i),fn(i)
1999 continue
400 format(i4,f6.1,f6.2)
c
2001 continue
c
c write out the parameters of the model
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
go to 10
c
end