; Grid MOPITT data for multiple days, arbitrary grid ; Reads L2 MOPITT V7 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 ;modified March to add chi2 ;------------------------------------------------------- ;------------------------------------------------------- ; 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 V7 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_all, date_beg,date_end, ialt, xbin, ybin, daynt, latbin, lonbin, deg_fr,nblon,nblat, altlab,nco,nptd,ndays pro grid_mop_v7_Aug23 ,T,match,nco,co,deg_fr,xbin,ybin,nblon,nblat,lonbin,latbin,ret_it,coError,qai2 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' ; 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 = fltarr(nblon,nblat) deg_fr=fltarr(nblon,nblat) ret_it=fltarr(nblon,nblat) coError=fltarr(nblon,nblat) qai2=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' 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 print, 'read Degrees of Retrieval iteration' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/RetrievalIterations') RI = H5D_Read(dataset_id) H5D_CLOSE, dataset_id print,'signalCh2' dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/SignalChi2') chi2 = H5D_Read(dataset_id) H5D_CLOSE, dataset_id if (ialt ge 1) then begin ; read retrieved profiles into variable 'rtvprofl' 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 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 dataset_id = H5D_OPEN(file_id, '/HDFEOS/SWATHS/MOP02/Data Fields/RetrievedCOTotalColumn') co_col = H5D_Read(dataset_id) H5D_CLOSE, dataset_id if (N_elements(co_col) le 1) then goto,skipfile co_data = reform(co_col[0,*]) co_error=reform(co_col[1,*]) print,'total co column' help,co_data 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) ;*************************************For loop***************************** for idat = 0L,N_elements(co_data)-1 do begin ;for idat = 22200L,22208 do 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 co[ilon,ilat] = co[ilon,ilat] + co_data[idat] nco[ilon,ilat] = nco[ilon,ilat]+1 deg_fr[ilon,ilat]=deg_fr[ilon,ilat]+DF[idat] coError[ilon,ilat]=coError[ilon,ilat]+co_error[idat] ret_it[ilon,ilat]=ret_it[ilon,ilat]+ RI[idat] qai2[ilon,ilat]=qai2[ilon,ilat]+ chi2[idat] ;qai square endif endfor ;**************close idat loop endif else $ print,'no MOPITT file: ',match skipfile: print,max(nco) for ilon=0,n_elements(lonbin)-1 do begin for ilat=0,n_elements(latbin)-1 do begin if nco(ilon,ilat) gt 0. then begin deg_fr(ilon,ilat)=deg_fr(ilon,ilat)/nco(ilon,ilat) coError(ilon,ilat)=coError(ilon,ilat)/nco(ilon,ilat) ret_it(ilon,ilat)=ret_it(ilon,ilat)/nco(ilon,ilat) qai2(ilon,ilat)=qai2(ilon,ilat)/nco(ilon,ilat) endif endfor endfor print,'max',max(deg_fr),max(coError),max(ret_it) end