function [data] = read_mod06_var(filenames, varToRead)

% Read one variable from MODIS Level2 HDF data files, and
% zoom Cloud Phase, Cloud Top Pressure and Effective Cloud Emissivity
% from 5km to 1km resolution
%
% Input
%   filenames -- Data file names
%   varToRead -- the variable to be read
%                Cloud_Phase
%                Cloud_Top_Pressure
%                Effective_Cloud_Emissivity
%                Effective_Particle_Radius
%                Cloud_Optical_Thickness
%                Water_Path
%                Cirrus_Reflectance
%                Cloud_Top_Temperature
%                Retrieval_Phase
%
% Output: data -- cloud product at 1km resolution
%                 zoom 5km to 1km:
%                   Cloud_Phase 
%                   Cloud_Top_Pressure
%                   Effective_Cloud_Emissivity
% 
% Merge IR Cloud Phase, 8 possible values
%     0 = clear  (from cloud mask)
%     1 = water
%     2 = ice
%     3 = mixed phase
%     4 = ice
%     5 = water
%     6 = uncertain
%   127 = missing data
% To 5 possible values:
%     1 -- clear
%     2 -- water cloud
%     3 -- ice cloud
%     4 -- mixed cloud
%     5 -- undecided
%
% Fill Value: -9999
%

disp(['Read variable: ' varToRead])

fill_value = -9999;

data = [];
for n = 1:length(filenames)
  tmp = modis_mod06_read(filenames{n}, varToRead);
  tmp_data = tmp.pimage;
  if strcmp(varToRead, 'Cloud_Phase') == 1 || ...
     strcmp(varToRead, 'Cloud_Top_Pressure') == 1 || ...
     strcmp(varToRead, 'Effective_Cloud_Emissivity') == 1
     t = modis_mod06_read(filenames{n}, 'Water_Path');
     dim = size(t.pimage);
     tmp_data = zoomImage(tmp_data, dim);
  end
  data = [data tmp_data];
end

if strcmp(varToRead, 'Cloud_Phase') == 1
  cphase = 5 * ones(size(data));
  ind_clear = find(data == 0);
  ind_water = find(data == 1 | data == 5);
  ind_ice = find(data == 2 | data == 4);
  ind_mixed = find(data == 3);
  cphase(ind_clear) = 1;
  cphase(ind_water) = 2;
  cphase(ind_ice) = 3;
  cphase(ind_mixed) = 4;
  data = cphase;
end

if strcmp(varToRead, 'Effective_Cloud_Emissivity') == 1
  ind = find(data > 1);
else
  ind = find(data < 0);
end
data(ind) = fill_value;
