PROGRAM READBUFR_TC C----------------------------------------------------------------------- C C THIS CODE WILL READ INFORMATION OUT OF NPP OMPS Totol Column BUFR OZONE FILES C AND BIN OZONE TO 1 X 1 DEG LAT/LON GRIDS C----------------------------------------------------------------------- CHARACTER SUBSET CHARACTER*8 ymd integer ksatid,syear,smonth,sday,shour,sminute,ssec real*8 slon0,slat0,rsat real*8 toz,sza,toq,viq,req real*8 hdroz(11) real*8 slat(250000),slon(250000),to3(250000) dimension toto3(360,181),ho3(1080) dimension x(364),y(364),u(360) integer nbr(360,181) DATA LUBFR/10/,IREC/0/ write(6,*) 'enter yyyymmdd:' read(5,'(a8)') ymd open(lubfr,file='/dcom/us007003/'//ymd//'/b008/xx016', .form='unformatted') open(7,file='../../omps/omps_tc_'//ymd//'.out') open(8,file='../../omps/omps_tc_'//ymd//'.b',form='unformatted') write(7,*)' year mon day hour min sec lat', &' lon toz' n0=0 n1=0 C OPEN BUFR FILE - READ IN FIRST MESSAGE C -------------------------------------------------------------- CALL DATELEN(10) CALL OPENBF(LUBFR,'IN',LUBFR) PRINT 100, LUBFR 100 FORMAT(/5X,'===> OZONE BUFR DATA SET IN UNIT',I3, $ ' SUCCESSFULLY OPENED FOR INPUT') C READ IN FIRST DATA MESSAGE AND GET DATE C --------------------------------------- CALL READMG(LUBFR,SUBSET,IDATE,IRET) IF(IRET.NE.0) GO TO 999 ! problem CALL UFBCNT(LUBFR,IREC,ISUB) PRINT 101, IDATE 101 FORMAT(/' OZONE BUFR FILE HAS DATE: ',I10///) 10 CONTINUE n0=n0+1 C READ A SUBSET (REPORT) IN MESSAGE C --------------------------------- CALL READSB(LUBFR,IRET) IF(IRET.EQ.0) GO TO 77 C ALL SUBSETS IN THIS MESSAGE PROCESSED, READ IN NEXT DATA MESSAGE C (IF ALL MESSAGES READ IN NO MORE DATA TO PROCESS) C --------------------------------------------------------------- CALL READMG(LUBFR,SUBSET,IDATE,IRET) IF(IRET.NE.0) GO TO 99 CALL UFBCNT(LUBFR,IREC,ISUB) GO TO 10 77 CONTINUE C DECODE REQUESTED INFORMATION OUT OF SUBSET (REPORT) C -------------------------------------------------------- CALL UFBINT(LUBFR,HDROZ,16,1,ILVL, $ 'SAID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SOZA ORBN') IF(ILVL.NE.1) THEN PRINT *,' *****HDROZ ARRAY NOT RETURNING 1 LEVEL - STOP 99' goto 999 END IF c print *, 'HEADER: ',n c print *, 'satellite id = ',HDROZ(1) c print *, 'year = ',HDROZ(2) c print *, 'month = ',HDROZ(3) c print *, 'day = ',HDROZ(4) c print *, 'hour = ',HDROZ(5) c print *, 'minute = ',HDROZ(6) c print *, 'second = ',HDROZ(7) c print *, 'latitude = ',HDROZ(8) c print *, 'longitude = ',HDROZ(9) c print *, 'solar zemuth = ',HDROZ(10) c print *, 'orbit number = ',HDROZ(11) rsat = hdroz(1)-191.0; ksatid=rsat slat0= hdroz(8) slon0= hdroz(9) syear= hdroz(2) smonth= hdroz(3) sday= hdroz(4) shour= hdroz(5) sminute= hdroz(6) ssec= hdroz(7) sza=hdroz(10) slat(n0)=slat0 slon(n0)=slon0 if (slon(n0) .lt. 0) slon(n0)=slon(n0)+360. CALL UFBINT(LUBFR,TOZ,1,1,ILVL,'OZON') IF(ILVL.NE.1) THEN PRINT *,' *****TOTOZ ARRAY NOT RETURNING 1 LEVEL - STOP 99' goto 999 END IF c print *, 'total ozone = ',TOZ to3(n0)=toz c CALL UFBINT(LUBFR,toq,1,1,ILVL,'OTCQF') c CALL UFBINT(LUBFR,viq,1,1,ILVL,'VIIRGQ') c CALL UFBINT(LUBFR,req,1,1,ILVL,'RETRQ') c c print *, 'quality flag= ',toq,viq,req if (toz .lt. 800 .and. toz .gt. 50) then n1=n1+1 write(7,'(i8,6i5,3f9.2)') . n1,syear,smonth,sday,shour,sminute,ssec, . slat0,slon0,toz endif GO TO 10 99 CONTINUE C ALL REPORTS HAVE BEEN DECODED AND PROCESSED C ------------------------------------------- PRINT 303, IREC 303 FORMAT(/10X,'==> ALL BUFR MESSAGES READ IN AND PROCESSED: TOTAL ', $ 'NUMBER OF DATA MESSAGES IN OZONE BUFR FILE IS:',I6//) CALL closbf(lubfr) print *, 'n0= ',n0,' n1=',n1 c c .... bin toz to 1x1 lat/lon grid c do i=1,360 do j=1,181 toto3(i,j)=0. nbr(i,j)=0 enddo enddo do n=1,n0 if (to3(n).gt.750.) goto 88 if (to3(n).lt.50.) goto 88 if (slat(n) .gt. -45. .and. to3(n).lt.150.) goto 88 i=int(slon(n)+0.5)+1 if(i .gt. 360) i=i-360 if(slat(n) .ge. 0) then j=91-int(slat(n)+0.5) else j=91-int(slat(n)-0.5) endif toto3(i,j)=toto3(i,j)+to3(n) nbr(i,j)=nbr(i,j)+1 88 continue enddo do i=1,360 do j=1,181 if(nbr(i,j) .gt. 0) then toto3(i,j)=toto3(i,j)/nbr(i,j) else toto3(i,j)=-9.999e20 endif enddo enddo c c .... drop out outliers c do j=1,181 do i=1,360 a1=0. na1=0. do k=1,61 m=i+k-31 if (m .lt. 1) m=m+360 if (m .gt. 360) m=m-360 if (toto3(m,j) .gt. 0) then a1=a1+toto3(m,j) na1=na1+1 endif enddo if (na1 .gt. 0) a1=a1/na1 if (a1 .gt. 0 .and. toto3(i,j) .gt. 0) then d1=abs(toto3(i,j)-a1) else d1=0. endif if (d1 .gt. 30.) toto3(i,j)=-9.999e20 enddo enddo write(8) toto3 c c .... interpolation of longitude gaps c do j=2,180 n=0 m=0 do i=1,360 if (toto3(i,j).gt.50 .and. toto3(i,j).lt.800) then n=n+1 x(n+2)=float(i) y(n+2)=toto3(i,j) else m=m+1 u(m)=float(i) endif enddo x(n+3)=x(3)+360 x(n+4)=x(4)+360 x(1)=x(n+1)-360 x(2)=x(n+2)-360 y(n+3)=y(3) y(n+4)=y(4) y(1)=y(n+1) y(2)=y(n+2) if (m .gt. 0 .and. m .lt. 180) then do k=1,m im=int(u(k)) call qip(n+4,x,y,u(k),f) toto3(im,j)=f enddo endif enddo c c ... latitude-dependent running mean smooth c do j=2,180 if(j .gt. 30 .and. j .lt. 152) m=3 do i=1,360 ho3(i)=toto3(i,j) ho3(i+360)=toto3(i,j) ho3(i+720)=toto3(i,j) enddo if (j .gt. 20 .or. j .lt. 161) m=5 if (j .gt. 15 .and. j .le. 20) m=9 if (j .ge. 161 .and. j .lt. 166) m=9 if (j .le. 15 .or. j .ge. 166) m=15 mh=(m-1)/2 do i=mh+1,mh+360 tmp=0. nsm=0 do k=1,m k1=k-mh-1 if (ho3(i+k1).gt.50 .and. ho3(i+k1).lt.800) then tmp=tmp+ho3(i+k1) nsm=nsm+1 endif enddo if (nsm .gt. 0) toto3(i-mh,j)=tmp/nsm enddo 333 continue enddo write(8) toto3 STOP 999 CONTINUE CALL closbf(lubfr) close(7) close(8) END subroutine qip(n,x,y,u,f) c ************************* c .... given y(n)=F(x(n)), interpolate y value(f) at x=u c dimension x(n),y(n) lin=0 nm1=n-1 do 10 i=2,nm1 if(u.le.x(i)) goto 20 10 continue i=nm1 goto 30 20 if(i.eq.2) goto 30 if(u-x(i-1).lt.x(i)-u) i=i-1 30 x1=x(i-1) x2=x(i) x3=x(i+1) d1=x2-x1 d2=x3-x2 if (d1 .gt. 3*d2 .or. d2 .gt. 3*d1) then lin=1 goto 60 endif if (x1.ne.x2 .and. x1.ne.x3 .and. x2.ne.x3) then a1=(u-x2)*(u-x3)/((x1-x2)*(x1-x3)) a2=(u-x3)*(u-x1)/((x2-x3)*(x2-x1)) a3=(u-x1)*(u-x2)/((x3-x1)*(x3-x2)) endif f=a1*y(i-1)+a2*y(i)+a3*y(i+1) if (f .lt. 0) lin=1 60 continue c c .... use linear interpolation instead c if(lin .eq. 1) then if(u.lt.x2 .and. x1.ne.x2) . f=(u-x2)/(x1-x2)*y(i-1)+(u-x1)/(x2-x1)*y(i) if(u.gt.x2 .and. x2.ne.x3) . f=(u-x3)/(x2-x3)*y(i)+(u-x2)/(x3-x2)*y(i+1) endif return end