;$Id$
pro sc_prob_calipso_prior

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_clear_thresh = 0.1
cod_ice_water_thresh = 0.1
count_min = 10.0

count_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0)
sc_cld_frac_sum_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0)
sc_cld_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

print, 'doing ', i_file, "of ", n_files

   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_5km', 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_cod_ice', codi_1d, error_status
   read_hdf, files[i_file],'closest_calipso_cod_water', codw_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
   read_hdf, files[i_file],'closest_calipso_cloud_mid_temperature', temp_2d, error_status
   tc_1d = reform(temp_2d[0,*])

   n = n_elements(num_cld_layers_1d)

   sc_1d = phase_1d
   sc_1d[*] = -1
   idx = where(phase_1d ne 2 or (tc_1d gt 0.0 and tc_1d gt 273.15),cc) & if (cc gt 0) then sc_1d[idx] = 0
   idx = where(phase_1d eq 2 and tc_1d le 273.15 and tc_1d gt 0.0,cc) & if (cc gt 0) then sc_1d[idx] = 1

   ;---- do some filtering on phase_1d
   idx = where(cod_1d lt cod_clear_thresh and cod_1d gt 0.0,cc) 
   if (cc gt 0) then begin
       sc_1d[idx] = 0
   endif

   for m = 0, n - 1 do begin

      if (phase_1d[m] lt 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 (sc_1d[m] eq 1) then begin  
         sc_cld_frac_sum_table[i,j,k,0] = sc_cld_frac_sum_table[i,j,k,0] + 1.0
         sc_cld_frac_sum_table[i,j,k,d] = sc_cld_frac_sum_table[i,j,k,d] + 1.0
      endif

      skip:

   endfor
cycle:
endfor


idx = where(count_table gt count_min,cc)
if (cc gt 0) then begin
   sc_cld_frac_table[idx] = sc_cld_frac_sum_table[idx] / count_table[idx]
endif

sc_cld_frac_table_smoothed = smooth(sc_cld_frac_table,[3,3,2,1],/Edge_Mirror,Missing=missing_output)

load_color_table, 33, 0, 1.0, 249

x0 = sc_cld_frac_table
idx = where(x0 eq missing_output, cc) & if (cc gt 0) then x0[idx] = 0.0

window, 0, xsize = 800, ysize = 600
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(mean(reform(x0[*,*,*,0]), dimension=3)), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Sc Cloud Fraction'
map_continents, color = 3, thick = 2, /lores

window, 1, xsize = 800, ysize = 600
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(mean(reform(x0[*,*,*,1]), dimension=3)), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Sc Cloud Fraction'
map_continents, color = 3, thick = 2, /lores

window, 2, xsize = 800, ysize = 600
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(mean(reform(x0[*,*,*,2]), dimension=3)), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Sc Cloud Fraction'
map_continents, color = 3, thick = 2, /lores

;---------------------------------
nc_filename = 'sc_prob_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, "timestamp",timestamp(), /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, 'sc_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, sc_cld_frac_table

;--- close netcdf files
ncdf_close, sd_id

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