program rd_ahf C 03/11: rd_ahf.f for AHF 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),moda1(2),moda2(2) character*2 dir ird = 11 C ahf05.asc contains 2005 AHF FPI data from Jul 9-11 and Aug 31 - Sep 01 only open (ird,file='ahf05.asc',status='old') c moda1(1) = 709 c moda2(1) = 711 c moda1(2) = 831 c moda2(2) = 901 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 AHF write (6,"(' year dayn moda uth Tn (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(2,j) .lt. -32760) go to 1675 uth = ipar2d(1,j)*0.001 C Set Un,Vn,Wn to missing if cloud cover > 4 if (ipar2d(mpar,j) .gt. 4) then do 1625 i=6,11 ipar2d(i,j) = -32767 1625 continue endif wi = ipar2d(10,j)*0.1 erwi = ipar2d(11,j)*0.1 C 9 direction scan 0,315,270,225,180,135,90,45 at 20deg elevation and vertical glat = ipar1d(3)*0.01 glon = ipar1d(4)*0.01 C Re = 6378 km + 250 km = 6628 km; circumference = 4piR or glon at glat = 2piRsin(90-|glat|) C distance at 20 deg elevation angle is 587 km from station at (77.83S,166.66E) C 6378*2pi=39748/360 = 110 km per deg latitude at ground or 587/110 = 5.33 deg lat (like 5.275 deg) C 5.275 deg implies 111.28 km per deg lat or R of 6376 (so is lat on ground, not at 250 km) C 2piRsin(90-glat) = 40062*cos(72.535,77.83,83.105 = .3001,.2108,.1075) = 12022.6,8445.1,4306.7 C / 360 deg = 33.396km/dlon,23.459km/dlon,11.96 km/dlon or in 587 km is 17.58,25.02,49.08dlon C cos(45) = 0.7071 *587km = 415.07km /111.07 km = 3.74 deg lat C 2piResin(90-glat) or 2piRecos(glat). For glat=77.8296S,+3.74=81.57,+5.23=83.06,-3.74=74.09,-5.23=72.60, C cos(glat) = 0.2198,0.1466,0.1208,0.2741,0.2990, so 2piRe(=6376,40,061.6)/360deg(=111.07km/deg)*cos(glat) C = 24.41km/dlon, 16.28km/dlon, 13.42km/dlon, 30.44km/dlon, 33.21km/dlon so 415.07 at 81.57S and 74.09S is C 25.50dlon and 21.53dlon while 587 km at 77.83S is 24.05dlon. The 9 locations become glat/glon pairs of: C N 72.60S,166.67E C NW 74.09S,145.14E C W 77.83S,142.62E C SW 81.57S,141.17E C S 83.06S,166.67E C Z 77.83S,166.67E C SE 81.57S,192.72E C E 77.83S,190.72E C NE 74.09S,188.20E if (ipar2d(16,j) .eq. 9000) then dir = ' Z' go to 1650 endif if (ipar2d(15,j) .eq. 0) then dir = ' N' glat = -72.60 endif if (ipar2d(15,j) .eq. 4500) then dir = 'NE' glat = -74.09 glon = 188.20 endif if (ipar2d(15,j) .eq. 9000) then dir = ' E' glon = 190.72 endif if (ipar2d(15,j) .eq. 13500) then dir = 'SE' glat = -81.57 glon = 192.72 endif if (ipar2d(15,j) .eq. 18000) then dir = ' S' glat = -83.06 endif if (ipar2d(15,j) .eq. 22500) then dir = 'SW' glat = -81.57 glon = 141.17 endif if (ipar2d(15,j) .eq. 27000) then dir = ' W' glon = 142.62 endif if (ipar2d(15,j) .eq. 31500) then dir = 'NW' glat = 74.09 glon = 145.14 endif 1650 continue write (6,"(1x,i4,2i5,f6.2,6i8,2f8.1,i3,2f8.2,1x,a2)") | ipar1d(1),ipar1d(2),ibdtt,uth,ipar2d(2,j),ipar2d(3,j), | (ipar2d(n2,j),n2=6,9),wi,erwi,ipar2d(mpar,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