program ionmod_dud
c
c prog to determine Dudeney 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 29-may-87
c modified 15-jul-91 to accept full unformatted maps of fof2
c--- last modified 15-jul-91
c
c LINK with files ips,layers,MSIS83,ilios,spoints,dudsub,
c dud_jones,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 dudsub dimsol,[hmax,ymf2 - not used]
c dud_jones dudney [only for looking at the N(h) profile]
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
c--- output data .. in the file IONMOD_DUD.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,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 (no need for ymf1 & hmf1)
c k=9 foF2 F2 layer
c k=10 ymf2
c k=11 hmf2
c k=12 h0 base of ionosphere -- DUDENEY parameters
c k=13 ad
c k=14 tau
c k=15 hd
c k=16 beta
c k=17 ht
c k=18 ft
c k=19 at
c
dimension p(50,5,20),ap(7),flatt(50),flonn(50),nday(12)
dimension fn(3,5,275)
character*10 munth(12)
character*8 modx
LOGICAL LCOND
REAL M3000
real*8 X,freq,hh
COMMON/IONDAT/FOF2, FOF1, FOE, HFF2, M3000, HME, YME,
* HMF2, YMF2, H0, XD, DV
COMMON/PARMS/BETA, HT, FT, At, HF, AD, TAU, HD, LCOND
common / XX / modx(2),X, plus, other, stuff
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 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_DUD.DAT
open (unit=1,file='ionmod_dud.dat',status='old')
c grid-point data file on logical unit 2 - see later for OPEN
cc open (unit=2,file=munth(mon),form='unformatted',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,20
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 the virtual height of the base of the F1/F2 layer is not used
hff2 = 0.
c
write(6,600)
600 format(//5x,'Prog to set up the profile parameters for',/
a ' a 3-D ionospheric model, using the DUDENEY 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),form='unformatted',
a 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
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 use default flux values if not supplied (set to 0.)
if(f107.eq.0.) f107 = flux12
if(f107a.eq.0.) f107a = flux12
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
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 layer and Dudeney 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,
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,mon,uthour,hourlt,idayno,flat,flon,
a foe,tindex,f107a,f107,ap,fof2,ymf2,hmf2)
c
c**********************************
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
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) = fof2
p(ii,jj,10) = ymf2
p(ii,jj,11) = hmf2
c
c calculate the Dudeney profile parameters at each location
c base of layer
h0 = 80.
c xd -- D and E layers touch at xd * foe
xd = 0.5
call DIMSOL (ierror)
if(ierror.eq.0) go to 1200
write(6,1150) ierror,flat,flon
1150 format(///5x,' *** error in DIMSOL ',i4, 'lat&lon =',2f6.2)
stop
1200 continue
c
c store the Dudeney parameters
p(ii,jj,12) = h0
p(ii,jj,13) = ad
p(ii,jj,14) = tau
p(ii,jj,15) = hd
p(ii,jj,16) = beta
p(ii,jj,17) = ht
p(ii,jj,18) = ft
p(ii,jj,19) = at
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
c we use s/r DUDNEY in file DUD_JONES, which was not written for the
c present purposes
c set frequency = 1., to get fn, not X
freq = 1.d0
ihs = 80
ihf = 350
c we don't want to set up the profile parameters again
ipass = 1
do 180 ih = ihs,ihf
iih = ih - ihs + 1
h = ih
hh = dble (h)
call dudney (ipass,freq,hh)
c result is returned in labelled common / xx/
fni = sngl(x)
fni = sqrt(fni)
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,20)
334 format(5e16.8)
c
type*,' Length of circuit & bearing ',distsv,bearsv
type*,' Special bearing case ',ispec
c
go to 10
end