program ionmod_IRI
c
c prog to determine IRI profile parameters at a grid of
c locations, for use in 3-D raytracing using the Jones program
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 10-jul-87
c--- last modified 11-aug-87
c
c LINK with files ips,fof2m3_ips,ilios,spoints,irinhpro,
c irinhsubs,xesubs,ggm,solar
c ips loadgp,f2m3,tcon,gyro
c fof2m3_ips fof2m3_ips
c ilios
c spoints spoints,distbear,latlon
c irinhpro irinhpro, iondem (currently a DUMMY file)
c irinhsubs all the IRI s/rs required for n(h)
c ggm
c solar
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 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
c--- output data .. in the file IONMOD_IRI.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,30 where
c k=1 colatitude of point (degrees)
c k=2 geographic latitude (degrees)
c k=3 geomagnetic colatitude (degrees)
c k=4 geomagnetic longitude
c k=5 hmf2 --- the IRI parameters
c k=6 nmf2
c k=7 hmf1
c k=8 b0
c k=9 b1
c k=10 c1
c k=11 hz
c k=12 t
c k=13 hst
c k=14 hdx
c k=15 hme
c k=16 nme
c k=17 hmd
c k=18 nmd
c k=19 hef
c k=20 d1
c k=21 xkk
c k=22 fp30
c k=23 fpu
c k=24 fp1
c k=25 fp2
c k=26 night
c k=27 e(1)
c k=28 e(2)
c k=29 e(3)
c k=30 e(4)
c
dimension p(50,5,30),flatt(50),flonn(50),nday(12),e(4)
dimension fn(3,5,275)
character*10 munth(12)
character*8 modx
LOGICAL LCOND,night
REAL M3000,nmd,nme,nmf2
c IRI labelled common
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR
& /BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,XKK,FP30,FP3U,FP1,FP2
& /BLOCK5/ZX,TNX,ATN,CTN/BLOCK6/NIGHT,E/BLOTE/AHH(7),ATE1,
& STTE(6),DTE(5)/CONST/UMR
& /BLOCK8/HS,TNHS,XSM,MM,MXSM/BLO10/BETA,ETA,DELTA,ZETA
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 = 8.0
c lat and lon spacing of grid for special cases
stplon = 2.0
stplat = 2.0
c
dv = 0.1
c
c size of array dimension
ndim = 50
c set the year in concrete
lyear = 87
c
c output is to file IONMOD_IRI.DAT
open (unit=1,file='ionmod_iri.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,30
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 IRI86 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,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
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
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
type*,' flat1,flon1,flat2,flon2 ', flat1,flon1,flat2,flon2
call spoints (flat1,flon1,flat2,flon2,nip,flatt,flonn)
type*,' flat1,flon1,flat2,flon2 ', flat1,flon1,flat2,flon2
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 IRI 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 F2-layer parameters - fof2 and m(3000)f2
call fof2m3_ips
a (lunmap,mon,uthour,flat,flon,tindex,fof2,m3000)
c
c********************************** can't do this for IRI
c keep a pretend F1 layer to avoid a very wide E-F step when the
c calculated value of foF1 is zero
ccc if(fof1.eq.0.) fof1 = 1.4 * foe
c make sure that fof2 exceeds fof1
ccc if(fof1.ge.fof2) fof1 = (foe + fof2) / 2.
c**********************************
c
c calculate the IRI profile parameters at each location
axhi = -100.0
c hmaxf2 is a supplied value (optional)
hmaxf2 = 0.
call irinhpro
a (flat,flon,0,r12,lmon,lday,hourlt,fof2,m3000,axhi,hmaxf2)
c
c store the IRI parameters
p(ii,jj,5) = hmf2
p(ii,jj,6) = nmf2
p(ii,jj,7) = hmf1
p(ii,jj,8) = b0
p(ii,jj,9) = b1
p(ii,jj,10) = c1
p(ii,jj,11) = hz
p(ii,jj,12) = t
p(ii,jj,13) = hst
p(ii,jj,14) = hdx
p(ii,jj,15) = hme
p(ii,jj,16) = nme
p(ii,jj,17) = hmd
p(ii,jj,18) = nmd
p(ii,jj,19) = hef
p(ii,jj,20) = d1
p(ii,jj,21) = xkk
p(ii,jj,22) = fp30
p(ii,jj,23) = fp3u
p(ii,jj,24) = fp1
p(ii,jj,25) = fp2
p(ii,jj,26) = 0.
if (night) p(ii,jj,26) = 1.
p(ii,jj,27) = e(1)
p(ii,jj,28) = e(2)
p(ii,jj,29) = e(3)
p(ii,jj,30) = e(4)
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
if(nop.lt.npro) npro = nop
if(ii.gt.npro) go to 250
ihs = 80
ihf = 350
do 180 ih = ihs,ihf
iih = ih - ihs + 1
h = ih
ed = xe(h)
fni = sqrt(ed/1.24e10)
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
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,30)
334 format(5e16.8)
c
type*,' Length of circuit & bearing ',distsv,bearsv
type*,' Special bearing case ',ispec
c
go to 10
end