基于GSI读取bufr文件

张__QI

1. #! /bin/sh
   
2. #
   
3. #extract necp bufr information and data
   
4. #
   
5. 
   
6. bufr_path="/home/qzhang/2011_07_18/gdas.airsev.t00z.20110718.bufr"
   
7. exe_path="/home/qzhang/airsev"
   
8. include_path="/home/qzhang/GSI/comGSI_v3.3/include"
   
9. bufr_lib_path="/home/qzhang/GSI/comGSI_v3.3/lib"
   
10. 
    
11. mkdir $exe_path
    
12. cat << EOF > $exe_path/extract.f90
    
13. program bufr_decode_radiance
    
14. !
    
15. ! read all radaince observations out from bufr.
    
16. ! read bufr table from prepbufr file
    
17. !
    
18. implicit none
    
19. 
    
20. integer, parameter :: mxmn=35, mxlv=250
    
21. character(80):: hdstr= &
    
22.    'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS'
    
23. character(80):: hdr2b='SAZA SOZA BEARAZ SOLAZI'
    
24. character(80):: obstr='TMBR'
    
25. ! character(80):: obstr='TMBRST'
    
26. real(8) :: hdr(mxmn),hdr2(mxmn),obs(mxlv)
    
27. 
    
28. INTEGER        :: ireadmg,ireadsb
    
29. 
    
30. character(8)   :: subset
    
31. integer        :: unit_in=10,idate,nmsg,ntb
    
32. integer,parameter:: max_sat_type=20
    
33. integer        :: nsat_type(max_sat_type),nsat_num(max_sat_type)
    
34. integer        :: i,k,iret,ksatid,nchanl,num_sat_type,ii
    
35. !
    
36. nchanl=15
    
37. nsat_num=0
    
38. nsat_type=0
    
39. !
    
40. open(24,file='bufr.table')
    
41. open(unit_in,file='$bufr_path',form='unformatted',status='old')
    
42. call openbf(unit_in,'IN',unit_in)
    
43. call dxdump(unit_in,24)
    
44. call datelen(10)
    
45.    nmsg=0
    
46.    ntb = 0
    
47.    num_sat_type = 0
    
48.    msg_report: do while (ireadmg(unit_in,subset,idate) == 0)
    
49.      nmsg=nmsg+1
    
50.      write(*,*)
    
51.      write(*,'(3a,i10)') 'subset=',subset,' cycle time =',idate
    
52.      sb_report: do while (ireadsb(unit_in) == 0)
    
53.        ntb = ntb+1
    
54.        call ufbint(unit_in,hdr ,mxmn,1  ,iret,hdstr)
    
55.        call ufbint(unit_in,hdr2,mxmn,1  ,iret,hdr2b)
    
56. 
    
57.        call ufbrep(unit_in,obs ,1   ,nchanl,iret,obstr)
    
58.        ksatid=nint(hdr(1))
    
59.        write(*,*) ntb,iret,ksatid
    
60.        write(*,*),(hdr(i),i=3,5)
    
61.        write(*,*),(hdr(i),i=6,8)
    
62.        write(*,*),(hdr(i),i=9,10)
    
63.        !write(*,'(a10,15f7.1)') 'obs=',(obs(i),i=1,iret)
    
64. 
    
65. ! satellite type inventory
    
66.        if(num_sat_type == 0 ) then
    
67.           num_sat_type=1
    
68.           nsat_type(num_sat_type)=ksatid
    
69.           nsat_num(num_sat_type)= 1
    
70.        else
    
71.           ii=0
    
72.           DO i=1,num_sat_type
    
73.              if(nsat_type(i) == ksatid) ii=i
    
74.           ENDDO
    
75.           if( ii > 0 .and. ii <=num_sat_type) then
    
76.              nsat_num(ii)=nsat_num(ii) + 1
    
77.           else
    
78.              num_sat_type=num_sat_type+1
    
79.              if(num_sat_type > max_sat_type) then
    
80.                 write(*,*) 'Error: too many satellite types'
    
81.                 write(*,*) 'Need to increase max_sat_type'
    
82.                 stop 1234
    
83.              endif
    
84.              nsat_type(num_sat_type)=ksatid
    
85.              nsat_num(num_sat_type)=1
    
86.           endif
    
87.        endif
    
88. 
    
89.      enddo sb_report
    
90.    enddo msg_report
    
91. call closbf(unit_in)
    
92. 
    
93. write(*,*) 'message=',nmsg,'  subset=',ntb
    
94. DO i=1,num_sat_type
    
95.    write(*,'(i4,2i10)') i,nsat_type(i),nsat_num(i)
    
96. ENDDO
    
97. 
    
98. 
    
99. end program
    
100. EOF
     
101. 
     
102. cd $exe_path
     
103. ifort  -O2   -I $include_path -c extract.f90
     
104. ifort  -o extract.exe -O2   extract.o -L $bufr_lib_path -lbufr_i4r8
     
105. ./extract.exe &> info.data

复制代码

老板脑洞开的有点大，让我把卫星辐射数据每条记录的经纬度信息画出来看看。没办法只能照做。

献上bash脚本一个，为大家提供灌水的素材。

如果要把bufr文件里的辐射数据输出出来，只需要在脚本最内层循环里面加上这句话：write(*,*) 'obs=',(obs(i),i=1,iret)即可。

脚本地址献上：

https://github.com/smft/ncep_bufr_reader_based_on_gsi

截个图，证明 it really works！

非对称

大哥，发错地方拉，这是python板块，不是fortran板块

ginanony

的确。。。。这是fortran吧。。。。{:eb315:}{:eb315:}

张__QI

不想去fortran

ginanony

哈哈哈，这样的话，大家可以试着把这段代码实现的功能改成python的？嘿嘿，都是个学习的过程呀
{:eb301:}{:eb301:}

zx9981

支持牛人！

虫儿飞

最主要的问题是有一个Fortran版本的prepbufr库，不知道Python有没有对应版本，要有这个才行

张__QI

虫儿飞 发表于 2015-11-16 08:49
最主要的问题是有一个Fortran版本的prepbufr库，不知道Python有没有对应版本，要有这个才行

已测，没有

张__QI

虫儿飞 发表于 2015-11-16 08:49
最主要的问题是有一个Fortran版本的prepbufr库，不知道Python有没有对应版本，要有这个才行

已测，没有。。。。。。。

novelbean

感谢分享啦