program rd_son5521 C 04/11: rd_son5521.f for Sondrestrom ISR ~4.5min 3-position (Z, NE+NW 75elev) 70-563/578km C Pick every Z position (cannot average w 75elev since at diff alts) from 100-578km 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), | htjm(43),corrne(43),ercorrne(43),rnelog10(43),h(3),p(3), | a20ne(42,5),a20erne(42,5),a20te(42,5),a20erte(42,5), | a20ti(42,5),a20erti(42,5),avne(42),averne(42),avte(42), | averte(42),avti(42),averti(42),a20ut(5),a20day(5),a20moda(5) character*2 dir 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 C son050514g.101 open (ird,file='son050514g.001',status='old') C Set previous minute to same as first time read mp = 3 nm = 0 open (iw1,file='son050514_15_nmax') write (iw1,"(' year dayn moda uth glat glon', | ' hmax foF2MHz Nmax(m-3) TEC(m-2)')") open (iw2,file='son050514_15_profiles') glat = 66.99 glon = -50.95 c write (iw2,"('Sondrestrom kindat=5521 ~4.5 min 3 position scan' c | /' Use Zenith position only 32 hts, glat glon =',2f8.2)") c | glat,glon write (iw2,"(' year dayn uth alt(km) Ti(K) ', | 'errTi Te(K) errTe lg10Ne(m-3) errNe')") open (i20w1,file='son050514_15_20m_nmax') write (i20w1,"(' year dayn moda uth glat glon', | ' hmax foF2MHz Nmax(m-3) TEC(m-2)')") open (i20w2,file='son050514_15_20m_profiles') write (i20w2,"(' year dayn uth alt(km) Ti(K) ', | 'errTi Te(K) errTe lg10Ne(m-3) errNe')") 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 ipar1d(2) > 8990 (elev*100) elev = ipar1d(2)*0.01 if (elev .gt. 89.) then ihr = ibhmt/100 imn = ibhmt - ihr*100 uth = ihr + float(imn)/60. glat = ipar1d(5)*0.01 glon = ipar1d(6)*0.01 C Set altitude and Ne and fill in gaps in Ne above 100 km (jne) jne = 11 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(jne+1,j) .gt. ipar2d(jne,j)) then ipar2d(jne,j) = missing ipar2d(jne+1,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(jne+1,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,11,-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(11) .lt. 0.) then do 1615 j=12,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(11) = rnelog10(j1) - diffne*(htjm(j1)-htjm(11)) endif C write (6,"(1x,' rnelog10(11) =',e12.4)") rnelog10(11) C Fill in missing rnelog10 values as linear interpolations do 1625 j=12,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. do 1626 j=26,nrows-1 C Find max Ne from 154.6-536.8 km if (corrne(j) .gt. xnemax .and. htjm(j) .gt. 150.) then xnemax = corrne(j) jx = j endif 1626 continue h(2) = htjm(jx) p(2) = corrne(jx) hmf2 = htjm(jx) xnmax = corrne(jx) C Check to see if max Ne next to missing data or not do 1627 j=jx-1,26,-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) 1631 fof2 = sqrt(xnmax/1.24e+10) C Find tec from non-missing data tec = 0. htjm1 = htjm(10) do 1635 j=11,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 tecn = 0. do 1640 j=11,nrows fillne(j) = 10.**(rnelog10(j)*0.001) tecn = tecn + fillne(j)*(htjm(j)-htjm(j-1))*1.e+3 1640 continue write (6,"(1x,' tec,n nmax hmax =',3e12.4,f7.2)") | 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 (j=11) and up do 1675 j=11,nrows C Delete if Te, Ti error bars > 200K if (ipar2d(14,j) .gt. 200.) then ipar2d(13,j) = missing ipar2d(14,j) = missing endif if (ipar2d(16,j) .gt. 200.) then ipar2d(15,j) = missing ipar2d(16,j) = missing endif write (iw2,"(1x,i4,i5,f6.2,f7.2,4i7,2f8.3)") ibyrt,ndayn,uth, | htjm(j),ipar2d(15,j),ipar2d(16,j),ipar2d(13,j),ipar2d(14,j), | ipar2d(11,j)*0.001,ipar2d(12,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) 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 1705 j=1,nrows nav = 0 avne(j) = 0. averne(j) = 0. do 1701 n=1,nm if (a20ne(j,n) .gt. 0.) then nav = nav + 1 avne(j) = avne(j) + a20ne(j,n) averne(j) = averne(j) + a20erne(j,n) endif 1701 continue if (nav .gt. 0) then avne(j) = avne(j)/nav averne(j) = averne(j)/nav else avne(j) = missing averne(j) = missing endif nav = 0 avte(j) = 0. averte(j) = 0. do 1702 n=1,nm if (a20te(j,n) .gt. 0.) then nav = nav + 1 avte(j) = avte(j) + a20te(j,n) averte(j) = averte(j) + a20erte(j,n) endif 1702 continue if (nav .gt. 0) then avte(j) = avte(j)/nav averte(j) = averte(j)/nav else avte(j) = missing averte(j) = missing endif nav = 0 avti(j) = 0. averti(j) = 0. do 1703 n=1,nm if (a20ti(j,n) .gt. 0.) then nav = nav + 1 avti(j) = avti(j) + a20ti(j,n) averti(j) = averti(j) + a20erti(j,n) endif 1703 continue if (nav .gt. 0) then avti(j) = avti(j)/nav averti(j) = averti(j)/nav else avti(j) = missing averti(j) = missing endif if (j .gt. 10) | write (i20w2,"(1x,i4,i5,f6.2,f7.2,4f9.1,2e12.4)") ibyrt, | mdayn,avut,htjm(j),avti(j),averti(j),avte(j),averte(j), | avne(j),averne(j) 1705 continue C Calculate hmF2, NmF2, and TEC (1 TECU=10^16elec/m2) xnemax = 0. do 1726 j=26,nrows-1 C Find max Ne from 154.6-536.8 km if (avne(j) .gt. xnemax .and. htjm(j) .gt. 150.) 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,26,-1 if (avne(j) .gt. 0.) 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.) 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(10) do 1735 j=11,nrows if (avne(j) .gt. 0.) 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 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) = ipar2d(13,j) a20erte(j,nm) = ipar2d(14,j) a20ti(j,nm) = ipar2d(15,j) a20erti(j,nm) = ipar2d(16,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)) write(6,"(1x,'hnmf2 nmf2 =',f7.2,e12.4)") hmf2,nmf2 end