;$Id$ ;Variable Name = passive_pixel_element ;Variable Name = closest_calipso_line ;Variable Name = closest_calipso_latitude ;Variable Name = closest_calipso_longitude ;Variable Name = closest_calipso_time ;Variable Name = closest_calipso_dem_elev ;Variable Name = closest_calipso_time_difference ;Variable Name = closest_calipso_distance ;Variable Name = closest_calipso_number_of_cloud_layers ;Variable Name = closest_calipso_cloud_fraction_x3 ;Variable Name = closest_calipso_cloud_fraction_x5 ;Variable Name = closest_calipso_cloud_fraction_x7 ;Variable Name = closest_calipso_cloud_mid_temperature ;Variable Name = closest_calipso_cod ;Variable Name = closest_calipso_multilayer_flag ;Variable Name = closest_calipso_top_alt_eff_geo ;Variable Name = closest_calipso_top_alt_eff_opd ;Variable Name = closest_calipso_top_pres_eff_geo ;Variable Name = closest_calipso_tropopause_height ;Variable Name = closest_calipso_cloud_phase_profile ;Variable Name = closest_calipso_cloud_phase_top_layer ;Variable Name = closest_calipso_aerosol_profile ;Variable Name = closest_calipso_cloud_top_profile ;Variable Name = closest_calipso_cloud_base_profile ;Variable Name = closest_calipso_cloud_top_pres_profile ;Variable Name = closest_calipso_cloud_base_pres_profile pro map_calipso input_path = "./input/" dlon = 10.0 dlat = 10.0 lat_min = -90.0 lat_max = 90.0 lon_min = -180.0 lon_max = 180.0 nlon = (lon_max-lon_min)/dlon nlat = (lat_max-lat_min)/dlat nmonth = 12 ntimes = 3 missing_input = -8888.0 missing_output = -999.0 cod_thresh = 0.1 count_min = 10.0 count_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) cld_frac_sum_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) cld_frac_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE=missing_output) emiss_sum_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) emiss_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE=missing_output) multi_sum_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) multi_frac_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE=missing_output) files = file_search(input_path,'*.hdf') n_files = n_elements(files) missing_input = -8888.0 for i_file = 0, n_files - 1 do begin temp = strsplit(files[i_file], '.', /extract) year = fix(strmid(temp[1], 1,4)) doy = fix(strmid(temp[1], 5,3)) compute_month_day, doy, leapyr(year), month, day ;--- read calipso data read_hdf, files[i_file],'passive_pixel_element', intersection_pixel_element, error_status if (error_status ne 0) then begin print, 'error reading this matchup file' goto, cycle endif read_hdf, files[i_file],'closest_calipso_number_of_cloud_layers', num_cld_layers_1d, error_status read_hdf, files[i_file],'closest_calipso_cod', cod_1d, error_status read_hdf, files[i_file],'closest_calipso_latitude', lat_1d, error_status read_hdf, files[i_file],'closest_calipso_longitude', lon_1d, error_status read_hdf, files[i_file],'closest_calipso_time', time_1d, error_status read_hdf, files[i_file],'closest_calipso_cloud_phase_top_layer', phase_1d, error_status read_hdf, files[i_file],'closest_calipso_multilayer_flag', multi_1d, error_status n = n_elements(num_cld_layers_1d) for m = 0, n - 1 do begin if (cod_1d[m] lt 0.0) then goto, skip i = (lon_1d[m] - lon_min) / dlon j = (lat_1d[m] - lat_min) / dlat k = month - 1 local_time = time_1d[m] + lon_1d[m] / 15.0 if (local_time lt 0) then local_time = local_time + 24.0 if (local_time gt 24) then local_time = local_time - 24.0 d = -1 if (local_time gt 6 and local_time lt 18) then d = 1 ;day if (local_time lt 6 or local_time gt 18) then d = 2 ;night ;print, time_1d[m], lon_1d[m], local_time, d count_table[i,j,k,0] = count_table[i,j,k,0] + 1.0 count_table[i,j,k,d] = count_table[i,j,k,d] + 1.0 if (cod_1d[m] lt cod_thresh) then goto, skip cld_frac_sum_table[i,j,k,0] = cld_frac_sum_table[i,j,k,0] + 1.0 multi_sum_table[i,j,k,0] = multi_sum_table[i,j,k,0] + multi_1d[m] cld_frac_sum_table[i,j,k,d] = cld_frac_sum_table[i,j,k,d] + 1.0 multi_sum_table[i,j,k,d] = multi_sum_table[i,j,k,d] + multi_1d[m] skip: endfor cycle: endfor idx = where(count_table gt count_min,cc) if (cc gt 0) then begin cld_frac_table[idx] = cld_frac_sum_table[idx] / count_table[idx] multi_frac_table[idx] = multi_sum_table[idx] / count_table[idx] endif cld_frac_table_smoothed = smooth(cld_frac_table,[3,3,2,1],/Edge_Mirror,Missing=missing_output) multi_frac_table_smoothed = smooth(multi_frac_table,[3,3,1,1],/Edge_Mirror,Missing=missing_output) window, 0, xsize = 800, ysize = 600 load_color_table, 33, 0, 1.0, 249 erase, color = 3 map_set, 0.0, 0.0, limit=[-90,-180,90,180], position = [0.05,0.15,0.95,0.95], /noerase view, reform(cld_frac_table[*,*,0,2]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Night Total Cloud Fraction', /current map_continents, color = 3, thick = 2, /lores stop view, reform(cld_frac_table[*,*,0,0]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Total Cloud Fraction' view, reform(cld_frac_table[*,*,0,1]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day Total Cloud Fraction' view, reform(cld_frac_table[*,*,0,2]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Night Total Cloud Fraction' view, reform(multi_frac_table[*,*,0,0]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Multilayer Fraction' view, reform(multi_frac_table[*,*,0,1]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day Multilayer Fraction' view, reform(multi_frac_table[*,*,0,2]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Night Multilayer Fraction' ;--------------------------------- nc_filename = 'nb_cloud_mask_calipso_prior.nc' header = 'Andy Heidinger made this' sd_id = ncdf_create(nc_filename,/CLOBBER) ncdf_attput,sd_id, "description", header, /CHAR, /global ncdf_attput,sd_id, "number_longitudes", float(nlon), /FLOAT, /global ncdf_attput,sd_id, "number_latitudes", float(nlat), /FLOAT, /global ncdf_attput,sd_id, "number_months", float(nmonth), /FLOAT, /global ncdf_attput,sd_id, "number_times", float(ntimes), /FLOAT, /global ncdf_attput,sd_id, "latitude_min", float(lat_min), /FLOAT, /global ncdf_attput,sd_id, "latitude_max", float(lat_max), /FLOAT, /global ncdf_attput,sd_id, "longitude_min", float(lon_min), /FLOAT, /global ncdf_attput,sd_id, "longitude_max", float(lon_max), /FLOAT, /global ncdf_attput,sd_id, "longitude_spacing", float(dlon), /FLOAT, /global ncdf_attput,sd_id, "latitude_spacing", float(dlat), /FLOAT, /global ;--- get dimensions nx_id = NCDF_DIMDEF(sd_id, 'longitude', nlon) ny_id = NCDF_DIMDEF(sd_id, 'latitude', nlat) nm_id = NCDF_DIMDEF(sd_id, 'month', nmonth) nd_id = NCDF_DIMDEF(sd_id, 'diurnal', ntimes) ;--- create variable sds_id1 = NCDF_VARDEF(sd_id, 'cloud_fraction_table', [nx_id,ny_id,nm_id,nd_id], /FLOAT) NCDF_CONTROL, sd_id, /ENDEF ; Put the file into data mode ;---------------------------------------------- ; write to netcdf file ;---------------------------------------------- NCDF_VARPUT, sd_id, sds_id1, cld_frac_table ;--- close netcdf files ncdf_close, sd_id stop end ;=================================================================== ; routine to read hdf sds's ;=================================================================== pro read_hdf, filename, sds_name, sds_data, Error_Status catch, Error_Status if (Error_Status ne 0) then begin print, "error in read_hdf, no ",sds_name sds_data = -1 return endif sd_id = hdf_sd_start(filename, /READ) sds_id = hdf_sd_select(sd_id, hdf_sd_nametoindex(sd_id,sds_name)) hdf_sd_getdata, sds_id, sds_data hdf_sd_endaccess, sds_id hdf_sd_end, sd_id catch, /CANCEL end