##// END OF EJS Templates
DATABASES » madrigal
RSS
Update configs
Update configs
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r9:2e29893b10f5
Show More
iritec.f
278 lines | 8.3 KiB | text/x-fortran | FortranFixedLexer
/ source / madf / geolib / iritec.f
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 C =====================================================================
C Downloaded from ftp://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/
C iri2001/fortran_code/
C by B. Rideout on May 4, 2007. Dated 2/25/2005 on ftp site.
C A description of this code can be found
C at http://modelweb.gsfc.nasa.gov/ionos/iri.htm
C
C $Id: iritec.f 7348 2021-03-31 15:03:00Z brideout $
C
c iritec.for, version number can be found at the end of this comment.
c-----------------------------------------------------------------------
C
C contains IRIT13, IONCORR, IRI_TEC subroutines to computed the
C total ionospheric electron content (TEC)
C
c-----------------------------------------------------------------------
C Corrections
C
C 3/25/96 jmag in IRIT13 as input
C 8/31/97 hu=hr(i+1) i=6 out of bounds condition corrected
C 9/16/98 JF(17) added to input parameters; OUTF(11,50->100)
C ?/ ?/99 Ne(h) restricted to values smaller than NmF2 for topside
C 11/15/99 JF(20) instead of JF(17)
C 10/16/00 if hr(i) gt hend then hr(i)=hend
C 12/14/00 jf(30),outf(20,100),oarr(50)
C
C Version-mm/dd/yy-Description (person reporting correction)
C 2000.01 05/07/01 current version
c 2000.02 10/28/02 replace TAB/6 blanks, enforce 72/line (D. Simpson)
c 2000.03 11/08/02 common block1 in iri_tec with F1reg
c-----------------------------------------------------------------------
C
C
subroutine IRIT13(ALATI,ALONG,jmag,jf,iy,md,hour,hbeg,hend,
& tec,tecb,tect)
c-----------------------------------------------------------------------
c Program for numerical integration of IRI-94 profiles from h=100km
C to h=alth.
C
C INPUT: ALATI,ALONG LATITUDE NORTH AND LONGITUDE EAST IN DEGREES
C jmag =0 geographic =1 geomagnetic coordinates
C jf(1:30) =.true./.false. flags; explained in IRISUB.FOR
C iy,md date as yyyy and mmdd (or -ddd)
C hour decimal hours LT (or UT+25)
c hbeg,hend upper and lower integration limits in km
C
C OUTPUT: TEC Total Electron Content in M-2
C tecb,tect percentage of bottomside and topside content
c-----------------------------------------------------------------------
dimension outf(20,100),oarr(50)
logical jf(30)
c
c select various options and choices for IRI-94
c
tec = -111.
tect= -111.
tecb= -111.
c
c initialize IRI parameter in COMMON blocks
c
abeg=hbeg
aend=hend
astp=hend-hbeg
call IRI_SUB(JF,JMAG,ALATI,ALONG,IY,MD,HOUR,
& abeg,aend,astp,OUTF,OARR)
c
c calculate total electron content (TEC) in m-2 using the
c stepsize selection 2 (highest accuracy)
c
call iri_tec (hbeg,hend,2,tec,tect,tecb)
return
end
c
c
real function ioncorr(tec,f)
c-----------------------------------------------------------------------
c computes ionospheric correction IONCORR (in m) for given vertical
c ionospheric electron content TEC (in m-2) and frequency f (in Hz)
c-----------------------------------------------------------------------
ioncorr = 40.3 * tec / (f*f)
return
end
c
c
subroutine iri_tec (hstart,hend,istep,tectot,tectop,tecbot)
c-----------------------------------------------------------------------
C subroutine to compute the total ionospheric content
C INPUT:
C hstart altitude (in km) where integration should start
C hend altitude (in km) where integration should end
C istep =0 [fast, but higher uncertainty <5%]
C =1 [standard, recommended]
C =2 [stepsize of 1 km; best TEC, longest CPU time]
C OUTPUT:
C tectot total ionospheric content in tec-units (10^16 m^-2)
C tectop topside content (in %)
C tecbot bottomside content (in %)
C
C The different stepsizes for the numerical integration are
c defined as follows (h1=100km, h2=hmF2-10km, h3=hmF2+10km,
c h4=hmF2+150km, h5=hmF2+250km):
C istep h1-h2 h2-h3 h3-h4 h4-h5 h5-hend
C 0 2.0km 1.0km 2.5km exponential approximation
C 1 2.0km 1.0km 2.5km 10.0km 30.0km
C 2 1.0km 0.5km 1.0km 1.0km 1.0km
C
c-----------------------------------------------------------------------
logical expo
dimension step(5),hr(6)
common /block1/hmf2,xnmf2,hmf1,f1reg
logical f1reg
ctest
save
expo = .false.
numstep = 5
xnorm = xnmf2/1000.
hr(1) = 100.
hr(2) = hmf2-10.
hr(3) = hmf2+10.
hr(4) = hmf2+150.
hr(5) = hmf2+250.
hr(6) = hend
do 2918 i=2,6
if (hr(i).gt.hend) hr(i)=hend
2918 continue
if (istep.eq.0) then
step(1)=2.0
step(2)=1.0
step(3)=2.5
step(4)=5.
if (hend.gt.hr(5)) expo=.true.
endif
if (istep.eq.1) then
step(1)=2.0
step(2)=1.0
step(3)=2.5
step(4)=10.0
step(5)=30.0
endif
if (istep.eq.2) then
step(1)=1.0
step(2)=0.5
step(3)=1.0
step(4)=1.0
step(5)=1.0
endif
sumtop = 0.0
sumbot = 0.0
C
C find the starting point for the integration
C
i=0
ia=1
3 i=i+1
h=hr(i)
if(hstart.gt.h) then
hr(i)=hstart
ia=i
goto 3
endif
C
C start the numerical integration
C
i=ia
h=hr(i)
hu=hr(i+1)
delx = step(i)
1 h = h + delx
hh = h
if (h.ge.hu) then
delx = hu - h + delx
hx = hu - delx/2.
YNE = XE_1(hx)
if((hx.gt.hmf2).and.(yne.gt.xnmf2)) yne=xnmf2
yyy = yne * delx / xnorm
i=i+1
if(i.lt.6) then
h = hr(i)
hu = hr(i+1)
delx = step(i)
endif
else
hx = h - delx/2.
YNE = XE_1(hx)
if((hx.gt.hmf2).and.(yne.gt.xnmf2)) yne=xnmf2
yyy = yne * delx / xnorm
endif
if (hx.le.hmf2) then
sumbot = sumbot + yyy
else
sumtop = sumtop + yyy
endif
if (expo.and.(hh.ge.hr(4))) goto 5
if (hh.lt.hend.and.i.lt.6) goto 1
zzz = sumtop + sumbot
tectop = sumtop / zzz * 100.
tecbot = sumbot / zzz * 100.
tectot = zzz * xnmf2
return
5 num_step = 3
hei_top = hr(4)
hei_end = hend
top_end = hei_end - hei_top
del_hei = top_end / num_step
xntop = xe_1(hei_end)/xnmf2
if(xntop.gt.0.9999) then
ss_t = top_end
goto 2345
endif
hei_2 = hei_top
hei_3 = hei_2 + del_hei
hei_4 = hei_3 + del_hei
hei_5 = hei_end
hss = top_end / 4.
C hss = 360.
xkk = exp ( - top_end / hss ) - 1.
x_2 = hei_2
x_3 =hei_top-hss*alog(xkk*(hei_3 - hei_top)/top_end + 1.)
x_4 =hei_top-hss*alog(xkk*(hei_4 - hei_top)/top_end + 1.)
x_5 = hei_end
ed_2 = xe_1(x_2)/xnmf2
if(ed_2.gt.1.) ed_2=1.
ed_3 = xe_1(x_3)/xnmf2
if(ed_3.gt.1.) ed_3=1.
ed_4 = xe_1(x_4)/xnmf2
if(ed_4.gt.1.) ed_4=1.
ed_5 = xntop
if(ed_3.eq.ed_2) then
ss_2 = ed_3 * (x_3 - x_2)
else
ss_2=( ed_3 - ed_2 ) * ( x_3 - x_2 ) / alog ( ed_3 / ed_2 )
endif
if(ed_4.eq.ed_3) then
ss_3 = ed_4 * (x_4 - x_3)
else
ss_3=( ed_4 - ed_3 ) * ( x_4 - x_3 ) / alog ( ed_4 / ed_3 )
endif
if(ed_5.eq.ed_4) then
ss_4 = ed_5 * (x_5 - x_4)
else
ss_4=( ed_5 - ed_4 ) * ( x_5 - x_4 ) / alog ( ed_5 / ed_4 )
endif
ss_t = ss_2 + ss_3 + ss_4
2345 sumtop = sumtop + ss_t * 1000.
zzz = sumtop + sumbot
tectop = sumtop / zzz * 100.
tecbot = sumbot / zzz * 100.
tectot = zzz * xnmf2
RETURN
END