program rd_mlh13300 C 04/11: rd_mlh13300.f for Millstone Hill ISR gridded vertical 15-min data C deduced from other elevation angles using B splines etc for GEM-CEDAR Challenge 2011 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) character*2 dir 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 mlh050831.101 and the reduced version mlh050901.asc contain MLH ISR data from Sep 1, 2005 C (mlh050901.asc deleted the initial prior guess (MSIS) Tn only from MSIS from 0000-0700 UT) C mlh061214_15.asc contains MLH ISR data from Dec 14-15, 2006 open (ird,file='mlh050901.asc',status='old') open (iw1,file='mlh050901_nmax') open (iw2,file='mlh050901_profiles') C open (ird,file='mlh061214_15.asc',status='old') C open (iw1,file='mlh061214_15_nmax') C open (iw2,file='mlh061214_15_profiles') write (iw1,"(' year dayn moda uth glat glon', | ' hmax foF2MHz Nmax(m-3) TEC(m-2)')") write (iw2,"(' year dayn uth alt(km) Tn(K) err_Tn Ti(K) ', | 'err_Ti Te(K) err_Te Ne(m-3) err_Ne')") 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 Print out Ne(uncorr), and calc Ne(corr)~Ne(uncorr)*(1+Te/Ti)/2, Te, Ti, and Tn data C from MLH ISR at all heights in one file, and TEC, hmax, Nmax, and foF2 in another file C write (iw1,"(' year dayn moda uth glat glon ', C | ' hmax(km) foF2 (MHz) Nmax(m-3) TEC(m-2)')") C Calculate Nmax=1.24e+10*foF2*foF2 ihr = ibhmt/100 imn = ibhmt - ihr*100 uth = ihr + float(imn)/60. glat = ipar1d(1)*0.01 glon = ipar1d(2)*0.01 fof2 = ipar1d(4)*0.01 tec = ipar1d(3)*1.e+15 xnmax = 1.24e+10*fof2*fof2 write (iw1,"(1x,i4,2i5,f6.2,2f7.2,i4,f6.2,2e12.4)") ibyrt, | ndayn,ibdtt,uth,glat,glon,ipar1d(5),fof2,xnmax,tec C write(iw2,"(' year dayn uth alt(km) Tn(K) Ti(K) Te(K) Ne(m-3)')") do 1675 j=1,nrows C Delete <=108 km if values of Te, Ti, or Ne > values at 112 km (j=4) do 1630 n=3,7,2 if (j .lt. 4 .and. ipar2d(n,j).gt.0. .and. ipar2d(n,4).gt.0. | .and. ipar2d(n,j).gt.ipar2d(n,4)) then do 1625 m=3,8 ipar2d(m,j) = -32767. 1625 continue endif 1630 continue C Delete if Te, Ti or Ne error bars > values do 1640 n=3,7,2 if (ipar2d(n+1,j).gt.ipar2d(n,j)) then do 1635 m=3,8 ipar2d(m,j) = -32767. 1635 continue endif 1640 continue C Find corrected Ne from uncorrected if have Te and Ti itn = ipar2d(9,j) iti = ipar2d(5,j) ite = ipar2d(7,j) itier = ipar2d(6,j) iteer = ipar2d(8,j) corrne = -32767. ercorrne = -32767. ercorrne = -32767. if (ipar2d(3,j) .gt. 0) then uncorrne = 10.**(ipar2d(3,j)*0.001) eruncorrne = 10.**(ipar2d(4,j)*0.001) if (iti .gt. 0 .and. ite .gt. 0) then corrne = uncorrne * (1. + float(ite)/float(iti))/2. ercorrne = eruncorrne * (1. + float(ite)/float(iti))/2. endif endif write (iw2,"(1x,i4,i5,f6.2,i4,5i7,2e12.4)") ibyrt,ndayn, | uth,ipar2d(1,j),itn,iti,itier,ite,iteer,corrne,ercorrne 1675 continue 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