; Grid MOPITT data for multiple days, arbitrary grid ; Reads L2 MOPITT V5 retrievals from NCAR archi ;objective is to grid the monthly co mopitt v5 data ;Heba Marey may 2012 ;modified june 2012 to be used in other pro to creat climatological maps of seasons ;------------------------------------------------------- ;------------------------------------------------------- ; Get the Vdata information from an HDF file Function get_vd,file_id,name ;Attempt to catch some errors CATCH, Error_status If Error_status NE 0 then begin Print, 'Error index: ', Error_status Print, 'Error message: ', !ERR_STRING return, -99999. endif vd_ref = hdf_vd_find(file_id, strtrim(name,2)) vdata=hdf_vd_attach(file_id,vd_ref) nread=hdf_vd_read(vdata,var) hdf_vd_detach, vdata return,var end ;------------------------------------------------------- ; Get the Scientific Data information from an HDF file Function get_sd, filename, name ;Attempt to catch some errors CATCH, Error_status If Error_status NE 0 then begin Print, 'Error index: ', Error_status Print, 'Error message: ', !ERR_STRING var = -99999. return, var endif sd_id = hdf_sd_start(filename,/read) index = hdf_sd_nametoindex(sd_id, name) sds_id = hdf_sd_select(sd_id,index) hdf_sd_getdata, sds_id, var hdf_sd_endaccess, sds_id hdf_sd_end, sd_id return, var end ;------------------------------------------------------- ; Get sw information from an HDF-EOS file function get_sw, filename, attrname ;Attempt to catch some errors CATCH, Error_status If Error_status NE 0 then begin Print, 'Error index: ', Error_status Print, 'Error message: ', !ERR_STRING var = -99999. return, var endif file_id = eos_sw_open(filename, /read) sw_id = eos_sw_attach(file_id, 'MOP02') sws_id = eos_sw_readfield(sw_id, attrname, databuffer) status = eos_sw_detach(sw_id) status = eos_sw_close(sw_id) return, databuffer end ;------------------------------------------------------- ;------------------------------------------------------- ; Grid MOPITT V5 data ; INPUTS: date_beg: start date ; date_end: end date ; ialt: altitude index ; xbin, ybin: Lon,Lat bin size ; daynt: day/night index ; aplimit: limit for % apriori ; moppath: path to MOPITT L2 data ; (e.g., '/DUFF/RR/Archive/L2/' or '/MOPITT/OPS/Archive/L2/') ; OUTPUT: latbin, lonbin, co ;------------------------------------------------------- pro grid_mop_v7_FAug23, T,match, xbin, ybin, daynt, latbin, lonbin, co,deg_fr,nblon,nblat, altlab,nco,f ; date_beg='20160601' ; date_end='20160601' ;f=5 ;ialt=-1 ;daynt=0 ;xbin=1 ;ybin=1 ;nblat = Fix(180./ybin) ;nblon = Fix(360./xbin) ;lonbin = findgen(nblon)*xbin - 180 ;new array of lon ,ex -180.,179.5,..180 ;latbin = findgen(nblat)*ybin - 90. ialt=-1 if (ialt eq 0) then altlab = 'Surface' else $ if (ialt eq 1) then altlab = '900hPa' else $ if (ialt eq 2) then altlab = '800hPa' else $ if (ialt eq 3) then altlab = '700hPa' else $ if (ialt eq 4) then altlab = '600hPa' else $ if (ialt eq 5) then altlab = '500hPa' else $ if (ialt eq 6) then altlab = '400hPa' else $ if (ialt eq 7) then altlab = '300hPa' else $ if (ialt eq 8) then altlab = '200hPa' else $ if (ialt eq 9) then altlab = '100hPa' else altlab = 'Column' if (daynt eq 0) then daylab = '(day)' else $ if (daynt eq 1) then daylab = '(night)' print,'Using values: ' ; Determine grid numbers and arrays binlabel = String(format='(f4.1,"!Eo!N Lon by",f4.1,"!Eo!N Lat")',xbin,ybin) ; define arrays for binned data co = fltarr(nblon,nblat) nco = intarr(nblon,nblat) deg_fr=fltarr(nblon,nblat) ;--------------------------------- g = 9.81 ;m/s2 mw = 28.966E-3/6.022E23 ;kg/molec of dry air units = 1.E4/100. ; [cm2/m2]/[pa/hPa] ;------------------------ ; Loop through days ;------------------------ ;-------------------------- ; find file for the day ;-------------------------- mopfiles = Findfile(match,count=nf) if (nf ge 1) then begin ;Open HDF file filename = mopfiles[0] print,filename file_id =H5F_OPEN(filename) print, 'read Solar Zenith Angle' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/SolarZenithAngle') sza = H5D_Read(dataset_id) H5D_CLOSE, dataset_id print, 'read Degrees of Freedom for Signal' ; dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/DegreesofFreedomforSignal') ;DF = H5D_Read(dataset_id) ;H5D_CLOSE, dataset_id ; read averaging kernel matrices into variable 'avkrn' print, 'read AK matrix' ; dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/RetrievalAveragingKernelMatrix') ;AKMatrix = H5D_Read(dataset_id) ;H5D_CLOSE, dataset_id print, 'Seconds in Day' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Geolocation Fields/SecondsinDay') secs = H5D_Read(dataset_id) H5D_CLOSE, dataset_id print, 'read latitudes' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Geolocation Fields/Latitude') lat_data = H5D_Read(dataset_id) H5D_CLOSE, dataset_id ; read longitudes into variable 'moplon' print, 'read longitudes' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Geolocation Fields/Longitude') lon_data = H5D_Read(dataset_id) H5D_CLOSE, dataset_id ; read cloud flags into variable 'moplon' print, 'read cloud flags' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/CloudDescription') cloud_data = H5D_Read(dataset_id) H5D_CLOSE, dataset_id if (ialt ge 1) then begin ; read retrieved profiles into variable 'rtvprofl' print, 'read retrieved CO profiles' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/RetrievedCOMixingRatioProfile') mixing_ratio = H5D_Read(dataset_id) H5D_CLOSE, dataset_id ; mixing_ratio = get_sd(filename, 'Retrieved CO Mixing Ratio Profile') if (N_elements(mixing_ratio) le 1) then goto,skipfile co_data = reform(mixing_ratio[0,ialt-1,*]) endif else if (ialt eq 0) then begin print, 'read retrieved CO surface' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/RetrievedCOSurfaceMixingRatio') mixing_ratio = H5D_Read(dataset_id) H5D_CLOSE, dataset_id ;mixing_ratio = get_sd(filename, 'Retrieved CO Surface Mixing Ratio') co_data = reform(mixing_ratio[0,*]) endif else begin print, 'Retrieved CO Total Column' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/RetrievedCOTotalColumn') co_col = H5D_Read(dataset_id) H5D_CLOSE, dataset_id ;co_col = get_sd(filename, 'Retrieved CO Total Column') if (N_elements(co_col) le 1) then goto,skipfile co_data = reform(co_col[0,*]) endelse ;Close HDF file H5F_CLOSE, file_id print,'Binning data ... ' ;lon_data runs from -180 to 180 print,'N_elements(co_data)',N_elements(co_data) print,max(nco) ;*************************************For loop***************************** for idat = 0L,N_elements(co_data)-1 do begin ; for idat = 122125L, 122130 do begin if (co_data[idat] gt 0 ) then begin ;if (co_data[idat] gt 0 ) then begin ilat = Floor((lat_data[idat] +90) /ybin) ;index of new gridded lat ilon = Floor((lon_data[idat] +180) /xbin) ;Floor round the number ;to the least if (ilon eq nblon) then ilon=0 ; the integer part,ex ;Floor(1.4 or 1.6 )is 1 ;because we want all the data in the ;grid bin to be for one index ,ex,if ;ybin is .25 , so, the data -89.75 to ;89.5 will be belong to index 1 ;DAY ; if (daynt eq 0 and sza[idat] le 80. and cloud_data[idat] eq f) then begin ;if (daynt eq 0 and sza[idat] le 80. ) then begin if (cloud_data[idat] eq f ) then begin ;convert one dim co to two dim (grided co) co[ilon,ilat] = co[ilon,ilat] + co_data[idat] nco[ilon,ilat] = nco[ilon,ilat]+1 endif $ ;NIGHT else if (daynt eq 1 and sza[idat] gt 80.) then begin co[ilon,ilat] = co[ilon,ilat] + co_data[idat] deg_fr[ilon,ilat]=deg_fr[ilon,ilat]+DF[idat] nco[ilon,ilat] = nco[ilon,ilat]+1 endif endif endfor ;**************close idat loop endif else $ print,max(nco) skipfile: ;print, 'nco',nco ;endfor ;close the day loop*************************************** end