program rd_eis6800 C 04/11: rd_eis6800.f for EISCAT ISR tau2 6800 kindat SW 77.5 elev 1 min 76-732km (166 km) C skipping kindat 6801 power profile C Program to read the ascii format data submitted to the CEDAR DB C .asc not available from DB access system directly, but can be converted from .cbf character*80 char80 dimension lprol(16),icod1d(100),ipar1d(100),icod2d(100), | ipar2d(100,900),itime(4), fillne(43),te(43),erte(43), | htjm(43),corrne(43),ercorrne(43),rnelog10(43),h(3),p(3), | a20ne(42,21),a20erne(42,21),a20te(42,21),a20erte(42,21), | a20ti(42,21),a20erti(42,21),avne(42),stdne(42),avte(42), | avti(42),a20ut(21),a20day(21),a20moda(21) missing = -32767 ird = 11 C Open 2 files - one for 1D TEC, hmax, Nmax, and FoF2, and the other for 2D ht profiles iw1 = 12 iw2 = 13 C Make 2 extra files of medians every 20 min i20w1 = 14 i20w2 = 15 nm = 0 glat = 69.58 glon = 19.23 ida = 14 C ida = 15 if (ida .eq. 14) then open (ird,file='NCAR_2006-12-14_tau2pl_uhf.bin',status='old') C Set previous minute (mp) and utp to same as first time read mp = 1 utp = 0. open (iw1,file='eis061214_nmax') open (iw2,file='eis061214_profiles') open (i20w1,file='eis061214_20m_nmax') open (i20w2,file='eis061214_20m_profiles') else open (ird,file='NCAR_2006-12-15_tau2pl_uhf.bin',status='old') C Set previous minute (mp) and utp to same as first time read mp = 1 utp = 11. open (iw1,file='eis061215_nmax') open (iw2,file='eis061215_profiles') open (i20w1,file='eis061215_20m_nmax') open (i20w2,file='eis061215_20m_profiles') endif write (iw1,"(' year dayn moda uth glat glon', | ' hmax foF2MHz Nmax(m-3) TEC(m-2)')") write (iw2,"(' year dayn uth alt(km) Ti(K) ', | 'errTi Te(K) errTe lg10Ne(m-3) errNe')") write (i20w1,"(' year dayn moda uth glat glon', | ' hmax foF2MHz Nmax(m-3) TEC(m-2)')") write (i20w2,"(' year dayn uth alt(km) #Ti Ti(K) ', | 'stdTi #Te Te(K) stdTe #Ne lg10Ne(m-3) stdNe')") C Ne,Ti,Te/Ti, errs at 6,7,8, 11,12,13 jne = 6 jerne = 11 jti = 7 jerti = 12 jteti = 8 jerteti = 13 rewind ird nrecrd = 0 1000 continue c write (6,"(1x,'Read nrecrd =',i6)") nrecrd C Read Prologue nrecrd = nrecrd + 1 read (ird,"(16i6)",end=2000) lprol C Read (and write) text portions (headers 3101 or catalogs 2101) if (lprol(2) .eq. 3101 .or. lprol(2) .eq. 2101) then if (nrecrd.eq.1) write (6,"(1x,16i6)") lprol do 1100 n=2,lprol(1) read (ird,"(a80)") char80 c if (nrecrd.eq.1) write (6,"(1x,a80)") char80 1100 continue go to 1000 endif kinst = lprol(3) kindat = lprol(4) ibyrt = lprol(5) ibdtt = lprol(6) ibhmt = lprol(7) ibcst = lprol(8) ieyrt = lprol(9) iedtt = lprol(10) iehmt = lprol(11) iecst = lprol(12) jpar = lprol(14) mpar = lprol(15) nrows = lprol(16) C Convert year mmdd hhmm csec to daynumber ndayn itime(1) = ibyrt itime(2) = ibdtt itime(3) = ibhmt itime(4) = ibcst call cvt2dn (6,itime,ndayn) C Read JPAR 1-D parameters, and NROWS of MPAR 2-D parameters from instrument KINST C between begin date IBYRT,IBDTT,IBHMT,IBCST and end date IEYRT,IEDTT,IEHMT,IECST C (year, monthday, UThourmin, centisecond) c write (6,"(1x,'ndayn kinst kindat begin,end_time 1-D 2-D rows=', c | 14i6)") ndayn, (lprol(i1),i1=3,12),(lprol(i2),i2=14,16) if (jpar .gt. 100 .or. mpar .gt. 100 .or. nrows .gt. 900) stop C Are there any JPAR (single valued parameters)? if (jpar .gt. 0) then C Read 1-D parameters C Read 1-D parameter codes read (ird,"(20i6)") (icod1d(i3),i3=1,jpar) c write (6,"(1x,'1D codes =',100i6)") (icod1d(n),n=1,jpar) C Read 1-D parameter values read (ird,"(20i6)") (ipar1d(i4),i4=1,jpar) endif C Are there any MPAR (multiple valued parameters)? if (mpar .gt. 0) then C Read 2-D parameters C Read 2-D parameter codes read (ird,"(20i6)") (icod2d(i5),i5=1,mpar) c write (6,"(1x,'2D codes =',100i6)") (icod2d(n),n=1,mpar) C Read 2-D parameter values do 1600 j=1,nrows read (ird,"(20i6)") (ipar2d(i6,j),i6=1,mpar) 1600 continue C Find if kindat=6800 if (kindat .eq. 6800) then ihr = ibhmt/100 imn = ibhmt - ihr*100 uth = ihr + float(imn)/60. C Set altitude and Ne and fill in gaps in Ne (jne) above 100 km (j100) j100 = 9 C Look at peak >100km instead of >150km j150 = j100 do 1602 j=1,nrows C Set htjm altitude htjm(j) = ipar2d(1,j)+ipar2d(2,j)*0.0001 C Delete if Ne error bar > Ne if (ipar2d(jerne,j) .gt. ipar2d(jne,j)) then ipar2d(jne,j) = missing ipar2d(jerne,j) = missing endif C Set Ne in m-3 rnelog10(j) = ipar2d(jne,j) corrne(j) = missing ercorrne(j) = missing if (ipar2d(jne,j) .gt. 0.) then corrne(j) = 10.**(ipar2d(jne,j)*0.001) ercorrne(j) = 10.**(ipar2d(jerne,j)*0.001) endif 1602 continue C write (6,"(1x,' ipar2d(jne,j) =',20i6)") (ipar2d(jne,j),j=1,nrows) C Set top rnelog10 if missing if (rnelog10(nrows) .lt. 0.) then do 1610 j=nrows-1,j100,-1 if (rnelog10(j) .gt. 0. .and. rnelog10(j-1) .gt. 0. | .and. rnelog10(j-1) .gt. rnelog10(j)) then diffne = (rnelog10(j-1)-rnelog10(j))/(htjm(j)-htjm(j-1)) C write (6,"(1x,'diffne rnelog10(j-1,j) htjm(j-1,j,nrows) =', C | 3e12.4,3f7.2)") diffne,rnelog10(j-1),rnelog10(j),htjm(j-1), C | htjm(j),htjm(nrows) do 1605 j1=j+1,nrows if (rnelog10(j1) .lt. 0.) then rnelog10(j1) = rnelog10(j) - diffne*(htjm(j1)-htjm(j)) C write (6,"(1x,'rnelog10 htjm =',e12.4,f7.2)") rnelog10(j1), C | htjm(j1) endif 1605 continue endif 1610 continue endif C write (6,"(1x,' rnelog10(nrows) =',e12.4)") rnelog10(nrows) C Set "bottom" rnelog10 if missing if (rnelog10(j100) .lt. 0.) then do 1615 j=j100+1,nrows if (rnelog10(j) .gt. 0.) go to 1616 1615 continue 1616 j1 = j do 1617 j=j1,nrows if (rnelog10(j) .gt. 0.) go to 1618 1617 continue 1618 j2 = j diffne = (rnelog10(j2)-rnelog10(j1))/(htjm(j2)-htjm(j1)) rnelog10(j100) = rnelog10(j1) - diffne*(htjm(j1)-htjm(j100)) endif C write (6,"(1x,' rnelog10(11) =',e12.4)") rnelog10(11) C Fill in missing rnelog10 values as linear interpolations do 1625 j=j100+1,nrows-1 if (rnelog10(j) .lt. 0.) then do 1619 j1=j,nrows-1 if (rnelog10(j1) .gt. 0.) go to 1620 1619 continue 1620 j2 = j1 diffne = (rnelog10(j1)-rnelog10(j-1))/(htjm(j1)-htjm(j-1)) rnelog10(j) = rnelog10(j-1) + diffne*(htjm(j)-htjm(j-1)) endif 1625 continue C Calculate hmF2, NmF2, and TEC (1 TECU=10^16elec/m2) xnemax = 0. jx = 0 do 1626 j=j150,nrows-1 C Find max Ne from 100-674 km if (corrne(j) .gt. xnemax) then xnemax = corrne(j) jx = j endif 1626 continue if (jx .eq. 0) then hmf2 = missing xnmax = missing fof2 = missing go to 1632 endif h(2) = htjm(jx) p(2) = corrne(jx) hmf2 = htjm(jx) xnmax = corrne(jx) C TEMP if (abs(uth-6.45) .lt. 0.01 .or. abs(uth-6.72) .lt. 0.01 | .or. abs(uth-14.20) .lt. 0.01) then write (6,"(1x,' uth nmax hmax =',f6.2,e12.4,f10.2)") | uth,xnmax,hmf2 endif C Check to see if max Ne next to missing data or not do 1627 j=jx-1,j150,-1 if (corrne(j) .gt. 0.) go to 1628 1627 continue go to 1631 1628 h(1) = htjm(j) p(1) = corrne(j) do 1629 j=jx+1,nrows if (corrne(j) .gt. 0.) go to 1630 1629 continue go to 1631 1630 h(3) = htjm(j) p(3) = corrne(j) if (j .eq. nrows .and. p(3) .gt. p(2)) go to 1631 C Use polynomial to find hmF2 and NmF2 call hnmf2 (h,p,hmf2,xnmax) C TEMP if (abs(uth-6.45) .lt. 0.01 .or. abs(uth-6.72) .lt. 0.01 | .or. abs(uth-14.20) .lt. 0.01) then write (6,"(1x,' nmax p hmax h =',4e12.4,4f7.2)") | xnmax,p,hmf2,h endif 1631 fof2 = sqrt(xnmax/1.24e+10) 1632 continue C Find tec from non-missing data tec = 0. htjm1 = htjm(j100-1) do 1635 j=j100,nrows if (corrne(j) .gt. 0.) then tec = tec + corrne(j)*(htjm(j)-htjm1)*1.e+3 htjm1 = htjm(j) endif 1635 continue C Find tec from filled in data c tecn = 0. c do 1640 j=j100,nrows c fillne(j) = 10.**(rnelog10(j)*0.001) c tecn = tecn + fillne(j)*(htjm(j)-htjm(j-1))*1.e+3 c1640 continue c write (6,"(1x,' tec,n nmax hmax =',3e12.4,f7.2)") c | tec,tecn,xnmax,hmf2 C Use polynomial to find hmF2 and NmF2 C call hnmf2 (htjm,corrne,hmf2,xnmax,jx) C fof2 = sqrt(xnmax/1.24e+10) write (iw1,"(1x,i4,2i5,f6.2,3f7.2,f6.2,2e12.4)") ibyrt, | ndayn,ibdtt,uth,glat,glon,hmf2,fof2,xnmax,tec C Save only for 100 km (j100) and up do 1675 j=j100,nrows te(j) = missing erte(j) = missing C Delete if Ti error bar>200K, or if Te/Ti<0.2 or >5.0 if (ipar2d(jerti,j) .gt. 200. .or. ipar2d(jteti,j) .lt. 200 | .or. ipar2d(jteti,j) .gt. 5000) then ipar2d(jerti,j) = missing ipar2d(jti,j) = missing te(j) = missing erte(j) = missing endif C Delete if Ti<500K when ht>300km if (ipar2d(jti,j) .lt. 500. .and. htjm(j) .gt. 300.) then ipar2d(jerti,j) = missing ipar2d(jti,j) = missing te(j) = missing erte(j) = missing endif if (ipar2d(jti,j) .gt. 0) then te(j) = ipar2d(jteti,j)*0.001*ipar2d(jti,j) if (ipar2d(jerteti,j) .gt. 0) then erte(j) = ipar2d(jerteti,j)*0.001*ipar2d(jti,j) C Delete if Te error bar > 200K if (erte(j) .gt. 200.) then ipar2d(jerti,j) = missing ipar2d(jti,j) = missing te(j) = missing erte(j) = missing endif endif endif write (iw2,"(1x,i4,i5,f6.2,f7.2,2i7,2f8.0,2f8.3)") ibyrt, | ndayn,uth,htjm(j),ipar2d(jti,j),ipar2d(jerti,j),te(j), | erte(j),ipar2d(jne,j)*0.001,ipar2d(jerne,j)*0.001 1675 continue C Gather data for medians every 20 min at 00,20,40 min if (imn .ge. 50 .or. imn .lt. 10) mm = 1 if (imn .ge. 10 .and. imn .lt. 30) mm = 2 if (imn .ge. 30 .and. imn .lt. 50) mm = 3 if (mm .eq. mp .and. abs(uth-utp) .lt. 0.4) then nm = nm + 1 else C Calculate aves from nm values at mp, then reset nm nav = 0 avut = 0. do 1699 n=1,nm avut = avut + a20ut(n) 1699 continue avut = avut/nm mdayn = a20day(nm) moda = a20moda(nm) if (a20day(1) .lt. a20day(nm)) then avut = a20ut(nm) do 1700 n=1,nm-1 if (a20day(n) .lt. a20day(nm)) a20ut(n) = a20ut(n)-24. avut = avut + a20ut(n) 1700 continue avut = avut/nm if (avut .lt. 0.) then avut = avut + 24. mdayn = a20day(1) moda = a20moda(1) endif endif do 1710 j=1,nrows nav = 0 avne(j) = 0. do 1701 n=1,nm if (a20ne(j,n) .gt. 0.) then nav = nav + 1 avne(j) = avne(j) + a20ne(j,n) endif 1701 continue if (nav .gt. 0) then avne(j) = avne(j)/nav else avne(j) = missing endif nne = nav nav = 0 avte(j) = 0. do 1702 n=1,nm if (a20te(j,n) .gt. 0.) then nav = nav + 1 avte(j) = avte(j) + a20te(j,n) endif 1702 continue if (nav .gt. 0) then avte(j) = avte(j)/nav else avte(j) = missing endif nte = nav nav = 0 avti(j) = 0. do 1703 n=1,nm if (a20ti(j,n) .gt. 0.) then nav = nav + 1 avti(j) = avti(j) + a20ti(j,n) endif 1703 continue if (nav .gt. 0) then avti(j) = avti(j)/nav else avti(j) = missing endif nti = nav C Find standard deviations for the error bar stdne(j) = 0. if (nne .gt. 0) then do 1706 n=1,nm if (a20ne(j,n) .gt. 0.) | stdne(j) = stdne(j) + (avne(j)-a20ne(j,n))**2 1706 continue stdne(j) = sqrt(stdne(j)/nne) endif stdte = 0. if (nte .gt. 0) then do 1707 n=1,nm if (a20te(j,n) .gt. 0.) | stdte = stdte + (avte(j)-a20te(j,n))**2 1707 continue stdte = sqrt(stdte/nte) endif stdti = 0. if (nti .gt. 0) then do 1708 n=1,nm if (a20ti(j,n) .gt. 0.) | stdti = stdti + (avti(j)-a20ti(j,n))**2 1708 continue stdti = sqrt(stdti/nti) endif if (j .gt. 10) | write (i20w2,"(1x,i4,i5,f6.2,f7.2,i4,2f9.1,i4,2f9.1,i4,2e12.4)") | ibyrt,mdayn,avut,htjm(j),nti,avti(j),stdti,nte,avte(j), | stdte,nne,avne(j),stdne(j) 1710 continue C Calculate hmF2, NmF2, and TEC (1 TECU=10^16elec/m2) xnemax = 0. do 1726 j=j150,nrows-1 C Find max Ne from 154.6-536.8 km for periods with more than 1 values in bin if (avne(j) .gt. xnemax .and. stdne(j) .gt. 0.01) then xnemax = avne(j) jx = j endif 1726 continue h(2) = htjm(jx) p(2) = avne(jx) hmf2 = htjm(jx) xnmax = avne(jx) C Check to see if max Ne next to missing data or not do 1727 j=jx-1,j150,-1 if (avne(j) .gt. 0. .and. stdne(j) .gt. 0.01) go to 1728 1727 continue go to 1731 1728 h(1) = htjm(j) p(1) = avne(j) do 1729 j=jx+1,nrows if (avne(j) .gt. 0. .and. stdne(j) .gt. 0.01) go to 1730 1729 continue go to 1731 1730 h(3) = htjm(j) p(3) = avne(j) if (j .eq. nrows .and. p(3) .gt. p(2)) go to 1731 C Use polynomial to find hmF2 and NmF2 call hnmf2 (h,p,hmf2,xnmax) 1731 fof2 = sqrt(xnmax/1.24e+10) C Find tec from non-missing data tec = 0. htjm1 = htjm(j100-1) do 1735 j=j100,nrows if (avne(j) .gt. 0. .and. stdne(j) .gt. 0.01) then tec = tec + avne(j)*(htjm(j)-htjm1)*1.e+3 htjm1 = htjm(j) endif 1735 continue write (i20w1,"(1x,i4,2i5,f6.2,3f7.2,f6.2,2e12.4)") ibyrt, | mdayn,moda,avut,glat,glon,hmf2,fof2,xnmax,tec nm = 1 endif mp = mm utp = uth a20day(nm) = ndayn a20moda(nm) = ibdtt a20ut(nm) = uth do 1750 j=1,nrows a20ne(j,nm) = corrne(j) a20erne(j,nm) = ercorrne(j) a20te(j,nm) = te(j) a20erte(j,nm) = erte(j) a20ti(j,nm) = ipar2d(jti,j) a20erti(j,nm) = ipar2d(jerti,j) 1750 continue endif endif go to 1000 2000 continue C write (6,"(1x,'Stop after reading ',i6,' records')") nrecrd stop end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE CVT2DN(MSGUN,ITIME,NDA) C This converts the standard I.S. format time ITIME(4) C yyyy, mmdd, hhmm & csec to day number of year NDA. DIMENSION ITIME(4),NUMD(12) INTEGER DD SAVE NUMD DATA NUMD/0,31,28,31,30,31,30,31,31,30,31,30/ MMDD=ITIME(2) MM=MMDD/100 DD=MMDD - MM*100 C Leap year adjustment NUMD(3)=28 IF(MOD(ITIME(1),4) .EQ. 0)NUMD(3)=29 NDA = 0 DO 100 I=1,MM 100 NDA = NDA+NUMD(I) NDA = NDA+DD IF(NDA .LT. 1 .OR. NDA.GT.366) + WRITE(MSGUN,'('' CVT2DN: Bogus day number'',I9,'', computed'', + '' from '',4I5)') NDA,(ITIME(I),I=1,4) RETURN END cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C subroutine hnmf2(ht,rne,hmf2,nmf2,kx) subroutine hnmf2(h,p,hmf2,nmf2) C Find the NmF2 and hmF2 based on the electron profile. C btf: This is taken from an idl procedure provided by Liying. real ht(43),rne(43),nmf2,hmf2 real h(3),p(3),h12,h22,h32,deltx,atx,ax,btx,bx,ctx,cx c h(1) = ht(kx-1) c h(2) = ht(kx) c h(3) = ht(kx+1) c p(1) = rne(kx-1) c p(2) = rne(kx) c p(3) = rne(kx+1) h12 = h(1)*h(1) h22 = h(2)*h(2) h32 = h(3)*h(3) deltx=h12*h(2)+h22*h(3)+h32*h(1)-h32*h(2)-h12*h(3)-h22*h(1) atx=p(1)*h(2)+p(2)*h(3)+p(3)*h(1)-h(2)*p(3)-h(3)*p(1)-h(1)*p(2) ax=atx/deltx btx=h12*p(2)+h22*p(3)+h32*p(1)-h32*p(2)-h12*p(3)-h22*p(1) bx=btx/deltx ctx=h12*h(2)*p(3)+h22*h(3)*p(1)+h32*h(1)*p(2)-h32*h(2)*p(1)- | h12*h(3)*p(2)-h22*h(1)*p(3) cx=ctx/deltx hmf2=-(bx/(2.*ax)) nmf2=-((bx*bx-4.*ax*cx)/(4.*ax)) c write(6,"(1x,'hnmf2 nmf2 =',f7.2,e12.4)") hmf2,nmf2 end