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 2918 if (hr(i).gt.hend) hr(i)=hend 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