读TRMM的3B43和3B42可以使用相同的程序.doc

读TRMM的3B43和3B42可以使用相同的程序，进行时间统计计算时注意它们是3小时和日值数据就行。附件3B43文档中有以IDL语言编写的读3B43的HDF文件的程序，用到3B42的HDF文件是没有问题的。建议先装IDL语言（软件），或ENVI（基于IDL平台，因此有IDL环境），然后把数据输出到TXT文件中（ASCII码格式），再用GRADS画图。或用IDL进行数据处理、画图，均可。SATE_TRMM_3B43_CASE.pro程序文件pro SATE_TRMM_3B43_CASE; 1*1440*400FILENAME=SATE_L3_TRM_MUTDS_MWB_3B43_GLB_V6-199807.HDFprint,*,FILENAMESD_ID=HDF_SD_START(FILENAME,/READ)P=precipitation ; range:0-10, unit:mm hr-1SDS_INDEX_P=HDF_SD_NAMETOINDEX(SD_ID,P)SDS_ID_P=HDF_SD_SELECT(SD_ID,SDS_INDEX_P)HDF_SD_GETDATA,SDS_ID_P,DATA_PHDF_SD_ENDACCESS,SDS_ID_PHDF_SD_END,SD_IDPRINT,DATA_P(0,0)print,size(DATA_P)OPENW,LUN,0_50N_0_180E_199807_output.txt,/GET_LUN ; 0-50N,0-180E east asiaFOR J=0,399 DO BEGINFOR I=0,1439 DO BEGINIF J gt 200 THEN BEGIN ;0-50N (J200) IF I gt 720 THEN BEGIN ;0-180E (1440I -1800, America RINTF,LUN,(J+1)/4.0-50.125,(I+1)/4.0-180.125,DATA_P(J,I,0), $ ;format=(2(f6.2,2x),f6.4,2x);ENDELSE;1-1440, I, - 0.125-359.875, Longitude;1-720 - (Europea, America), 0.125-179.875W, Lon=I/4.0-180.125 (Asia, Australia)0.125-179.875E, Lon=I/4.0-180.12549.875N-49.875S, Latitude;1-200 -Lat=J/4.0-50.125 Lat=J/4.0-50.1250, 0.125N-49.875NENDFORENDFORCLOSE,LUNprint,endEND附件中文档所带的FORTRAN程序未用过。但下面的程序读BIN格式的3B4X数据应该没问题。PROGRAM TRMMreal*4 data1(1440,400)real data(1440,400)real lon(1440)real lat(400)open(10,file=3B42.000101.00z.6.precipitation.bin,+ form=unformatted,+ access=DIRECT,status=OLD,recl=1440*400*4)open(12,file=rain.out,form=unformatted,+ access=DIRECT,status=unknown,recl=1440*400*4)cread(10,rec=1)(data(ii,jj),ii=1,1440),jj=1,400)cdo 10 ii=1,1440do 10 jj=1,400lon(ii) = -179.875+0.25*(ii-1)lat(jj) = -49.875+0.25*(jj-1)C write(*,*)lon(ii),lat(jj),data(ii,jj)data1(ii,jj)=data(ii,jj)10 continuewrite(12,rec=1)(data1(ii,jj),ii=1,1440),jj=1,400)cclose(10)cend=pro Main_ps, year, month, day, cs_granule infile_cloudsat = 2013290031616_39743_CS_2B-CWC-RO_GRANULE_P_R04_E06.hdf pattern = F:thesisfrancisco+infile_cloudsatread_Cloudsat_2B_CWC, pattern, Profile_time, UTC_start, TAI_start, Latitude, Longitude, Height, $ Data_quality, Data_status,RO_ice_effective_radius,RO_ice_number_conc,RO_ice_water_content,RO_liq_effective_radius,RO_liq_number_conc,RO_liq_water_content,$ Data_targetID,DEM_elevation,FlagCounts,BinCounts,Status;= max_Cloudsat = n_elements(Latitude)-1 R_2 = dblarr(max_Cloudsat,125) NB_2 = dblarr(max_Cloudsat,125) CON_2 = dblarr(max_Cloudsat,125) LER_2 = dblarr(max_Cloudsat,125) LNC_2 = dblarr(max_Cloudsat,125) LWC_2 = dblarr(max_Cloudsat,125) Height_2 = dblarr(max_Cloudsat,125) for i=0.0,124.0 do begin for j=0.0,max_Cloudsat-1 do begin R_2j,i =RO_Ice_Effective_Radiusi,j NB_2j,i = RO_Ice_Number_Conci,j CON_2j,i = RO_Ice_Water_Contenti,j LNC_2j,i= RO_Liq_Number_Conci,j LER_2j,i= RO_Liq_Effective_Radiusi,j LWC_2j,i= RO_Liq_Water_Contenti,j Height_2j,i = Heighti,j endfor endfornn=1139;524.0 for a =10334.0, max_Cloudsat-1.0, nn do begin;7684 is lat 16, 7684+1140 is lat 26 for ll = 0, 124 do begin R_2_temp = fltarr(n_elements(R_2*,ll) NB_2_temp = fltarr(n_elements(R_2*,ll) CON_2_temp = fltarr(n_elements(R_2*,ll) LER_2_temp = fltarr(n_elements(R_2*,ll) LNC_2_temp = fltarr(n_elements(R_2*,ll) LWC_2_temp = fltarr(n_elements(R_2*,ll) R_2_temp = r_2*,ll/10. NB_2_temp = nb_2*,ll/10. CON_2_temp = con_2*,ll LER_2_temp = Ler_2*,ll/10. LNC_2_temp = Lnc_2*,ll/10. LWC_2_temp = Lwc_2*,ll ; Mask out all the missing value, which is less than 0.0 index = where( (R_2_temp le 0.0) ) if(index0 ne -1) then R_2_tempindex = 0.0 R_2a:a+nn,ll = R_2_tempa:a+nn; Mask out all the missing value, which is less than 0.0 index = where( (NB_2_temp le 0.0) ) if(index0 ne -1) then NB_2_tempindex = 0.0 NB_2a:a+nn,ll = NB_2_tempa:a+nn; Mask out all the missing value, which is less than 0.0 index = where( (CON_2_temp le 0.0) ) if(index0 ne -1) then CON_2_tempindex = 0.0 CON_2a:a+nn,ll = CON_2_tempa:a+nn; Mask out all the missing value, which is less than 0.0 index = where( (LER_2_temp le 0.0) ) if(index0 ne -1) then LER_2_tempindex = 0.0 LER_2a:a+nn,ll = LER_2_tempa:a+nn ; Mask out all the missing value, which is less than 0.0 index = where( (LNC_2_temp LE 0.0) if(index0 ne -1) then LNC_2_tempindex = 0.0 LNC_2a:a+nn,ll = LNC_2_tempa:a+nn; Mask out all the missing value, which is less than 0.0 index = where( (LWC_2_temp le 0.0) ) if(index0 ne -1) then LWC_2_tempindex = 0.0 LWC_2a:a+nn,ll = LWC_2_tempa:a+nnendforH=fltarr(nn+1,125)H0:nn,* = Height_2a:a+nn,*/1000. print,max(r_2) print,max(nb_2) print,max(con_2) print,max(ler_2) print,max(lnc_2) print,max(lwc_2) print,r_2(11235,68) set_plot,ps Device, filename=F:thesisfranciscomain.ps,/color,bits_per_pixel=8;,xsize=15,ysize=10 max_value=200 levels=INDGEN(50)*4;+1;there r 50 colors here,so need to change the value of 4 or color number for different max_value nlevels=n_elements(levels) ncolors=nlevels+1 c_levels=levels c_colors=INDGEN(ncolors)+1;+3 LoadCT,33,ncolors =ncolors,file=F:thesisfranciscoCloudsat.tbl tvlct,rr,gg,bb,/get rgb = bytarr(256,3) rgb*,0=rr* rgb*,1=gg* rgb*,2=bb* iblack = 0 rgb(iblack,*) = 0,0,0 iwhite = 1 rgb(iwhite,*) = 255,255,255 tvlct, rgb0:49,* nct00=Latitudea+nn-Latitudea nct1=nct00/6.0 xstr=strarr(7) nct000=Longitudea+nn-Longitudea nct2=nct000/6.0 xstr=strarr(7) for i=0,6 do begin xstri = strtrim(string(Latitudea+i*nct1, form=(f7.2), 2) +N +!c + $ strtrim(string(Longitudea+i*nct2, form=(f7.2), 2) + E endforPRINT,xstr ; contour, ler_2a:a+nn,*, Latitudea:a+nn, H0:nn,*, /Fill, c_colors=c_colors, levels=levels,$ ; background=255,nlevels=nlevels, position=0.1,0.15,0.9,0.5,xstyle=1,thick=6.0,xrange=6,15,yrange=0,20,$ ; xticks=6,title=Ice_Effective_Radius (um!u !n),ytitle=Altitude(Km),CHARSIZE=1.2,charthick=2.0,$ ; xminor=5,xthick=6.0,ythick=6.0 ,/norm,xtickname=xstr ; colorbar, ncolors =ncolors,bottom=0, divisions=5, range=0,max_value,CHARSIZE=0.9,charthick=2.0,format=(i4),vertical=1,right=1,$ ; position=0.93,0.15,0.98,0.5,/norm ; XYoutS,0.48,0.45,/normal,Eye,charsize=1.2,Color=iblack,charthick=2.0 ; xyouts,0.45,0.03,Lat|Lon,charsize=1.3,charthick=2.0,/norm, color=iblack;0 contour, r_2a:a+nn,*, Latitudea:a+nn, H0:nn,*, /Fill, c_colors=c_colors, levels=levels,$ background=255,nlevels=nlevels, position=0.1,0.15,0.9,0.5,xstyle=1,thick=6.0,xrange=6,16,yrange=0,20,$ xticks=6,title=RO_Ice_Effective_Radius (um!u !n),ytitle=Altitude(Km),CHARSIZE=1.2,charthick=2.0,$ xminor=5,xthick=6.0,ythick=6.0 ,/norm,xtickname=xstr colorbar, ncolors =ncolors,bottom=0, divisions=5, range=0,max_value,CHARSIZE=0.9,charthick=2.0,format=(i4),vertical=1,right=1,$ position=0.93,0.15,0.98,0.5,/norm XYoutS,0.45,0.47,/normal,Eye,charsize=1,Color=iblack,charthick=2.0 xyouts,0.45,0.03,Lat|Lon,charsize=1.3,charthick=2.0,/norm, color=iblack;0 ;RO_Ice_Effective_Radius=max_value=800 levels=INDGEN(50)*16;+1;there r 50 colors here,so need to change the value of 4 or color number for different max_value nlevels=n_elements(levels) ncolors=nlevels+1 c_levels=levels c_colors=INDGEN(ncolors)+1;+3 LoadCT,33,ncolors =ncolors,file=F:thesisfranciscoCloudsat.tbl tvlct,rr,gg,bb,/get rgb = bytarr(256,3) rgb*,0=rr* rgb*,1=gg* rgb*,2=bb* iblack = 0 rgb(iblack,*) = 0,0,0 iwhite = 1 rgb(iwhite,*) = 255,255,255 tvlct, rgb0:49,* contour, nb_2a:a+nn,*, Latitudea:a+nn, H0:nn,*, /Fill, c_colors=c_colors, levels=levels,$ background=255,nlevels=nlevels, position=0.1,0.15,0.9,0.5,xstyle=1,thick=6.0,xrange=6,16,yrange=0,20,$ xticks=6,title=RO_Ice_Number_Conc (L!u-1!n),ytitle=Altitude(Km),CHARSIZE=1.2,charthick=2.0,$ xminor=5,xthick=6.0,ythick=6.0 ,/norm,xtickname=xstr colorbar, ncolors =ncolors,bottom=0, divisions=5, range=0,max_value,CHARSIZE=0.9,charthick=2.0,format=(i4),vertical=1,right=1,$ position=0.93,0.15,0.98,0.5,/normXYoutS,0.45,0.47,/normal,Eye,charsize=1,Color=iblack,charthick=2.0 xyouts,0.45,0.03,Lat|Lon,charsize=1.3,charthick=2.0,/norm, color=iblack;0;RO_Ice_Number_Conc= max_value=1600 levels=INDGEN(50)*32;+1;there r 50 colors here,so need to change the value of 4 or color number for different max_value nlevels=n_elements(levels) ncolors=nlevels+1 c_levels=levels c_colors=INDGEN(ncolors)+1;+3 LoadCT,33,ncolors =ncolors,file=F:thesisfranciscoCloudsat.tbl tvlct,rr,gg,bb,/get rgb = bytarr(256,3) rgb*,0=rr* rgb*,1=gg* rgb*,2=bb* iblack = 0 rgb(iblack,*) = 0,0,0 iwhite = 1 rgb(iwhite,*) = 255,255,255 tvlct, rgb0:49,* contour, con_2a:a+nn,*, Latitudea:a+nn, H0:nn,*, /Fill, c_colors=c_colors, levels=levels,$ background=255,nlevels=nlevels, position=0.1,0.15,0.9,0.5,xstyle=1,thick=6.0,xrange=6,16,yrange=0,20,$ xticks=6,title=RO_Ice_Water_Content (mg*m!u-3!n),ytitle=Altitude(Km),CHARSIZE=1.2,charthick=2.0,$ xminor=5,xthick=6.0,ythick=6.0 ,/norm,xtickname=xstr colorbar, ncolors =ncolors,bottom=0, divisions=5, range=0,max_value,CHARSIZE=0.9,charthick=2.0,format=(i4),vertical=1,right=1,$ position=0.93,0.15,0.98,0.5,/normXYoutS,0.45,0.47,/normal,Eye,charsize=1,Color=iblack,charthick=2.0 xyouts,0.45,0.03,Lat|Lon,charsize=1.3,charthick=2.0,/norm, color=iblack;0 ;RO_Ice_Water_Content=max_value=100 levels=INDGEN(50)*2;+1;there r 50 colors here,so need to change the value of 4 or color number for different max_value nlevels=n_elements(levels) ncolors=nlevels+1 c_levels=levels c_colors=INDGEN(ncolors)+1;+3 LoadCT,33,ncolors =ncolors,file=F:thesisfranciscoCloudsat.tbl tvlct,rr,gg,bb,/get rgb = bytarr(256,3) rgb*,0=rr* rgb*,1=gg* rgb*,2=bb* iblack = 0 rgb(iblack,*) = 0,0,0 iwhite = 1 rgb(iwhite,*) = 255,255,255 tvlct, rgb0:49,* contour, lnc_2a:a+nn,*, Latitudea:a+nn, H0:nn,*, /Fill, c_colors=c_colors, levels=levels,$ background=255,nlevels=nlevels, position=0.1,0.15,0.9,0.5,xstyle=1,thick=6.0,xrange=6,16,yrange=0,20,$ xticks=6,title=RO_Liq_Number_Conc (L!u-1!n),ytitle=Altitude(Km),CHARSIZE=1.2,charthick=2.0,$ xminor=5,xthick=6.0,ythick=6.0 ,/norm,xtickname=xstr colorbar, ncolors =ncolors,bottom=0, divisions=5, range=0,max_value,CHARSIZE=0.9,charthick=2.0,format=(i4),vertical=1,right=1,$ position=0.93,0.15,0.98,0.5,/normXYoutS,0.45,0.47,/normal,Eye,charsize=1,Color=iblack,charthick=2.0 xyouts,0.45,0.03,Lat|Lon,charsize=1.3,charthick=2.0,/norm, color=iblack;0 ;RO_Liq_Number_Conc=max_value=25 levels=INDGEN(50)*0.5;+1;there r 50 colors here,so need to change the value of 4 or color number for different max_value nlevels=n_elements(levels) ncolors=nlevels+1 c_levels=levels c_colors=INDGEN(ncolors)+1;+3 LoadCT,33,ncolors =ncolors,file=F:thesisfranciscoCloudsat.tbl tvlct,rr,gg,bb,/get rgb = bytarr(256,3) rgb*,0=rr* rgb*,1=gg* rgb*,2=bb* iblack = 0 rgb(iblack,*) = 0,0,0 iwhite = 1 rgb(iwhite,*) = 255,255,255 tvlct, rgb0:49,* contour, ler_2a:a+nn,*, Latitudea:a+nn, H0:nn,*, /Fill, c_colors=c_colors, levels=levels,$ background=255,nlevels=nlevels, position=0.1,0.15,0.9,0.5,xstyle=1,thick=6.0,xrange=6,16,yrange=0,20,$ xticks=6,title=RO_Liq_Effective_Radius (um!u !n),ytitle=Altitude(Km),CHARSIZE=1.2,charthick=2.0,$ xminor=5,xthick=6.0,ythick=6.0 ,/norm,xtickname=xstr colorbar, ncolors =ncolors,bottom=0, divisions=5, range=0,max_value,CHARSIZE=0.9,charthick=2.0,format=(i4),vertical=1,right=1,$ position=0.93,0.15,0.98,0.5,/normXYoutS,0.45,0.47,/normal,Eye,charsize=1,Color=iblack,charthick=2.0 xyouts,0.45,0.03,La