program rd_rfp C 04/11: rd_rfp.f for RFP red-line FPI data 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 ahf05.asc contains RFP FPI data from Dec 14-15, 2006 open (ird,file='rfp061214_15.asc',status='old') 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 c 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 wind and Tn data from RFP (cc in tenths instead of octas or eighths) C year dayn moda uth RELTn errTn Un errUn Vn errVn Wn errWn cc glat glon dir are 2-d params C 1 2 3 4 6 7 10 11 12 13 14 15 25 (from az#23 and elev#24) write (6,"(' year dayn moda uth relTn K err Tn Un m/s err Un Vn', | ' m/s err Vn Wn m/s err Wn cc glat glon dir')") do 1675 j=1,nrows C Skip if Tn missing if (ipar2d(6,j) .lt. -32760) go to 1675 uth = ipar2d(4,j)*0.001 C Set Un,Vn,Wn to missing if cloud cover > 4/10s if (ipar2d(25,j) .gt. 4) then do 1625 i=10,15 ipar2d(i,j) = -32767 1625 continue endif wi = ipar2d(14,j)*0.1 erwi = ipar2d(15,j)*0.1 C 4 direction scan 270,90,180,0 at 45deg elevation and vertical glat = ipar1d(1)*0.01 glon = ipar1d(2)*0.01 C RFP (Resolute Bay) station at (74.73N, 94.89W) C Re = 6376 km; circumference = 4piR or glon at glat = 2piRsin(90-|glat|) C 6376*2pi/360 = 111.28 km per deg latitude at ground. C distance at 45 deg elevation angle is (flat earth): tan45 = 250/x or x=250/tan45(=1)=250 km C 250 km / 11.28km/deg = 2.25 deg in latitude C OR 250km/(sin,cos(45)=0.7071) = 353.55km or 2piR(=353.55km)*45*2/360 = 555.36 km between N-S,E-W C (or 250/sin20=730.95 km * 2pi*140/360 = 1786.05 km for 20 deg elev instead of 1174 km?) C or 555.36/(2*29.31km) = 9.474 deg lon? something off here so use flat earth! C 2piR(6376km)sin(90-glat) = 40062*cos(74.73=0.2634) = 10552 km / 360 deg = 29.31 km/deglon so 250 km C so 250 km / 29.31km/deglon = 8.53 deg lon C Z 74.73N, 94.89W C W 74.73N, 103.42W C E 74.73N, 86.36W C N 76.98N, 94.89W C S 72.48N, 94.89W if (ipar2d(24,j) .gt. 8900) then dir = ' Z' go to 1650 endif if (abs(ipar2d(23,j)*0.01-0.) .lt. 0.1) then dir = ' N' glat = ipar1d(1)*0.01 + 2.25 endif if (abs(ipar2d(23,j)*0.01-90.) .lt. 0.1) then dir = ' E' glon = ipar1d(2)*0.01 + 8.53 endif if (abs(ipar2d(23,j)*0.01-180.) .lt. 0.1) then dir = ' S' glat = ipar1d(1)*0.01 - 2.25 endif if (abs(ipar2d(23,j)*0.01-270.) .lt. 0.1) then dir = ' W' glon = ipar1d(2)*0.01 - 8.53 endif 1650 continue C year dayn moda uth RELTn errTn Un errUn Vn errVn Wn errWn cc glat glon dir are 2-d params C 1 2 3 4 6 7 10 11 12 13 14 15 25 (from az#23 and elev#24) write (6,"(1x,i4,2i5,f6.2,6i8,2f8.1,i3,2f8.2,1x,a2)") | (ipar2d(j1,j),j1=1,3),uth,ipar2d(6,j),ipar2d(7,j), | (ipar2d(n2,j),n2=10,13),wi,erwi,ipar2d(25,j),glat,glon,dir 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