Plot MODIS L2 Gridded Data (all cloudy pixels) with python 3

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from pyhdf.SD import SD, SDC 

import pprint
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm

file_name = 'gridded_data_2008_01.txt'

input_data = np.loadtxt(file_name, skiprows=1)

#----------------------------------------------------------------------------------------#
# plot cloud fraction

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,1] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

fig = plt.figure()

m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90, 
			llcrnrlon=-180,urcrnrlon=180,resolution='c')

m.drawcoastlines()
m.drawparallels(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))

plt.title(r'$P(Cloudy | 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - January 2008', size=8)

cmap = [(0.75,0.75,0.75)] + [(cm.jet (i)) for i in range(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [-1.0,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(data, cmap=cmap, norm=norm, interpolation='none')

bounds = [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]

cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, 
				ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])

cb.ax.tick_params(labelsize=8) 

plt.savefig('global_map_monthly_modis_l2_gridded_data_cloud_mask.png', bbox_inches='tight', dpi=200)
plt.close()

#----------------------------------------------------------------------------------------#

How to plot a confusion matrix with matplotlib and seaborn

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

import seaborn as sns
import math

from mpl_toolkits.axes_grid1 import make_axes_locatable

import matplotlib as mpl

mpl.style.use('seaborn')

conf_arr = np.array([[88,14,4],[12,85,11],[5,15,91]])

sum = conf_arr.sum()

conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

df_cm = pd.DataFrame(conf_arr, 
  index = [ 'Dog', 'Cat', 'Rabbit'],
  columns = ['Dog', 'Cat', 'Rabbit'])

fig = plt.figure()

plt.clf()

ax = fig.add_subplot(111)
ax.set_aspect(1)

cmap = sns.cubehelix_palette(light=1, as_cmap=True)

res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)

res.invert_yaxis()

plt.yticks([0.5,1.5,2.5], [ 'Dog', 'Cat', 'Rabbit'],va='center')

plt.title('Confusion Matrix')

plt.savefig('confusion_matrix.png', dpi=100, bbox_inches='tight' )

plt.close()

Tracer une matrice de confusion avec matplotlib et seaborn

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

import seaborn as sns
import math

from mpl_toolkits.axes_grid1 import make_axes_locatable

import matplotlib as mpl

mpl.style.use('seaborn')

conf_arr = np.array([[88,14,4],[12,85,11],[5,15,91]])

sum = conf_arr.sum()

conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

df_cm = pd.DataFrame(conf_arr, 
  index = [ 'Dog', 'Cat', 'Rabbit'],
  columns = ['Dog', 'Cat', 'Rabbit'])

fig = plt.figure()

plt.clf()

ax = fig.add_subplot(111)
ax.set_aspect(1)

cmap = sns.cubehelix_palette(light=1, as_cmap=True)

res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)

res.invert_yaxis()

plt.yticks([0.5,1.5,2.5], [ 'Dog', 'Cat', 'Rabbit'],va='center')

plt.title('Confusion Matrix')

plt.savefig('confusion_matrix.png', dpi=100, bbox_inches='tight' )

plt.close()

Plot cldclass-lidar granule vertical profile using python 3

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from pyhdf.HDF import *
from pyhdf.VS import *

import pprint
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import matplotlib.cm as cm

#----------------------------------------------------------------------------------------#

file_path = ''
file_name = '2008183012329_11573_CS_2B-CLDCLASS-LIDAR_GRANULE_P_R04_E02.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files (VD data) Latitude & Longitude

f = HDF(file_path+file_name, SDC.READ) 
vs = f.vstart() 

Latitude = vs.attach('Latitude')
Longitude = vs.attach('Longitude')

latitude_table = Latitude[:]

Latitude.detach() 
Longitude.detach() 
vs.end()
f.close() 

#print(latitude_table[0:100])

#----------------------------------------------------------------------------------------#
# SDS Data

file = SD(file_path+file_name, SDC.READ)

file_info = file.info()
print(file_info)  # number of sds and metadata

print('---------- CloudLayerBase ----------')

sds_obj = file.select('CloudLayerBase') # select sds

CloudLayerBase = sds_obj.get() 

sds_info = sds_obj.info()

print(CloudLayerBase.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

for i in range(10):
	print([CloudLayerBase[i,j] for j in range(10)])

print('---------- CloudLayerTop ----------')

sds_obj = file.select('CloudLayerTop') # select sds

CloudLayerTop = sds_obj.get() 

sds_info = sds_obj.info()

print(CloudLayerTop.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

for i in range(10):
	print([CloudLayerTop[i,j] for j in range(10)])
	
print('---------- CloudPhase ----------')

sds_obj = file.select('CloudPhase') # select sds

CloudPhase = sds_obj.get() 

sds_info = sds_obj.info()

print(CloudPhase.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

for i in range(10):
	print([CloudPhase[i,j] for j in range(10)])


cldclass_lidar_start_idx = 000
cldclass_lidar_end_idx = 1000

tag = '2'
	
#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='cloud')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='lightblue', edgecolor='lightblue', label='clear')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)

plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Cloud Mask')

plt.savefig('plot_'+tag+'_cldclass_lidar_cloud_layers_method_01.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers + nb cloud layers (aggregation 1)

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			if nb_cloud_layer == 2:
				bar_color = 'blue'
				bar_edgecolor = 'blue'
			if nb_cloud_layer > 2:
				bar_color = 'salmon'
				bar_edgecolor = 'salmon'
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='nb cloud layer = 1')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='blue', edgecolor='blue', label='nb cloud layer = 2')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='salmon', edgecolor='salmon', label='nb cloud layer > 2')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)

plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Number of cloud layers (Method 1)')

plt.savefig('plot_'+tag+'_cldclass_lidar_cloud_layers_and_nb_cloud_layers_method_01.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers + nb cloud layers (aggregation 2: merge two layers if
# same cloud phase and separation distance < threshold (3 for instance) )

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	#for j in range(nb_cloud_layer):
	#	print(CloudLayerBase[i,j],CloudLayerTop[i,j])
	#print('-------')
	nb_cloud_layer_after_merging = nb_cloud_layer
	for j in range(nb_cloud_layer-1):	
			if ( CloudLayerBase[i,j+1] - CloudLayerTop[i,j] ) < 3.0:
				if CloudPhase[i,j+1] == CloudPhase[i,j]:
					nb_cloud_layer_after_merging = nb_cloud_layer_after_merging - 1
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			if nb_cloud_layer_after_merging == 2:
				bar_color = 'blue'
				bar_edgecolor = 'blue'
			if nb_cloud_layer_after_merging > 2:
				bar_color = 'salmon'
				bar_edgecolor = 'salmon'
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='nb cloud layer = 1')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='blue', edgecolor='blue', label='nb cloud layer = 2')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='salmon', edgecolor='salmon', label='nb cloud layer > 2')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)

plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Number of cloud layers (Method 2)')

plt.savefig('plot_'+tag+'_cldclass_lidar_cloud_layers_and_nb_cloud_layers_aggregated.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers + nb cloud phase

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			if CloudPhase[i,j] == 3:
				bar_color = 'blue'
				bar_edgecolor = 'blue'
			if CloudPhase[i,j] == 2:
				bar_color = 'salmon'
				bar_edgecolor = 'salmon'			
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='ice clouds')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='salmon', edgecolor='salmon', label='mixed clouds')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='blue', edgecolor='blue', label='liquid clouds')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)

plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Cloud thermodynamic phase')

plt.savefig( 'plot_'+tag+'_cldclass_lidar_cloud_layers_and_thermodynamic_phase_method_01.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

MODIS multilayer clouds vs CoudSat & CALIPSO (analysis main python code)

					
#!/usr/bin/env python

import numpy as np

import sample_homogenization_module 

import modis_cldclass_lidar_analysis_001 as analysis_001
import modis_cldclass_lidar_analysis_002 as analysis_002
import modis_cldclass_lidar_analysis_003 as analysis_003

#----------------------------------------------------------------------------------------#
# Settings

dict_settings = {}

dict_settings['input_path'] = '/Users/bmarcha1/Desktop/media/files/monthly_modis_cldclass_lidar/'

dict_settings['output_path'] = '../outputs/'
dict_settings['output_file_name_prefix'] = 'modis_cldclass_lidar_analyisis_'

month_list = [1,2,3,4,5,6,7,8,9,10,11,12]

for idx,month in enumerate(month_list):
	if idx > 0:
		data_1 = np.loadtxt(
		dict_settings['input_path'] + '2008_{:02d}_sample_01_modis_caliop_cldclass_lidar.txt'.format(month), skiprows=1)
		data = np.concatenate((data,data_1), axis=0)
	else:		
		data = np.loadtxt(
		dict_settings['input_path'] + '2008_{:02d}_sample_01_modis_caliop_cldclass_lidar.txt'.format(month), skiprows=1)

#----------------------------------------------------------------------------------------#
# sample homogeneization

for idx,month in enumerate(month_list):
	if idx > 0:
		data_1 = np.loadtxt(
		dict_settings['input_path'] + '2008_{:02d}_sample_00_modis_caliop_cldclass_lidar.txt'.format(month), skiprows=1)
		data_from_sample_0 = np.concatenate((data_from_sample_0,data_1), axis=0)
	else:		
		data_from_sample_0 = np.loadtxt(
		dict_settings['input_path'] + '2008_{:02d}_sample_00_modis_caliop_cldclass_lidar.txt'.format(month), skiprows=1)

prob_cloudy_bin = sample_homogenization_module.calculate_probabilities_from_sample_0(data_from_sample_0)

data = sample_homogenization_module.homogenization(data,prob_cloudy_bin)

#----------------------------------------------------------------------------------------#
# confusion matrices (multilayer clouds)

fig_count = 0

dict_settings['fig_title'] = 'modis c6 vs cldclass-lidar'
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 (without PH test) vs cldclass-lidar'
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['product_source'] = 3

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 vs caliop 5km'
dict_settings['fig_ylabel'] = 'caliop 5km'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 2

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 vs caliop 1km'
dict_settings['fig_ylabel'] = 'caliop 1km'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 1

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# confusion matrices (cloud phase)

dict_settings['fig_title'] = 'Cloud Thermodynamic Phase Comparisons \n modis c6 vs cldclass-lidar (nb cloud layer = 1) \n score = '
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = 0.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_001.cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# cldclass-lidar (nb cloud layer) -- Cloud Layer Vertical Profile Histogram

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60

fig_count = analysis_001.cldclass_lidar_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# cldclass-lidar (cloud phase) -- Cloud Layer Vertical Profile Histogram

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) > 0'

fig_count = analysis_001.cldclass_lidar_cloud_phase_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution \n (nb cloud layer = 1)'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 1'

fig_count = analysis_001.cldclass_lidar_cloud_phase_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution \n (nb cloud layer = 2)'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 2'

fig_count = analysis_001.cldclass_lidar_cloud_phase_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution \n (nb cloud layer > 1)'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) > 1'

fig_count = analysis_001.cldclass_lidar_cloud_phase_layer_vertical_distribution_bar(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# cldclass-lidar vs modis (nb cloud layer) -- Cloud Layer Vertical Profile Histogram

dict_settings['fig_title'] = 'modis (monolayer cloud) vs cldclass-lidar \n cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_multilayer_flag'] = 0
dict_settings['modis_without_pavolonis_heidinger_test'] = False

fig_count = analysis_001.cldclass_lidar_vs_modis_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis (multilayer cloud) vs cldclass-lidar \n cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_multilayer_flag'] = 1
dict_settings['modis_without_pavolonis_heidinger_test'] = False

fig_count = analysis_001.cldclass_lidar_vs_modis_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis (Undet. multilayer cloud flag (COT < 4)) vs cldclass-lidar \n cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_multilayer_flag'] = 2
dict_settings['modis_without_pavolonis_heidinger_test'] = False

fig_count = analysis_001.cldclass_lidar_vs_modis_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
#  cldclass-lidar vs modis (cloud phase) -- Cloud Layer Vertical Profile Histogram

dict_settings['fig_title'] = 'modis (ice phase) vs cldclass-lidar (nb cloud layer = 1) \n cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_cloud_phase_flag'] = 3
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 1'

fig_count = analysis_001.modis_vs_cldclass_lidar_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis (liquid phase) vs cldclass-lidar (nb cloud layer = 1) \n cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_cloud_phase_flag'] = 2
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 1'

fig_count = analysis_001.modis_vs_cldclass_lidar_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis (undet. phase) vs cldclass-lidar (nb cloud layer = 1) \n cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_cloud_phase_flag'] = 4
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 1'

fig_count = analysis_001.modis_vs_cldclass_lidar_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# matrices re diff

dict_settings['fig_title'] = 'modis c6 (liquid) cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 (ice) cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# modis cloud optical thickness vs cldclass-liar histograms

dict_settings['fig_title'] = 'modis c6 (liq) cloud optical thickness vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 cloud optical thickness'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_optical_thickness_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 (ice) cloud optical thickness vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 cloud optical thickness'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_optical_thickness_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# confusion matrices (homogeneous cloud phase, n = 2)

dict_settings['fig_title'] = 'Cloud Thermodynamic Phase Comparisons \n modis c6 vs cldclass-lidar (nb cloud layer = 2)'
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = 0.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_001.homogeneous_cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# confusion matrices (inhomogeneous cloud phase, n = 2)

dict_settings['fig_title'] = 'Cloud Thermodynamic Phase Comparisons \n modis c6 vs cldclass-lidar (nb cloud layer = 2)'
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = 0.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_001.inhomogeneous_cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# matrices re diff (cldclass-lidar 2 layers; inhomogeneous cloud phase)

dict_settings['fig_title'] = 'modis c6 (liquid) cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['fig_yticks'] = ['mixed / liquid', 'ice / mixed', 'ice / liquid']
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['product_source'] = 3
dict_settings['cldclass_lidar_profile'] = ((1,3),(1,2),(2,3))

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_heterogeneous_phase_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 (ice) cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['fig_yticks'] = ['mixed / liquid', 'ice / mixed', 'ice / liquid']
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['product_source'] = 3
dict_settings['cldclass_lidar_profile'] = ((1,3),(1,2),(2,3))

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_heterogeneous_phase_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# matrices re diff (cldclass-lidar 2 layers; homogeneous cloud phase)

dict_settings['fig_title'] = 'modis c6 (liquid) cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['fig_yticks'] = ['mixed / mixed', 'liquid / liquid', 'ice / ice']
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['product_source'] = 3
dict_settings['cldclass_lidar_profile'] = ((1,1),(3,3),(2,2))

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_heterogeneous_phase_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 (ice) cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['fig_yticks'] = ['mixed / mixed', 'liquid / liquid', 'ice / ice']
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['product_source'] = 3
dict_settings['cldclass_lidar_profile'] = ((1,1),(3,3),(2,2))

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_heterogeneous_phase_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# cldclass-lidar (nb cloud layer) -- Cloud Layer Vertical Profile Histogram

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60

fig_count = analysis_001.cldclass_lidar_and_caliop_5km_cot_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# Modis versus cldclass-lidar frequency table

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60

fig_count = analysis_002.modis_multilayer_cloud_vs_cldclass_lidar(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# zonal plots

dict_settings['output_file_name_prefix'] = 'modis_vs_cldclass_lidar_multilayer_cloud_zonal_histograms_'

dict_settings['fig_title'] = r'Zonal Histograms -  year 2008'
dict_settings['fig_y_label'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'Modis C6 multilayer clouds'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar multilayer clouds'
dict_settings['nb_layers_flag'] = True
dict_settings['distance_threshold'] = 0
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0

fig_count = analysis_001.zonal_plots(data,dict_settings,fig_count)

dict_settings['fig_title'] = r'Zonal Histograms ($\tau > 4.0$) -  year 2008'
dict_settings['fig_y_label'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'Modis C6 multilayer clouds'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar multilayer clouds'
dict_settings['nb_layers_flag'] = True
dict_settings['distance_threshold'] = 0
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0

fig_count = analysis_001.zonal_plots(data,dict_settings,fig_count)

dict_settings['fig_title'] = r'Zonal Histograms ($\tau > 4.0$) -  year 2008'
dict_settings['fig_y_label'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'Modis C6 multilayer clouds \n (without PH test)'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar multilayer clouds'
dict_settings['nb_layers_flag'] = True
dict_settings['distance_threshold'] = 0
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0

fig_count = analysis_001.zonal_plots(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# modis cloud effective radius histograms vs cldclass-lidar monolayer cloud

dict_settings['output_file_name_prefix'] = 'modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud_'

dict_settings['fig_title'] = ''
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($2.1\mu m$)'
dict_settings['fig_ylabel'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'cldclass-lidar: liquid monolayer cloud \n modis: liquid cloud'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar: ice monolayer cloud \n modis: liquid cloud'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count)

dict_settings['fig_title'] = ''
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($1.6\mu m$)'
dict_settings['fig_ylabel'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'cldclass-lidar: liquid monolayer cloud \n modis: liquid cloud'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar: ice monolayer cloud \n modis: liquid cloud'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count)

dict_settings['fig_title'] = ''
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['fig_ylabel'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'cldclass-lidar: liquid monolayer cloud \n modis: liquid cloud'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar: ice monolayer cloud \n modis: liquid cloud'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count)

dict_settings['fig_title'] = ''
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($2.1\mu m$)'
dict_settings['fig_ylabel'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'cldclass-lidar: ice monolayer cloud \n modis: ice cloud'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar: liquid monolayer cloud \n modis: ice cloud'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count)

dict_settings['fig_title'] = ''
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($1.6\mu m$)'
dict_settings['fig_ylabel'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'cldclass-lidar: ice monolayer cloud \n modis: ice cloud'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar: liquid monolayer cloud \n modis: ice cloud'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count)

dict_settings['fig_title'] = ''
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['fig_ylabel'] = 'Probability Density Function'
dict_settings['fig_bar_1_label'] = 'cldclass-lidar ice monolayer cloud \n modis: ice cloud'
dict_settings['fig_bar_2_label'] = 'cldclass-lidar: liquid monolayer cloud \n modis: ice cloud'
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# modis cloud effective radius vs cldclass-lidar histograms

dict_settings['output_file_name_prefix'] = 'modis_cloud_effective_radius_hist_vs_cldclass_lidar_'

dict_settings['fig_title'] = 'modis c6 cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 cloud effective radius vs cldclass-lidar'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# modis (cot < 4) cloud effective radius histograms vs cldclass-liar 

dict_settings['output_file_name_prefix'] = 'modis_cot_lt_4_cloud_effective_radius_hist_vs_cldclass_lidar_'

dict_settings['fig_title'] = 'Modis C6: Cloud Optical Thickness < 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0 and modis_cloud_optical_tickness < 4.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Cloud Optical Thickness < 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0 and modis_cloud_optical_tickness < 4.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Cloud Optical Thickness < 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0 and modis_cloud_optical_tickness < 4.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Cloud Optical Thickness < 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0 and modis_cloud_optical_tickness < 4.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Cloud Optical Thickness < 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0 and modis_cloud_optical_tickness < 4.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Cloud Optical Thickness < 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 0.0 and modis_cloud_optical_tickness < 4.0'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# modis (cot > 4) monolayer cloud effective radius histograms vs cldclass-liar 

dict_settings['output_file_name_prefix'] = 'modis_cot_gt_4_monolayer_cloud_effective_radius_hist_vs_cldclass_lidar_'

dict_settings['fig_title'] = 'Modis C6: Monolayer & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# modis (cot > 4) monolayer (PH turned off) cloud effective radius histograms vs cldclass-liar 

dict_settings['output_file_name_prefix'] = 'modis_cot_gt_4_monolayer_no_ph_cloud_effective_radius_hist_vs_cldclass_lidar_'

dict_settings['fig_title'] = 'Modis C6: Monolayer (whitout PH) & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer (whitout PH) & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer (whitout PH) & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 liquid cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_phase_label'] = 'liquid'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 2'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer (whitout PH) & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($2.1\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 10
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer (whitout PH) & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($1.6\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 11
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'Modis C6: Monolayer (whitout PH) & Cloud Optical Thickness > 4'
dict_settings['fig_xlabel'] = r'modis c6 ice cloud effective radius ($3.7\mu m$)'
dict_settings['modis_cloud_optical_thickness_expression'] = 'modis_cloud_optical_tickness > 4.0 and modis_multilayer_sds == 1'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['modis_cloud_phase_label'] = 'ice'
dict_settings['modis_cloud_phase_logical_expression'] = 'int(data[i,6]) == 3'
dict_settings['modis_cloud_effective_radius_idx'] = 12
dict_settings['product_source'] = 3

fig_count = analysis_001.modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# probability of modis ice phase given cldclass-lidar upeer layer cloud optical thickness and d

dict_settings['output_file_name_prefix'] = 'probability_of_modis_cloud_phase_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = ''
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_003.modis_ice_phase_probability_given_utau_and_d(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# probability of modis multilayer given cldclass-lidar upeer layer cloud optical thickness and d

dict_settings['output_file_name_prefix'] = 'probability_of_modis_multilayer_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = ''
dict_settings['modis_cloud_optical_thickness_threshold'] = -999.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['cldclass_lidar_merge_cloud_layers'] = False
dict_settings['cldclass_lidar_merge_cloud_layers_threshold'] = 0.0
dict_settings['product_source'] = 3

fig_count = analysis_003.distance_vs_upper_layer_cloud_optical_thickness(data,dict_settings,fig_count)

dict_settings['output_file_name_prefix'] = 'probability_of_modis_multilayer_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = '\n modis cloud optical thickness > 4.0 '
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['cldclass_lidar_merge_cloud_layers'] = False
dict_settings['cldclass_lidar_merge_cloud_layers_threshold'] = 0.0
dict_settings['product_source'] = 3

fig_count = analysis_003.distance_vs_upper_layer_cloud_optical_thickness(data,dict_settings,fig_count)

dict_settings['output_file_name_prefix'] = 'probability_of_modis_multilayer_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = '\n modis cloud optical thickness > 4.0 and without PH test '
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['cldclass_lidar_merge_cloud_layers'] = False
dict_settings['cldclass_lidar_merge_cloud_layers_threshold'] = 0.0
dict_settings['product_source'] = 3

fig_count = analysis_003.distance_vs_upper_layer_cloud_optical_thickness(data,dict_settings,fig_count)

dict_settings['output_file_name_prefix'] = 'probability_of_modis_multilayer_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = '\n modis cloud optical thickness > 4.0 and without PH test '
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['cldclass_lidar_merge_cloud_layers'] = True
dict_settings['cldclass_lidar_merge_cloud_layers_threshold'] = 3.0
dict_settings['product_source'] = 3

fig_count = analysis_003.distance_vs_upper_layer_cloud_optical_thickness(data,dict_settings,fig_count)

dict_settings['output_file_name_prefix'] = 'probability_of_modis_multilayer_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = '\n modis cloud optical thickness > 4.0'
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = False
dict_settings['cldclass_lidar_merge_cloud_layers'] = True
dict_settings['cldclass_lidar_merge_cloud_layers_threshold'] = 3.0
dict_settings['product_source'] = 3

fig_count = analysis_003.distance_vs_upper_layer_cloud_optical_thickness(data,dict_settings,fig_count)


dict_settings['output_file_name_prefix'] = 'probability_of_modis_multilayer_given_cldclass_lidar_utau_and_d_'
dict_settings['fig_title'] = '\n modis cloud optical thickness > 4.0 and without PH test '
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['cldclass_lidar_merge_cloud_layers'] = True
dict_settings['cldclass_lidar_merge_cloud_layers_threshold'] = 3.0
dict_settings['product_source'] = 3
dict_settings['cldclass_lidar_cloud_phase_flag'] = True
dict_settings['cldclass_lidar_cloud_phase_expr'] = 'cldclass_lidar_cloud_layer_phase_list[0] == 1 and cldclass_lidar_cloud_layer_phase_list[1] == 3' # ice overlapping liquid cloud

fig_count = analysis_003.distance_vs_upper_layer_cloud_optical_thickness(data,dict_settings,fig_count)

MODIS multilayer clouds vs CoudSat & CALIPSO (analysis 01 python code)

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

import seaborn as sns
import math

from mpl_toolkits.axes_grid1 import make_axes_locatable

import matplotlib as mpl

mpl.style.use('seaborn')

#----------------------------------------------------------------------------------------#
# Confusion Matrices (Multilayer Clouds)

def multilayer_cloud_confusion_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	conf_arr = np.zeros((3,2))

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
		
			if 'modis_without_pavolonis_heidinger_test' in dict_settings:
				if dict_settings['modis_without_pavolonis_heidinger_test']:
					if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
		
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
		
			modis_multilayer_flag = -1
			if modis_multilayer_sds == 1 : modis_multilayer_flag = 0 
			if modis_multilayer_sds > 1: modis_multilayer_flag = 1
			if modis_cloud_optical_tickness < 4.0: modis_multilayer_flag = 2
		
			cldclass_lidar_multilayer_flag = -1
			if cldclass_lidar_nb_cloud_layer == 1: cldclass_lidar_multilayer_flag = 0
			if cldclass_lidar_nb_cloud_layer > 1: cldclass_lidar_multilayer_flag = 1		
		
			if 	cldclass_lidar_multilayer_flag >= 0 and modis_multilayer_flag >= 0:
				conf_arr[modis_multilayer_flag,cldclass_lidar_multilayer_flag] = conf_arr[modis_multilayer_flag,cldclass_lidar_multilayer_flag] + 1
		
	sum = conf_arr.sum()
	
	score = ( conf_arr[0,0] + conf_arr[1,1] ) * 100.0 / ( conf_arr[0,0] + conf_arr[1,0] + conf_arr[0,1] + conf_arr[1,1] )
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = [ 'modis C6 Undet.', 'modis C6 monolayer cloud', 'modis C6 multilayer cloud'],
					  columns = [dict_settings['fig_ylabel']+'\n nb cloud layer = 1', dict_settings['fig_ylabel']+'\n nb cloud layer > 1'])

	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)

	cmap = sns.cubehelix_palette(light=1, as_cmap=True)

	res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], ['modis C6 \n Undet. (COT < 4.0) ', 'modis C6 \n multilayer cloud', 'modis C6 \n monolayer cloud'],va='center')

	plt.title(dict_settings['fig_title'])

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# Confusion Matrices (Cloud Thermodynamic Phase; cldclass-lidar monolayer)

def cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	conf_arr = np.zeros((3,3))

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]
		modis_cloud_phase_sds = data[i,6]		
		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
			
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					#print( data[i,cloud_layer_idx+4], data[i,cloud_layer_idx], data[i,cloud_layer_idx+1] )
					if int(data[i,cloud_layer_idx+4]) == 3: # source cldclass-lidar products
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5
			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
			
			#print(cldclass_lidar_phase_list,cldclass_lidar_nb_cloud_layer)
			#print('-------')
			
			#cldclass_lidar_multiple_phase = False
			#cloud_phase_layer_idx = 22
			#for j in range(int(cldclass_lidar_nb_cloud_layer)):
			#	if int(data[i,cloud_phase_layer_idx]) != int(cldclass_lidar_phase_sds):
			#		cldclass_lidar_multiple_phase = True
			#	cloud_phase_layer_idx = cloud_phase_layer_idx + 5
						
			if cldclass_lidar_nb_cloud_layer == 1:
			#if not cldclass_lidar_multiple_phase:
				
				cldclass_lidar_phase_sds = cldclass_lidar_phase_list[0]
				
				modis_cloud_phase_flag = -1
				if modis_cloud_phase_sds == 3: modis_cloud_phase_flag = 0 # ice
				if modis_cloud_phase_sds == 2: modis_cloud_phase_flag = 1 # liquid
				if modis_cloud_phase_sds == 4: modis_cloud_phase_flag = 2 # undetermined
		
				cldclass_lidar_cloud_phase_flag = -1
				if cldclass_lidar_phase_sds == 1: cldclass_lidar_cloud_phase_flag = 0 # ice
				if cldclass_lidar_phase_sds == 3: cldclass_lidar_cloud_phase_flag = 1 # liquid
				if cldclass_lidar_phase_sds == 2: cldclass_lidar_cloud_phase_flag = 2 # mixed	
		
				if 	modis_cloud_phase_flag >= 0 and cldclass_lidar_cloud_phase_flag >= 0:
					conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] = conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] + 1
		
	sum = conf_arr.sum()
	
	score = ( conf_arr[0,0] + conf_arr[1,1] ) * 100.0 / (conf_arr[0,0] + conf_arr[1,0] + conf_arr[0,1] + conf_arr[1,1])
		
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = ['modis C6 ice phase', 'modis C6 liquid phase', 'modis C6 undet. phase'],
					  columns = [dict_settings['fig_ylabel']+'\n ice phase', dict_settings['fig_ylabel']+'\n liquid phase', dict_settings['fig_ylabel']+'\n mixed phase'])

	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)

	cmap = sns.cubehelix_palette(light=1, as_cmap=True)

	res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], ['modis C6 \n undet. phase', 'modis C6 \n liquid phase', 'modis C6 \n ice phase'],va='center')

	plt.title(dict_settings['fig_title'] + str(round(score,2)))

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# Confusion Matrices (Modis cloud phase; cldclass-lidar two layers with homogeneous phase, nb -layers = 2)

def homogeneous_cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	conf_arr = np.zeros((3,3))
	
	nb_cloud_layer_2 = 0
	nb_homgeneous_cloud_phase = 0
	
	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]
		modis_cloud_phase_sds = data[i,6]		
		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
			
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					#print( data[i,cloud_layer_idx+4], data[i,cloud_layer_idx], data[i,cloud_layer_idx+1] )
					if int(data[i,cloud_layer_idx+4]) == 3: # source cldclass-lidar products
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5
			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
								
			if cldclass_lidar_nb_cloud_layer == 2:
				
				nb_cloud_layer_2 = nb_cloud_layer_2 + 1
				
				cldclass_lidar_phase_upper_layer = cldclass_lidar_phase_list[0]
				cldclass_lidar_phase_lower_layer = cldclass_lidar_phase_list[1]
				
				modis_cloud_phase_flag = -1
				if modis_cloud_phase_sds == 3: modis_cloud_phase_flag = 0 # ice
				if modis_cloud_phase_sds == 2: modis_cloud_phase_flag = 1 # liquid
				if modis_cloud_phase_sds == 4: modis_cloud_phase_flag = 2 # undetermined
		
				cldclass_lidar_cloud_phase_flag = -1
				if cldclass_lidar_phase_upper_layer == 1 and cldclass_lidar_phase_lower_layer == 1: cldclass_lidar_cloud_phase_flag = 0
				if cldclass_lidar_phase_upper_layer == 3 and cldclass_lidar_phase_lower_layer == 3: cldclass_lidar_cloud_phase_flag = 1
				if cldclass_lidar_phase_upper_layer == 2 and cldclass_lidar_phase_lower_layer == 2: cldclass_lidar_cloud_phase_flag = 2
				
				if cldclass_lidar_cloud_phase_flag >= 0:
					nb_homgeneous_cloud_phase = nb_homgeneous_cloud_phase + 1
						
				if 	modis_cloud_phase_flag >= 0 and cldclass_lidar_cloud_phase_flag >= 0:
					conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] = conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] + 1
	
	print('Homogeneous cloud phase, nb cloud layer = 2', nb_homgeneous_cloud_phase * 100.0 / nb_cloud_layer_2)
	
	sum = conf_arr.sum()
			
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = ['modis C6 ice phase', 'modis C6 liquid phase', 'modis C6 undet. phase'],
					  columns = [dict_settings['fig_ylabel']+'\n ice / ice ', dict_settings['fig_ylabel']+'\n liquid / liquid', dict_settings['fig_ylabel']+'\n mixed / mixed'])
	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)
	
	cmap = sns.cubehelix_palette(light=1, as_cmap=True)

	res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], ['modis C6 \n undet. phase', 'modis C6 \n liquid phase', 'modis C6 \n ice phase'],va='center')
	plt.xticks(rotation='vertical')

	plt.title(dict_settings['fig_title'])

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# Confusion Matrices (Modis cloud phase; cldclass-lidar two layers with inhomogeneous phase, nb -layers = 2)

def inhomogeneous_cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	conf_arr = np.zeros((3,6))

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]
		modis_cloud_phase_sds = data[i,6]		
		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
			
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					#print( data[i,cloud_layer_idx+4], data[i,cloud_layer_idx], data[i,cloud_layer_idx+1] )
					if int(data[i,cloud_layer_idx+4]) == 3: # source cldclass-lidar products
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5
			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
								
			if cldclass_lidar_nb_cloud_layer == 2:
				
				cldclass_lidar_phase_upper_layer = cldclass_lidar_phase_list[0]
				cldclass_lidar_phase_lower_layer = cldclass_lidar_phase_list[1]
				
				modis_cloud_phase_flag = -1
				if modis_cloud_phase_sds == 3: modis_cloud_phase_flag = 0 # ice
				if modis_cloud_phase_sds == 2: modis_cloud_phase_flag = 1 # liquid
				if modis_cloud_phase_sds == 4: modis_cloud_phase_flag = 2 # undetermined
		
				cldclass_lidar_cloud_phase_flag = -1
				if cldclass_lidar_phase_upper_layer == 1 and cldclass_lidar_phase_lower_layer == 3: cldclass_lidar_cloud_phase_flag = 0
				if cldclass_lidar_phase_upper_layer == 1 and cldclass_lidar_phase_lower_layer == 2: cldclass_lidar_cloud_phase_flag = 1
				if cldclass_lidar_phase_upper_layer == 3 and cldclass_lidar_phase_lower_layer == 1: cldclass_lidar_cloud_phase_flag = 2
				if cldclass_lidar_phase_upper_layer == 3 and cldclass_lidar_phase_lower_layer == 2: cldclass_lidar_cloud_phase_flag = 3
				if cldclass_lidar_phase_upper_layer == 2 and cldclass_lidar_phase_lower_layer == 1: cldclass_lidar_cloud_phase_flag = 4
				if cldclass_lidar_phase_upper_layer == 2 and cldclass_lidar_phase_lower_layer == 3: cldclass_lidar_cloud_phase_flag = 5
						
				if 	modis_cloud_phase_flag >= 0 and cldclass_lidar_cloud_phase_flag >= 0:
					conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] = conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] + 1
		
	sum = conf_arr.sum()
			
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = ['modis C6 ice phase', 'modis C6 liquid phase', 'modis C6 undet. phase'],
					  columns = [dict_settings['fig_ylabel']+'\n ice / liquid ', dict_settings['fig_ylabel']+'\n ice / mixed', dict_settings['fig_ylabel']+'\n liquid / ice',dict_settings['fig_ylabel']+'\n liquid / mixed', dict_settings['fig_ylabel']+'\n mixed / ice', dict_settings['fig_ylabel']+'\n mixed / liquid'])
	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)
	
	cmap = sns.cubehelix_palette(light=1, as_cmap=True)

	res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap, cbar_kws={"shrink": 0.5})
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], ['modis C6 \n undet. phase', 'modis C6 \n liquid phase', 'modis C6 \n ice phase'],va='center')
	plt.xticks(rotation='vertical')

	plt.title(dict_settings['fig_title'])

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# cldclass-lidar (nb cloud layer) Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def cldclass_lidar_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]
	
	width = 0.35  
	
	cloud_layer_count_bar_x = [i*0.1 + 0.05 for i in range(200)]
	cloud_layer_count_bar_height = [0 for i in range(200)]

	cloud_layer_mono_count_bar_height = [0 for i in range(200)]
	cloud_layer_multi_count_bar_height = [0 for i in range(200)]
	cloud_layer_two_count_bar_height = [0 for i in range(200)]	
		
	for i in range(data.shape[0]):
	#for i in range(100):
		
		nb_cloud_layer = int( data[i,15] )
		
		cloud_layer_parameters = np.zeros((10,5))
		
		m = 19
		if nb_cloud_layer <= 10:
			for j in range(nb_cloud_layer):
				cloud_layer_parameters[j,0] = data[i,m]
				cloud_layer_parameters[j,1] = data[i,m+1]
				cloud_layer_parameters[j,2] = data[i,m+2]
				cloud_layer_parameters[j,3] = data[i,m+3]
				cloud_layer_parameters[j,4] = data[i,m+4]
				m = m + 5
			
			#print('nb_cloud_layer --->', nb_cloud_layer)
			#for j in range(nb_cloud_layer):
			#	print( [ "{0:0.1f}".format(value) for value in list(cloud_layer_parameters[j,:]) ] )
			#print('------')
					
			for j in range(nb_cloud_layer):
				if cloud_layer_parameters[j,0] > 0.0 and cloud_layer_parameters[j,1] > 0.0:
					cloud_layer_top_idx = int(cloud_layer_parameters[j,0]/0.1)
					cloud_layer_base_idx = int(cloud_layer_parameters[j,1]/0.1)
					if int(cloud_layer_parameters[j,4]) == 3: # cldclass-lidar products
						for m in range(cloud_layer_base_idx,cloud_layer_top_idx+1):
							cloud_layer_count_bar_height[m] = cloud_layer_count_bar_height[m] + 1
							if int(data[i,18]) == 1: cloud_layer_mono_count_bar_height[m] = cloud_layer_mono_count_bar_height[m] + 1 
							if int(data[i,18]) > 1 : cloud_layer_multi_count_bar_height[m] = cloud_layer_multi_count_bar_height[m] + 1 
							if int(data[i,18]) == 2: cloud_layer_two_count_bar_height[m] = cloud_layer_two_count_bar_height[m] + 1 	
	
	plt.barh(cloud_layer_count_bar_x, cloud_layer_count_bar_height)
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	plt.barh(cloud_layer_count_bar_x, cloud_layer_multi_count_bar_height, width, color="lightcoral", alpha=1.0, align='center', label='nb cloud layer > 1')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_mono_count_bar_height, width, color="royalblue", alpha=1.0, align='center', label='nb cloud layer = 1')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_two_count_bar_height, width, color="wheat", alpha=0.8, align='center', label='nb cloud layer = 2')

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_two_count_bar_height[j],cloud_layer_two_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="wheat")

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_multi_count_bar_height[j],cloud_layer_multi_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="lightcoral")
		
	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_mono_count_bar_height[j],cloud_layer_mono_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="royalblue")				

	plt.legend()
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# cldclass-lidar (nb cloud layer) + caliop 5km cot: Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def cldclass_lidar_and_caliop_5km_cot_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]
	
	width = 0.35  
	
	cloud_layer_count_bar_x = [i*0.1 + 0.05 for i in range(200)]

	cldclass_lidar_two_cloud_layer_bar_height = [0 for i in range(200)]	
	cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height = [0 for i in range(200)]		
	
	caliop_5km_cloud_optical_thickness_hist = []
	
	for i in range(data.shape[0]):
	#for i in range(100):
		
		nb_cloud_layer = int( data[i,15] )
		
		caliop_5km_cloud_layer_top_list = []
		caliop_5km_cloud_layer_base_list = []
		caliop_5km_cloud_layer_cot_list = []
		caliop_5km_cloud_layer_phase_list = []

		cldclass_lidar_cloud_layer_top_list = []
		cldclass_lidar_cloud_layer_base_list = []
		cldclass_lidar_cloud_layer_phase_list = []
		
		cloud_layer_idx = 19
		for n in range(10):
			if int(data[i,cloud_layer_idx+4]) == 2: # caliop 5km
				caliop_5km_cloud_layer_top_list.append(data[i,cloud_layer_idx])
				caliop_5km_cloud_layer_base_list.append(data[i,cloud_layer_idx+1])
				caliop_5km_cloud_layer_cot_list.append(data[i,cloud_layer_idx+2])												
				caliop_5km_cloud_layer_phase_list.append(int(data[i,cloud_layer_idx+3]))
			if int(data[i,cloud_layer_idx+4]) == 3: # cldclass-lidar
				cldclass_lidar_cloud_layer_top_list.append(data[i,cloud_layer_idx])
				cldclass_lidar_cloud_layer_base_list.append(data[i,cloud_layer_idx+1])											
				cldclass_lidar_cloud_layer_phase_list.append(int(data[i,cloud_layer_idx+3]))
			cloud_layer_idx = cloud_layer_idx + 5			

		caliop_5km_nb_cloud_layer = len(caliop_5km_cloud_layer_phase_list)
		cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_cloud_layer_phase_list)
				
		if cldclass_lidar_nb_cloud_layer == 2:

			if cldclass_lidar_cloud_layer_top_list[0] > 0.0 and cldclass_lidar_cloud_layer_base_list[0] > 0.0:
				if cldclass_lidar_cloud_layer_top_list[1] > 0.0 and cldclass_lidar_cloud_layer_base_list[1] > 0.0: 

					cloud_upper_layer_top_idx = int( cldclass_lidar_cloud_layer_top_list[0]/0.1)
					cloud_upper_layer_base_idx = int( cldclass_lidar_cloud_layer_base_list[0]/0.1)

					cloud_lower_layer_top_idx = int( cldclass_lidar_cloud_layer_top_list[1]/0.1 )
					cloud_lower_layer_base_idx = int( cldclass_lidar_cloud_layer_base_list[1]/0.1 )				

					for m in range(cloud_upper_layer_base_idx,cloud_upper_layer_top_idx+1):
						cldclass_lidar_two_cloud_layer_bar_height[m] = cldclass_lidar_two_cloud_layer_bar_height[m] + 1 

					for m in range(cloud_lower_layer_base_idx,cloud_lower_layer_top_idx+1):
						cldclass_lidar_two_cloud_layer_bar_height[m] = cldclass_lidar_two_cloud_layer_bar_height[m] + 1 
					
					if caliop_5km_nb_cloud_layer > 0.0:
						if caliop_5km_cloud_layer_cot_list[0] > 0.0:			
							
							d2 = ( caliop_5km_cloud_layer_top_list[0] - cldclass_lidar_cloud_layer_top_list[0] )**2 + ( caliop_5km_cloud_layer_base_list[0] - cldclass_lidar_cloud_layer_base_list[0] )**2
							d = math.sqrt(d2)
							
							print(d)
							
							if d < 1.0:
								
								caliop_5km_cloud_optical_thickness_hist.append(caliop_5km_cloud_layer_cot_list[0])
								
								for m in range(cloud_upper_layer_base_idx,cloud_upper_layer_top_idx+1):
									cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height[m] = cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height[m] + 1 

								for m in range(cloud_lower_layer_base_idx,cloud_lower_layer_top_idx+1):
									cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height[m] = cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height[m] + 1 
							
	plt.barh(cloud_layer_count_bar_x, cldclass_lidar_two_cloud_layer_bar_height, width, color="lightcoral", alpha=1.0, align='center', label='nb cloud layer = 2')
	plt.barh(cloud_layer_count_bar_x, cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height, width, color="royalblue", alpha=1.0, align='center', label='nb cloud layer = 2 + caliop 5km cot')

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height[j],cldclass_lidar_two_cloud_layer_and_caliop_5km_cot_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="royalblue")
		
	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cldclass_lidar_two_cloud_layer_bar_height[j],cldclass_lidar_two_cloud_layer_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="lightcoral")				

	plt.legend()
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	bins = np.linspace(0, 5, 20)

	plt.hist(caliop_5km_cloud_optical_thickness_hist, bins=bins, color="royalblue", alpha=0.5, normed=1)

	plt.title('Caliop 5km Cloud Optical Thickness \n (cldclass-lidar nb cloud layer = 2)',fontsize=10)

	plt.xlabel('Cloud Optical Thichness',fontsize=10)
	plt.ylabel('Probability Density Function',fontsize=10)

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# cldclass-lidar (cloud phase) Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def cldclass_lidar_cloud_phase_layer_vertical_distribution_bar(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]
	
	width = 0.35  
	
	cloud_layer_count_bar_x = [i*0.1 + 0.05 for i in range(200)]
	cloud_layer_count_bar_height = [0 for i in range(200)]

	cloud_layer_ice_phase_count_bar_height = [0 for i in range(200)]
	cloud_layer_liq_phase_count_bar_height = [0 for i in range(200)]
	cloud_layer_mixed_phase_count_bar_height = [0 for i in range(200)]	
		
	for i in range(data.shape[0]):
	#for i in range(100):
		
		nb_cloud_layer = int( data[i,15] )
		
		cloud_layer_parameters = np.zeros((10,5))
		
		m = 19
		if nb_cloud_layer <= 10:
			for j in range(nb_cloud_layer):
				cloud_layer_parameters[j,0] = data[i,m]
				cloud_layer_parameters[j,1] = data[i,m+1]
				cloud_layer_parameters[j,2] = data[i,m+2]
				cloud_layer_parameters[j,3] = data[i,m+3]
				cloud_layer_parameters[j,4] = data[i,m+4]
				m = m + 5
			
			#print('nb_cloud_layer --->', nb_cloud_layer)
			#for j in range(nb_cloud_layer):
			#	print( [ "{0:0.1f}".format(value) for value in list(cloud_layer_parameters[j,:]) ] )
			#print('------')
					
			for j in range(nb_cloud_layer):
				if cloud_layer_parameters[j,0] > 0.0 and cloud_layer_parameters[j,1] > 0.0:
					cloud_layer_top_idx = int(cloud_layer_parameters[j,0]/0.1)
					cloud_layer_base_idx = int(cloud_layer_parameters[j,1]/0.1)
					if int(cloud_layer_parameters[j,4]) == 3: # cldclass-lidar products
						if eval(dict_settings['cloud_layer_logical_expression'] ):
							for m in range(cloud_layer_base_idx,cloud_layer_top_idx+1):
								cloud_layer_count_bar_height[m] = cloud_layer_count_bar_height[m] + 1
								if int(cloud_layer_parameters[j,3]) == 1: cloud_layer_ice_phase_count_bar_height[m] = cloud_layer_ice_phase_count_bar_height[m] + 1 # ice
								if int(cloud_layer_parameters[j,3]) == 3: cloud_layer_liq_phase_count_bar_height[m] = cloud_layer_liq_phase_count_bar_height[m] + 1 # liquid
								if int(cloud_layer_parameters[j,3]) == 2: cloud_layer_mixed_phase_count_bar_height[m] = cloud_layer_mixed_phase_count_bar_height[m] + 1 # mixed	
	
	plt.barh(cloud_layer_count_bar_x, cloud_layer_liq_phase_count_bar_height, width, color="lightcoral", alpha=1.0, align='center', label='liquid phase')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_ice_phase_count_bar_height, width, color="royalblue", alpha=1.0, align='center', label='ice phase')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_mixed_phase_count_bar_height, width, color="wheat", alpha=0.8, align='center', label='mixed phase')

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_mixed_phase_count_bar_height[j],cloud_layer_mixed_phase_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="wheat")

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_liq_phase_count_bar_height[j],cloud_layer_liq_phase_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="lightcoral")
		
	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_ice_phase_count_bar_height[j],cloud_layer_ice_phase_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="royalblue")				

	plt.legend()
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# cldclass-lidar (nb cloud layer) vs modis -- Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def cldclass_lidar_vs_modis_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]
	
	width = 0.35  
	
	cloud_layer_count_bar_x = [i*0.1 + 0.05 for i in range(200)]
	cloud_layer_count_bar_height = [0 for i in range(200)]

	cloud_layer_mono_count_bar_height = [0 for i in range(200)]
	cloud_layer_multi_count_bar_height = [0 for i in range(200)]
	cloud_layer_two_count_bar_height = [0 for i in range(200)]	
		
	for i in range(data.shape[0]):
	#for i in range(100):

		modis_multilayer_sds = data[i,0]
		modis_cloud_optical_tickness = data[i,13] * 0.01

		if 'modis_without_pavolonis_heidinger_test' in dict_settings:
			if dict_settings['modis_without_pavolonis_heidinger_test']:
				if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3		

		if modis_multilayer_sds == 1 : modis_multilayer_flag = 0 
		if modis_multilayer_sds > 1: modis_multilayer_flag = 1
		if modis_cloud_optical_tickness < 4.0: modis_multilayer_flag = 2
				
		nb_cloud_layer = int( data[i,15] )
		
		cloud_layer_parameters = np.zeros((10,5))
		
		m = 19
		if nb_cloud_layer <= 10:
			for j in range(nb_cloud_layer):
				cloud_layer_parameters[j,0] = data[i,m]
				cloud_layer_parameters[j,1] = data[i,m+1]
				cloud_layer_parameters[j,2] = data[i,m+2]
				cloud_layer_parameters[j,3] = data[i,m+3]
				cloud_layer_parameters[j,4] = data[i,m+4]
				m = m + 5
			
			#print('nb_cloud_layer --->', nb_cloud_layer)
			#for j in range(nb_cloud_layer):
			#	print( [ "{0:0.1f}".format(value) for value in list(cloud_layer_parameters[j,:]) ] )
			#print('------')
					
			for j in range(nb_cloud_layer):
				if cloud_layer_parameters[j,0] > 0.0 and cloud_layer_parameters[j,1] > 0.0:
					cloud_layer_top_idx = int(cloud_layer_parameters[j,0]/0.1)
					cloud_layer_base_idx = int(cloud_layer_parameters[j,1]/0.1)
					if int(cloud_layer_parameters[j,4]) == 3: # cldclass-lidar products
						if modis_multilayer_flag == dict_settings['modis_multilayer_flag']:
							for m in range(cloud_layer_base_idx,cloud_layer_top_idx+1):
								cloud_layer_count_bar_height[m] = cloud_layer_count_bar_height[m] + 1
								if int(data[i,18]) == 1: cloud_layer_mono_count_bar_height[m] = cloud_layer_mono_count_bar_height[m] + 1 
								if int(data[i,18]) > 1 : cloud_layer_multi_count_bar_height[m] = cloud_layer_multi_count_bar_height[m] + 1 
								if int(data[i,18]) == 2: cloud_layer_two_count_bar_height[m] = cloud_layer_two_count_bar_height[m] + 1 	
	
	plt.barh(cloud_layer_count_bar_x, cloud_layer_multi_count_bar_height, width, color="lightcoral", alpha=1.0, align='center', label='nb cloud layer > 1')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_mono_count_bar_height, width, color="royalblue", alpha=1.0, align='center', label='nb cloud layer = 1')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_two_count_bar_height, width, color="wheat", alpha=0.8, align='center', label='nb cloud layer = 2')

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_two_count_bar_height[j],cloud_layer_two_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="wheat")

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_multi_count_bar_height[j],cloud_layer_multi_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="lightcoral")
		
	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_mono_count_bar_height[j],cloud_layer_mono_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="royalblue")				

	plt.legend()
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# cldclass-lidar vs modis (cloud phase) Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def modis_vs_cldclass_lidar_cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]
	
	width = 0.35  
	
	cloud_layer_count_bar_x = [i*0.1 + 0.05 for i in range(200)]
	cloud_layer_count_bar_height = [0 for i in range(200)]

	cloud_layer_ice_phase_count_bar_height = [0 for i in range(200)]
	cloud_layer_liq_phase_count_bar_height = [0 for i in range(200)]
	cloud_layer_mixed_phase_count_bar_height = [0 for i in range(200)]	
		
	for i in range(data.shape[0]):
	#for i in range(100):

		modis_cloud_phase_sds = int( data[i,6] )
		
		nb_cloud_layer = int( data[i,15] )
		
		cloud_layer_parameters = np.zeros((10,5))
		
		m = 19
		if nb_cloud_layer <= 10 and modis_cloud_phase_sds == dict_settings['modis_cloud_phase_flag']:

			for j in range(nb_cloud_layer):
				cloud_layer_parameters[j,0] = data[i,m]
				cloud_layer_parameters[j,1] = data[i,m+1]
				cloud_layer_parameters[j,2] = data[i,m+2]
				cloud_layer_parameters[j,3] = data[i,m+3]
				cloud_layer_parameters[j,4] = data[i,m+4]
				m = m + 5
			
			#print('nb_cloud_layer --->', nb_cloud_layer)
			#for j in range(nb_cloud_layer):
			#	print( [ "{0:0.1f}".format(value) for value in list(cloud_layer_parameters[j,:]) ] )
			#print('------')
					
			for j in range(nb_cloud_layer):
				if cloud_layer_parameters[j,0] > 0.0 and cloud_layer_parameters[j,1] > 0.0:
					cloud_layer_top_idx = int(cloud_layer_parameters[j,0]/0.1)
					cloud_layer_base_idx = int(cloud_layer_parameters[j,1]/0.1)
					if int(cloud_layer_parameters[j,4]) == 3: # cldclass-lidar products
						if eval(dict_settings['cloud_layer_logical_expression'] ):
							for m in range(cloud_layer_base_idx,cloud_layer_top_idx+1):
								cloud_layer_count_bar_height[m] = cloud_layer_count_bar_height[m] + 1
								if int(cloud_layer_parameters[j,3]) == 1: cloud_layer_ice_phase_count_bar_height[m] = cloud_layer_ice_phase_count_bar_height[m] + 1 # ice
								if int(cloud_layer_parameters[j,3]) == 3: cloud_layer_liq_phase_count_bar_height[m] = cloud_layer_liq_phase_count_bar_height[m] + 1 # liquid
								if int(cloud_layer_parameters[j,3]) == 2: cloud_layer_mixed_phase_count_bar_height[m] = cloud_layer_mixed_phase_count_bar_height[m] + 1 # mixed	
	
	plt.barh(cloud_layer_count_bar_x, cloud_layer_liq_phase_count_bar_height, width, color="lightcoral", alpha=1.0, align='center', label='liquid phase')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_ice_phase_count_bar_height, width, color="royalblue", alpha=1.0, align='center', label='ice phase')
	plt.barh(cloud_layer_count_bar_x, cloud_layer_mixed_phase_count_bar_height, width, color="wheat", alpha=0.8, align='center', label='mixed phase')

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_mixed_phase_count_bar_height[j],cloud_layer_mixed_phase_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="wheat")

	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_liq_phase_count_bar_height[j],cloud_layer_liq_phase_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="lightcoral")
		
	for j in range(len(cloud_layer_count_bar_x)-1):
		xpt = (cloud_layer_ice_phase_count_bar_height[j],cloud_layer_ice_phase_count_bar_height[j+1])
		ypt = (cloud_layer_count_bar_x[j], cloud_layer_count_bar_x[j+1])
		plt.plot(xpt, ypt, color="royalblue")				

	plt.legend()
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#

np.set_printoptions(precision=2)

def modis_cloud_effective_radius_hist_vs_cldclass_lidar_monolayer_cloud(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	modis_cloud_effective_radius_hist_mono_same_phase = []
	modis_cloud_effective_radius_hist_mono_diff_phase = []

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
		
			if 'modis_without_pavolonis_heidinger_test' in dict_settings:
				if dict_settings['modis_without_pavolonis_heidinger_test']:
					if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
		
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
			
			if cldclass_lidar_nb_cloud_layer > 0:

				modis_cloud_phase = data[i,6]
				modis_cloud_effective_radius = data[i,dict_settings['modis_cloud_effective_radius_idx']] * 0.01
				
				if eval(dict_settings['modis_cloud_phase_logical_expression']):

					if modis_cloud_effective_radius > 0.0:
						if cldclass_lidar_nb_cloud_layer == 1:
							if cldclass_lidar_phase_list[0] == 1 and modis_cloud_phase == 3: # both ice
								modis_cloud_effective_radius_hist_mono_same_phase.append(modis_cloud_effective_radius)
							if cldclass_lidar_phase_list[0] == 3 and modis_cloud_phase == 2: # both liquid
								modis_cloud_effective_radius_hist_mono_same_phase.append(modis_cloud_effective_radius)							
							if cldclass_lidar_phase_list[0] == 3 and modis_cloud_phase == 3: # liquid vs ice
								modis_cloud_effective_radius_hist_mono_diff_phase.append(modis_cloud_effective_radius)
							if cldclass_lidar_phase_list[0] == 1 and modis_cloud_phase == 2: # ice vs liquid
								modis_cloud_effective_radius_hist_mono_diff_phase.append(modis_cloud_effective_radius)

	mono_mean = sum(modis_cloud_effective_radius_hist_mono_same_phase) * 1.0 / len(modis_cloud_effective_radius_hist_mono_same_phase)
	multi_mean = sum(modis_cloud_effective_radius_hist_mono_diff_phase) * 1.0 / len(modis_cloud_effective_radius_hist_mono_diff_phase)
					
	bins = np.linspace(0, 60, 60)

	plt.hist(modis_cloud_effective_radius_hist_mono_same_phase, bins=bins, color="royalblue", alpha=0.5, normed=1)
	plt.hist(modis_cloud_effective_radius_hist_mono_diff_phase, bins=bins, color="lightcoral", alpha=0.5, normed=1)

	plt.title(dict_settings['fig_title'],fontsize=10)

	plt.xlabel(dict_settings['fig_xlabel'],fontsize=8)
	plt.ylabel(dict_settings['fig_ylabel'],fontsize=8)


	labels= [dict_settings['fig_bar_1_label'],dict_settings['fig_bar_2_label']]

	plt.legend(labels)

	plt.xlim(0,60)

	plt.grid()		

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# modis cloud effective radius vs cldclass-lidar histograms

np.set_printoptions(precision=2)

def modis_cloud_effective_radius_vs_cldclass_lidar_histograms(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	modis_cloud_effective_radius_hist_mono = []
	modis_cloud_effective_radius_hist_multi = []

	nb_tot = 0

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01

		if 'modis_without_pavolonis_heidinger_test' in dict_settings:
			if dict_settings['modis_without_pavolonis_heidinger_test']:
				if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
			
		if eval(dict_settings['modis_cloud_optical_thickness_expression']):
			
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
			
			if cldclass_lidar_nb_cloud_layer > 0:

				modis_cloud_phase = data[i,6]
				modis_cloud_effective_radius = data[i,dict_settings['modis_cloud_effective_radius_idx']] * 0.01
				
				if eval(dict_settings['modis_cloud_phase_logical_expression']):
					nb_tot = nb_tot + 1
					if modis_cloud_effective_radius > 0.0:
						if cldclass_lidar_nb_cloud_layer == 1:
							if cldclass_lidar_phase_list[0] == 1 and modis_cloud_phase == 3: # both ice
								modis_cloud_effective_radius_hist_mono.append(modis_cloud_effective_radius)
							if cldclass_lidar_phase_list[0] == 3 and modis_cloud_phase == 2: # both liquid
								modis_cloud_effective_radius_hist_mono.append(modis_cloud_effective_radius)							
						if cldclass_lidar_nb_cloud_layer > 1:
							modis_cloud_effective_radius_hist_multi.append(modis_cloud_effective_radius)

	mono_mean = sum(modis_cloud_effective_radius_hist_mono) * 1.0 / len(modis_cloud_effective_radius_hist_mono)
	multi_mean = sum(modis_cloud_effective_radius_hist_multi) * 1.0 / len(modis_cloud_effective_radius_hist_multi)
					
	bins = np.linspace(0, 60, 60)

	plt.hist(modis_cloud_effective_radius_hist_mono, bins=bins, color="royalblue", alpha=0.5, normed=1)
	plt.hist(modis_cloud_effective_radius_hist_multi, bins=bins, color="lightcoral", alpha=0.5, normed=1)

	plt.title(dict_settings['fig_title'],fontsize=10)

	plt.xlabel(dict_settings['fig_xlabel'],fontsize=10)
	plt.ylabel('Probability Density Function',fontsize=10)

	labels= ['cldlclass-lidar monolayer clouds \n' + dict_settings['modis_cloud_phase_label'] +' cloud phase \n (mean: ' + str(round(mono_mean,2)) + ')',"cldlclass-lidar multilayer clouds \n (mean: " + str(round(multi_mean,2)) + ')']

	plt.legend(labels)

	plt.xlim(0,60)

	plt.grid()		

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()
	
	print('nb_tot ---->', nb_tot)
	
	return fig_count

#----------------------------------------------------------------------------------------#
# modis cloud optical thickness vs cldclass-lidar histograms

np.set_printoptions(precision=2)

def modis_cloud_optical_thickness_vs_cldclass_lidar_histograms(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	modis_cloud_optical_thickness_hist_mono = []
	modis_cloud_optical_thickness_hist_multi = []

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
		
			if 'modis_without_pavolonis_heidinger_test' in dict_settings:
				if dict_settings['modis_without_pavolonis_heidinger_test']:
					if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
		
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
			
			if cldclass_lidar_nb_cloud_layer > 0:

				modis_cloud_phase = data[i,6]
				
				if eval(dict_settings['modis_cloud_phase_logical_expression']):
					if modis_cloud_optical_tickness > 0.0:
						if cldclass_lidar_nb_cloud_layer == 1:
							if cldclass_lidar_phase_list[0] == 1 and modis_cloud_phase == 3: # both ice
								modis_cloud_optical_thickness_hist_mono.append(modis_cloud_optical_tickness)
							if cldclass_lidar_phase_list[0] == 3 and modis_cloud_phase == 2: # both liquid
								modis_cloud_optical_thickness_hist_mono.append(modis_cloud_optical_tickness)							
						if cldclass_lidar_nb_cloud_layer > 1:
							modis_cloud_optical_thickness_hist_multi.append(modis_cloud_optical_tickness)

	mono_mean = sum(modis_cloud_optical_thickness_hist_mono) * 1.0 / len(modis_cloud_optical_thickness_hist_mono)
	multi_mean = sum(modis_cloud_optical_thickness_hist_multi) * 1.0 / len(modis_cloud_optical_thickness_hist_multi)
					
	bins = np.linspace(0, 100, 50)

	plt.hist(modis_cloud_optical_thickness_hist_mono, bins=bins, color="royalblue", alpha=0.5, normed=1)
	plt.hist(modis_cloud_optical_thickness_hist_multi, bins=bins, color="lightcoral", alpha=0.5, normed=1)

	plt.title(dict_settings['fig_title'],fontsize=10)

	plt.xlabel(dict_settings['fig_xlabel'],fontsize=10)
	plt.ylabel('Probability Density Function',fontsize=10)


	labels= ['cldlclass-lidar monolayer clouds \n' + dict_settings['modis_cloud_phase_label'] +' cloud phase \n (mean: ' + str(round(mono_mean,2)) + ')',"cldlclass-lidar multilayer clouds \n (mean: " + str(round(multi_mean,2)) + ')']

	plt.legend(labels)

	plt.xlim(0,50)

	plt.grid()		

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# modis cloud effective radius diff (modis & cldclass-lidar multilayer flag matrices)

np.set_printoptions(precision=2)

def modis_cloud_effective_radius_vs_cldclass_lidar_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	modis_cloud_effective_radius_mean_mono = np.zeros((3,3))
	modis_cloud_effective_radius_mean_multi = np.zeros((3,3))

	modis_cloud_effective_radius_mean_mono_count = np.zeros((3,3))
	modis_cloud_effective_radius_mean_multi_count = np.zeros((3,3))

	conf_arr = np.zeros((3,3))

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
		
			if 'modis_without_pavolonis_heidinger_test' in dict_settings:
				if dict_settings['modis_without_pavolonis_heidinger_test']:
					if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
		
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
	
			if cldclass_lidar_nb_cloud_layer > 0:	
		
				modis_multilayer_flag = -1
				if modis_multilayer_sds == 1 : modis_multilayer_flag = 0 
				if modis_multilayer_sds > 1: modis_multilayer_flag = 1
				if modis_cloud_optical_tickness < 4.0: modis_multilayer_flag = 2

				modis_cloud_phase = data[i,6]
			
				modis_cloud_effective_radius_21 = data[i,10] * 0.01
				modis_cloud_effective_radius_16 = data[i,11] * 0.01
				modis_cloud_effective_radius_37 = data[i,12] * 0.01

				if eval(dict_settings['modis_cloud_phase_logical_expression']):
					if modis_multilayer_flag >= 0:
						for idx,modis_cloud_effective_radius_idx in enumerate([10,11,12]):
							modis_cloud_effective_radius = data[i,modis_cloud_effective_radius_idx] * 0.01					
							if modis_cloud_effective_radius > 0.0:
								if cldclass_lidar_nb_cloud_layer == 1:
									if cldclass_lidar_phase_list[0] == 1 and modis_cloud_phase == 3: # both ice
										modis_cloud_effective_radius_mean_mono[modis_multilayer_flag,idx] = modis_cloud_effective_radius_mean_mono[modis_multilayer_flag,idx] + modis_cloud_effective_radius
										modis_cloud_effective_radius_mean_mono_count[modis_multilayer_flag,idx] = modis_cloud_effective_radius_mean_mono_count[modis_multilayer_flag,idx] + 1
									if cldclass_lidar_phase_list[0] == 3 and modis_cloud_phase == 2: # both liquid
										modis_cloud_effective_radius_mean_mono[modis_multilayer_flag,idx] = modis_cloud_effective_radius_mean_mono[modis_multilayer_flag,idx] + modis_cloud_effective_radius
										modis_cloud_effective_radius_mean_mono_count[modis_multilayer_flag,idx] = modis_cloud_effective_radius_mean_mono_count[modis_multilayer_flag,idx] + 1						
								if cldclass_lidar_nb_cloud_layer > 1:
									modis_cloud_effective_radius_mean_multi[modis_multilayer_flag,idx] = modis_cloud_effective_radius_mean_multi[modis_multilayer_flag,idx] + modis_cloud_effective_radius
									modis_cloud_effective_radius_mean_multi_count[modis_multilayer_flag,idx] = modis_cloud_effective_radius_mean_multi_count[modis_multilayer_flag,idx] + 1

	modis_cloud_effective_radius_mean_mono[:,:] = modis_cloud_effective_radius_mean_mono[:,:] / modis_cloud_effective_radius_mean_mono_count[:,:]
	modis_cloud_effective_radius_mean_multi[:,:] = modis_cloud_effective_radius_mean_multi[:,:] / modis_cloud_effective_radius_mean_multi_count[:,:]
		
	conf_arr[:,:] = modis_cloud_effective_radius_mean_multi[:,:] - modis_cloud_effective_radius_mean_mono[:,:]
	
	print(conf_arr)
		
	df_cm = pd.DataFrame(conf_arr, 
					  index = [ 'modis C6 Undet.', 'modis C6 monolayer cloud', 'modis C6 multilayer cloud'],
					  columns = ['modis c6 \n' + r'CER ($2.1\mu m$)', 'modis c6 \n' + r'CER ($1.6\mu m$)', 'modis c6 \n' + r'CER ($3.7\mu m$)'])

	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)
	
	
	label_modis_mono_1 = [modis_cloud_effective_radius_mean_mono[0,0],modis_cloud_effective_radius_mean_mono[0,1],modis_cloud_effective_radius_mean_mono[0,2]]
	label_modis_mono_2 = [modis_cloud_effective_radius_mean_mono[1,0],modis_cloud_effective_radius_mean_mono[1,1],modis_cloud_effective_radius_mean_mono[1,2]]
	label_modis_mono_3 = [modis_cloud_effective_radius_mean_mono[2,0],modis_cloud_effective_radius_mean_mono[2,1],modis_cloud_effective_radius_mean_mono[2,2]]
	
	label_modis_mono_1 = [round(i,1) for i in label_modis_mono_1]
	label_modis_mono_2 = [round(i,1) for i in label_modis_mono_2]
	label_modis_mono_3 = [round(i,1) for i in label_modis_mono_3]

	label_modis_multi_1 = [modis_cloud_effective_radius_mean_multi[0,0],modis_cloud_effective_radius_mean_multi[0,1],modis_cloud_effective_radius_mean_multi[0,2]]
	label_modis_multi_2 = [modis_cloud_effective_radius_mean_multi[1,0],modis_cloud_effective_radius_mean_multi[1,1],modis_cloud_effective_radius_mean_multi[1,2]]
	label_modis_multi_3 = [modis_cloud_effective_radius_mean_multi[2,0],modis_cloud_effective_radius_mean_multi[2,1],modis_cloud_effective_radius_mean_multi[2,2]]
	
	label_modis_multi_1 = [round(i,1) for i in label_modis_multi_1]
	label_modis_multi_2 = [round(i,1) for i in label_modis_multi_2]
	label_modis_multi_3 = [round(i,1) for i in label_modis_multi_3]

	label_modis_1 = ['mono ' + str(label_modis_mono_1[idx]) + '\n multi ' + str(label_modis_multi_1[idx]) for idx in range(3) ]
	label_modis_2 = ['mono ' + str(label_modis_mono_2[idx]) + '\n multi ' + str(label_modis_multi_2[idx]) for idx in range(3) ]
	label_modis_3 = ['mono ' + str(label_modis_mono_3[idx]) + '\n multi ' + str(label_modis_multi_3[idx]) for idx in range(3) ]
	
	labels =  np.array([label_modis_1,label_modis_2,label_modis_3])

	if dict_settings['modis_cloud_phase_label'] == 'liquid':
		vmin = -4.0
		vmax = 4.0
	else:
		vmin = -10.0
		vmax = 10.0

	res = sn.heatmap(df_cm, annot = labels, vmin=vmin, vmax=vmax, fmt='', cmap="RdBu_r", cbar_kws={'label': 'CER (multi) - CER (mono)'})
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], ['modis C6 \n Undet. (COT < 4.0) ', 'modis C6 \n multilayer cloud', 'modis C6 \n monolayer cloud'],va='center')

	plt.title(dict_settings['fig_title'] + str())

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# modis cloud effective radius diff matrices (cldclass-lidar: 2 layers, heterogeneous phase)

np.set_printoptions(precision=2)

def modis_cloud_effective_radius_vs_cldclass_lidar_heterogeneous_phase_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	modis_cloud_effective_radius_mean_mono = np.zeros((3))
	modis_cloud_effective_radius_mean_multi = np.zeros((3,3))

	modis_cloud_effective_radius_mean_mono_count = np.zeros((3))
	modis_cloud_effective_radius_mean_multi_count = np.zeros((3,3))

	conf_arr = np.zeros((3,3))

	cldclass_lidar_profile = dict_settings['cldclass_lidar_profile']
	
	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
		
			if 'modis_without_pavolonis_heidinger_test' in dict_settings:
				if dict_settings['modis_without_pavolonis_heidinger_test']:
					if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
		
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)

			if cldclass_lidar_nb_cloud_layer > 0:	
				
				cldclass_lidar_cloud_phase_flag = -1
				if cldclass_lidar_nb_cloud_layer == 2:
				
					cldclass_lidar_phase_upper_layer = cldclass_lidar_phase_list[0]
					cldclass_lidar_phase_lower_layer = cldclass_lidar_phase_list[1]
				
					if cldclass_lidar_phase_upper_layer == cldclass_lidar_profile[0][0] and cldclass_lidar_phase_lower_layer == cldclass_lidar_profile[0][1]: cldclass_lidar_cloud_phase_flag = 0
					if cldclass_lidar_phase_upper_layer == cldclass_lidar_profile[1][0] and cldclass_lidar_phase_lower_layer == cldclass_lidar_profile[1][1]: cldclass_lidar_cloud_phase_flag = 1
					if cldclass_lidar_phase_upper_layer == cldclass_lidar_profile[2][0] and cldclass_lidar_phase_lower_layer == cldclass_lidar_profile[2][1]: cldclass_lidar_cloud_phase_flag = 2

				modis_cloud_phase = data[i,6]
			
				if eval(dict_settings['modis_cloud_phase_logical_expression']):
					if cldclass_lidar_nb_cloud_layer == 1:
						for idx,modis_cloud_effective_radius_idx in enumerate([10,11,12]):
							modis_cloud_effective_radius = data[i,modis_cloud_effective_radius_idx] * 0.01					
							if modis_cloud_effective_radius > 0.0:
								if cldclass_lidar_phase_list[0] == 1 and modis_cloud_phase == 3: # both ice
									modis_cloud_effective_radius_mean_mono[idx] = modis_cloud_effective_radius_mean_mono[idx] + modis_cloud_effective_radius
									modis_cloud_effective_radius_mean_mono_count[idx] = modis_cloud_effective_radius_mean_mono_count[idx] + 1
								if cldclass_lidar_phase_list[0] == 3 and modis_cloud_phase == 2: # both liquid
									modis_cloud_effective_radius_mean_mono[idx] = modis_cloud_effective_radius_mean_mono[idx] + modis_cloud_effective_radius
									modis_cloud_effective_radius_mean_mono_count[idx] = modis_cloud_effective_radius_mean_mono_count[idx] + 1


				if eval(dict_settings['modis_cloud_phase_logical_expression']):
					if cldclass_lidar_cloud_phase_flag >= 0:
						for idx,modis_cloud_effective_radius_idx in enumerate([10,11,12]):
							modis_cloud_effective_radius = data[i,modis_cloud_effective_radius_idx] * 0.01					
							if modis_cloud_effective_radius > 0.0:
								if cldclass_lidar_nb_cloud_layer == 2:
									modis_cloud_effective_radius_mean_multi[cldclass_lidar_cloud_phase_flag,idx] = modis_cloud_effective_radius_mean_multi[cldclass_lidar_cloud_phase_flag,idx] + modis_cloud_effective_radius
									modis_cloud_effective_radius_mean_multi_count[cldclass_lidar_cloud_phase_flag,idx] = modis_cloud_effective_radius_mean_multi_count[cldclass_lidar_cloud_phase_flag,idx] + 1

	modis_cloud_effective_radius_mean_mono[:] = modis_cloud_effective_radius_mean_mono[:] / modis_cloud_effective_radius_mean_mono_count[:]
	modis_cloud_effective_radius_mean_multi[:,:] = modis_cloud_effective_radius_mean_multi[:,:] / modis_cloud_effective_radius_mean_multi_count[:,:]

	print(modis_cloud_effective_radius_mean_mono)
	print(modis_cloud_effective_radius_mean_mono_count)
	
	print(modis_cloud_effective_radius_mean_multi)
	print(modis_cloud_effective_radius_mean_multi_count)
	
	print('--------')
	
	for j in range(3):
		conf_arr[:,j] = modis_cloud_effective_radius_mean_multi[:,j] - modis_cloud_effective_radius_mean_mono[j]
	
	print(conf_arr)
		
	df_cm = pd.DataFrame(conf_arr, 
					  index = [ 'modis C6 Undet.', 'modis C6 monolayer cloud', 'modis C6 multilayer cloud'],
					  columns = ['modis c6 \n' + r'CER ($2.1\mu m$)', 'modis c6 \n' + r'CER ($1.6\mu m$)', 'modis c6 \n' + r'CER ($3.7\mu m$)'])

	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)
	
	
	label_modis_mono_1 = [modis_cloud_effective_radius_mean_mono[0],modis_cloud_effective_radius_mean_mono[1],modis_cloud_effective_radius_mean_mono[2]]
	label_modis_mono_2 = [modis_cloud_effective_radius_mean_mono[0],modis_cloud_effective_radius_mean_mono[1],modis_cloud_effective_radius_mean_mono[2]]
	label_modis_mono_3 = [modis_cloud_effective_radius_mean_mono[0],modis_cloud_effective_radius_mean_mono[1],modis_cloud_effective_radius_mean_mono[2]]
	
	label_modis_mono_1 = [round(i,1) for i in label_modis_mono_1]
	label_modis_mono_2 = [round(i,1) for i in label_modis_mono_2]
	label_modis_mono_3 = [round(i,1) for i in label_modis_mono_3]

	label_modis_multi_1 = [modis_cloud_effective_radius_mean_multi[0,0],modis_cloud_effective_radius_mean_multi[0,1],modis_cloud_effective_radius_mean_multi[0,2]]
	label_modis_multi_2 = [modis_cloud_effective_radius_mean_multi[1,0],modis_cloud_effective_radius_mean_multi[1,1],modis_cloud_effective_radius_mean_multi[1,2]]
	label_modis_multi_3 = [modis_cloud_effective_radius_mean_multi[2,0],modis_cloud_effective_radius_mean_multi[2,1],modis_cloud_effective_radius_mean_multi[2,2]]
	
	label_modis_multi_1 = [round(i,1) for i in label_modis_multi_1]
	label_modis_multi_2 = [round(i,1) for i in label_modis_multi_2]
	label_modis_multi_3 = [round(i,1) for i in label_modis_multi_3]

	label_modis_1 = ['mono ' + str(label_modis_mono_1[idx]) + '\n multi ' + str(label_modis_multi_1[idx]) for idx in range(3) ]
	label_modis_2 = ['mono ' + str(label_modis_mono_2[idx]) + '\n multi ' + str(label_modis_multi_2[idx]) for idx in range(3) ]
	label_modis_3 = ['mono ' + str(label_modis_mono_3[idx]) + '\n multi ' + str(label_modis_multi_3[idx]) for idx in range(3) ]
	
	labels =  np.array([label_modis_1,label_modis_2,label_modis_3])

	if dict_settings['modis_cloud_phase_label'] == 'liquid':
		vmin = -4.0
		vmax = 4.0
	else:
		vmin = -10.0
		vmax = 10.0

	res = sn.heatmap(df_cm, annot = labels, vmin=vmin, vmax=vmax, fmt='', cmap="RdBu_r", cbar_kws={'label': 'CER (multi) - CER (mono)'})
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], dict_settings['fig_yticks'],va='center')

	plt.title(dict_settings['fig_title'] + str())

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# zonal plots

def zonal_plots(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	latitude_bin = [i for i in range(-90,91)]
	
	modis_and_cldclass_lidar_cloud_counts_by_latitude_bin = [0 for i in latitude_bin]
	
	modis_multilayer_cloud_counts_by_latitude_bin = [0 for i in latitude_bin]
	cldclass_lidar_multilayer_cloud_counts_by_latitude_bin = [0 for i in latitude_bin]
	
	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]
		modis_cloud_optical_tickness = data[i,13] * 0.01
		modis_latitude = data[i,7]	
		
		latitude_idx = int(modis_latitude + 90.0)
		
		if 'modis_without_pavolonis_heidinger_test' in dict_settings:
			if dict_settings['modis_without_pavolonis_heidinger_test']:
				if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3

		caliop_5km_cloud_layer_top_list = []
		caliop_5km_cloud_layer_base_list = []
		caliop_5km_cloud_layer_cot_list = []
		caliop_5km_cloud_layer_phase_list = []

		cldclass_lidar_cloud_layer_top_list = []
		cldclass_lidar_cloud_layer_base_list = []
		cldclass_lidar_cloud_layer_phase_list = []
	
		cloud_layer_idx = 19
		for n in range(10):
			if int(data[i,cloud_layer_idx+4]) == 2: # caliop 5km
				caliop_5km_cloud_layer_top_list.append(data[i,cloud_layer_idx])
				caliop_5km_cloud_layer_base_list.append(data[i,cloud_layer_idx+1])
				caliop_5km_cloud_layer_cot_list.append(data[i,cloud_layer_idx+2])												
				caliop_5km_cloud_layer_phase_list.append(int(data[i,cloud_layer_idx+3]))
			if int(data[i,cloud_layer_idx+4]) == 3: # cldclass-lidar
				cldclass_lidar_cloud_layer_top_list.append(data[i,cloud_layer_idx])
				cldclass_lidar_cloud_layer_base_list.append(data[i,cloud_layer_idx+1])											
				cldclass_lidar_cloud_layer_phase_list.append(int(data[i,cloud_layer_idx+3]))
			cloud_layer_idx = cloud_layer_idx + 5	

		cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_cloud_layer_phase_list)
					
		if cldclass_lidar_nb_cloud_layer > 0:

			modis_and_cldclass_lidar_cloud_counts_by_latitude_bin[latitude_idx] = modis_and_cldclass_lidar_cloud_counts_by_latitude_bin[latitude_idx] + 1

			if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:	 

				cloud_layer_separation_distance = 0.0
				if cldclass_lidar_nb_cloud_layer > 1:
					cloud_layer_separation_distance = cldclass_lidar_cloud_layer_base_list[0] - cldclass_lidar_cloud_layer_top_list[1]
							
				if modis_multilayer_sds > 1 : 
					modis_multilayer_cloud_counts_by_latitude_bin[latitude_idx] = modis_multilayer_cloud_counts_by_latitude_bin[latitude_idx] + 1 
			
				if cldclass_lidar_nb_cloud_layer > 1 and cloud_layer_separation_distance > dict_settings['distance_threshold']:
					cldclass_lidar_multilayer_cloud_counts_by_latitude_bin[latitude_idx] = cldclass_lidar_multilayer_cloud_counts_by_latitude_bin[latitude_idx] + 1

	bins = np.linspace(-90, 90, 90)
	
	for idx,i in enumerate(modis_and_cldclass_lidar_cloud_counts_by_latitude_bin):
		if i > 0.0:
			modis_multilayer_cloud_counts_by_latitude_bin[idx] = 1.0 * modis_multilayer_cloud_counts_by_latitude_bin[idx] / i
			cldclass_lidar_multilayer_cloud_counts_by_latitude_bin[idx] = 1.0 * cldclass_lidar_multilayer_cloud_counts_by_latitude_bin[idx] / i 
	
	plt.bar([i for i in latitude_bin], modis_multilayer_cloud_counts_by_latitude_bin, 1.0, color="lightcoral", alpha=0.5)

	plt.bar([i for i in latitude_bin], cldclass_lidar_multilayer_cloud_counts_by_latitude_bin, 1.0, color="royalblue", alpha=0.5)

	plt.xlim(-90,90)
	
	plt.xlabel('Latitude')
	plt.ylabel(dict_settings['fig_y_label'])
	
	plt.title(dict_settings['fig_title'])
	
	labels= [dict_settings['fig_bar_1_label'], dict_settings['fig_bar_2_label']]
	plt.legend(labels)
	
	#plt.ylim(0,1400*1.0)
	plt.ylim(0,1.0)
	
	plt.grid(True)
	
	fig_count = fig_count + 1
	
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

Weighted random sampling with a reservoir example in python

					
from random import uniform
from random import randint

import numpy as np
import matplotlib.pyplot as plt

# Reference: 
# Weighted random sampling with a reservoir; Pavlos S.Efraimidis; Paul G. Spirakis (2006)

#----------------------------------------------------------------------------------------#
# Define a weighting function (for example x**2)

def weight_function(x):
	return x**2

#----------------------------------------------------------------------------------------#
# Input: a population V of n weighted items

n = 100000

v_items = []
v_weight = []
for i in range(n):
	x = randint(1,100)
	v_items.append(x)
	v_weight.append(weight_function(x))

plt.hist(v_items, color="lightblue")

plt.grid()

plt.title('Weighted random sampling with a reservoir \n (initial population 1/2)')
plt.ylabel('Counts')
plt.xlim(0,100)

#plt.savefig("population.png", bbox_inches='tight')
plt.show()

#----------------------------------------------------------------------------------------#
# Weighted Random Sampling 

m = 1000 # m < n (sampling size)

reservoir = np.zeros((m,2))

for i in range(n):
	ui = uniform(0.0,1.0)
	if i < m:
		ki = ui**(1.0/v_weight[i])
		reservoir[i,0] = v_items[i]
		reservoir[i,1] = ki
	else:
		ki = ui**(1.0/v_weight[i])
		ki_min_idx = np.argmin(reservoir[:,1])
		ki_min = reservoir[ki_min_idx,1]
		if ki > ki_min:
			reservoir[ki_min_idx,0] = v_items[i]
			reservoir[ki_min_idx,1] = ki
			ki_min_idx = np.argmin(reservoir[:,1])
			ki_min = reservoir[ki_min_idx,1]
		
#----------------------------------------------------------------------------------------#
# Plot Results

plt.hist(reservoir[:,0], color="red")

plt.grid()

plt.title('Weighted random sampling with a reservoir \n (Sampling 2/2)')
plt.ylabel('Counts')
plt.xlim(0,100)

#plt.savefig("sampling.png", bbox_inches='tight')
plt.show()

Fortran 90 module to read MODIS MOD/MYD021 products

					
MODULE get_modis_myd021_data

implicit none

!----------------------------------------------------------------------------------------!
! Dynamic Allocation of Arrays

integer*2, dimension(:,:,:), allocatable :: EV_250_Aggr1km_RefSB !!!!! data type 23 => integer 16 bits unsigned !!!!!!
integer*2, dimension(:,:,:), allocatable :: EV_500_Aggr1km_RefSB !!!!! data type 23 => integer 16 bits unsigned !!!!!!
integer*2, dimension(:,:,:), allocatable :: EV_1KM_RefSB !!!!! data type 23 => integer 16 bits unsigned !!!!!!
integer*2, dimension(:,:,:), allocatable :: EV_1KM_Emissive !!!!! data type 23 => integer 16 bits unsigned !!!!!!

real(kind=4), dimension(2) :: EV_250_Aggr1km_RefSB_reflectance_offsets
real(kind=4), dimension(2) :: EV_250_Aggr1km_RefSB_reflectance_scales

real(kind=4), dimension(5) :: EV_500_Aggr1km_RefSB_reflectance_offsets
real(kind=4), dimension(5) :: EV_500_Aggr1km_RefSB_reflectance_scales

real(kind=4), dimension(15) :: EV_1KM_RefSB_reflectance_offsets
real(kind=4), dimension(15) :: EV_1KM_RefSB_reflectance_scales

real(kind=4), dimension(16) :: EV_1KM_Emissive_radiance_offsets
real(kind=4), dimension(16) :: EV_1KM_Emissive_radiance_scales

!----------------------------------------------------------------------------------------!
! User Inputs (user_selected_sds)

integer :: nb_user_selected_sds
parameter ( nb_user_selected_sds = 4 )

character*300 selected_sds( nb_user_selected_sds )      

data selected_sds / 'EV_250_Aggr1km_RefSB','EV_500_Aggr1km_RefSB', &
                  & 'EV_1KM_RefSB','EV_1KM_Emissive'/

!----------------------------------------------------------------------------------------!
! Subroutine read_cloudsat_2b_geoprof

CONTAINS

subroutine read_modis_myd021(modis_c6_myd021_file_name)

use read_write_hdf_file_variables

implicit none

integer :: i,j,k,l,m,n
integer :: AllocateStatus, DeAllocateStatus
integer :: ierr ! etat de reussite (si ==0) ou echec (si /= 0)
character*300 :: modis_c6_myd021_file_name

integer sfstart, sfselect, sfrdata, sfendacc, sfend , sffattr , sfrattr, var_id, attr_id, err_code
integer status, n_attrs
integer sffinfo , sfginfo

integer :: sfgainfo

integer :: attr_index, n_values

character*20 attr_name

sd_id = sfstart(modis_c6_myd021_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do index = 0, n_datasets - 1

	sds_id = sfselect(sd_id, index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )
    sds_name = name(1:len(name))

if ( ANY( selected_sds == sds_name ) ) then

	!do attr_index = 0, n_attrs - 1
   	! status = sfgainfo(sds_id, attr_index, attr_name, data_type, n_values)
    !	write(6,*) attr_name, data_type, n_values
	!end do

	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do
	
!----------------------------------------------------------------------------------------!
	
	if( sds_name == 'EV_250_Aggr1km_RefSB' )then
	
		allocate( EV_250_Aggr1km_RefSB(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus)
		IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus EV_250_Aggr1km_RefSB)***"
		status = sfrdata( sds_id, start, stride, edges, EV_250_Aggr1km_RefSB )
	
		attr_index = sffattr(sds_id, 'reflectance_offsets')
		status = sfrattr(sds_id, attr_index, EV_250_Aggr1km_RefSB_reflectance_offsets)

		attr_index = sffattr(sds_id, 'reflectance_scales')
		status = sfrattr(sds_id, attr_index, EV_250_Aggr1km_RefSB_reflectance_scales)	
	
	end if

!----------------------------------------------------------------------------------------!

	if( sds_name == 'EV_500_Aggr1km_RefSB' )then
	
		allocate( EV_500_Aggr1km_RefSB(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus)
		IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus EV_500_Aggr1km_RefSB)***"
		status = sfrdata( sds_id, start, stride, edges, EV_500_Aggr1km_RefSB )

		attr_index = sffattr(sds_id, 'reflectance_offsets')
		status = sfrattr(sds_id, attr_index, EV_500_Aggr1km_RefSB_reflectance_offsets)

		attr_index = sffattr(sds_id, 'reflectance_scales')
		status = sfrattr(sds_id, attr_index, EV_500_Aggr1km_RefSB_reflectance_scales)
	
	end if

!----------------------------------------------------------------------------------------!
	
	if( sds_name == 'EV_1KM_RefSB' )then
		
		allocate( EV_1KM_RefSB(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus)
		IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus EV_1KM_RefSB)***"
		status = sfrdata( sds_id, start, stride, edges, EV_1KM_RefSB )

		attr_index = sffattr(sds_id, 'reflectance_offsets')
		status = sfrattr(sds_id, attr_index, EV_1KM_RefSB_reflectance_offsets)

		attr_index = sffattr(sds_id, 'reflectance_scales')
		status = sfrattr(sds_id, attr_index, EV_1KM_RefSB_reflectance_scales)
	
	end if

!----------------------------------------------------------------------------------------!	
	
	if( sds_name == 'EV_1KM_Emissive' )then
	
		allocate( EV_1KM_Emissive(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus)
		IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus EV_1KM_Emissive)***"
		status = sfrdata( sds_id, start, stride, edges, EV_1KM_Emissive )

		attr_index = sffattr(sds_id, 'radiance_offsets')
		status = sfrattr(sds_id, attr_index, EV_1KM_Emissive_radiance_offsets)

		attr_index = sffattr(sds_id, 'radiance_scales')
		status = sfrattr(sds_id, attr_index, EV_1KM_Emissive_radiance_scales)
	
	end if	

!----------------------------------------------------------------------------------------!

end if

end do

!----------------------------------------------------------------------------------------!

end subroutine read_modis_myd021

!----------------------------------------------------------------------------------------!

END MODULE get_modis_myd021_data

Fortran 90 code to create MODIS CALIOP CLOUDSAT products samples

					
program create_modis_caliop_cloudsat_daily_samples

! Note:
!
! Compilation
!
! make -f makefile.modis_caliop_cloudsat_multilayer_clouds
! make clear -f makefile.modis_caliop_cloudsat_multilayer_clouds
!
! Execution:
! run_modis_caliop_cloudsat_multilayer_clouds
! /Users/marchant/Desktop/run_modis_caliop_cloudsat_multilayer_clouds
! /home/benjamin.hg.marchant/projects/evaluating_modis_multilayer_and_cloud_phase_algorithms_with_caliop_cloudsat/run_modis_caliop_cloudsat_multilayer_clouds
!
! python create_modis_caliop_cloudsat_monthly_samples.py 2008 07


use read_write_hdf_file_variables 
use get_modis_c6_myd03_data
use get_modis_c6_myd06_data
use get_caliop_v4_01kmCLay_data
use get_caliop_v4_05kmCLay_data
use get_cldclass_lidar_data
use get_modis_myd021_data

implicit none 

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,k,l,m,n
integer :: mod_file
integer :: year,month,day,hour,minute,julian_day
character*300 :: path_calmod
character*100 :: calmod_index,calmod_file_name
character*400 :: input_colocated_file_index,input_colocated_file_name
integer :: ierror, iazm, iszen, iview
integer :: string_len

character*4 :: hhmm
character*300 :: file_name_collocated_pixels, file_name_modis_linked
character*300 :: path_to_modis_caliop_collocated_pixels, path_to_modis_linked_files
character*300 :: path_to_modis_caliop_collocated_pixels_01, path_to_modis_linked_files_01

integer :: ierr ! successful (if == 0) ou failed (si /= 0)

character*50 :: machine_id = 'nasa_desktop' 
!character*50 :: machine_id = 'icare_server'

integer :: nb_colocated_data

logical :: found_colocated_data, missing_file

logical :: input_files_exist

integer, dimension(:,:), allocatable :: modis_index
integer, dimension(:), allocatable :: caliop_index, cloudsat_index

integer :: cldsat_idx

integer :: AllocateStatus, DeAllocateStatus

character*300 :: input_MYD021
character*300 :: input_MYD03
character*300 :: input_MYD06_c6
character*300 :: input_CAL_01km
character*300 :: input_CAL_05km
character*300 :: input_cloudsat_2b_geoprof
character*300 :: input_cldclass_lidar

character*16 :: file_common_path

character*11 :: output_file_prefix

integer :: s1,s2,s3,s4,s5,s6,s7,s8

integer, dimension(12) :: dpm_nyear, dpm_byear
data dpm_nyear /31,28,31,30,31,30,31,31,30,31,30,31/
data dpm_byear /31,29,31,30,31,30,31,31,30,31,30,31/

!----------------------------------------------------------------------------------------!
! New Variables

integer :: cloud_mask_flag
integer :: surface_flag
integer :: day_flag 
integer :: modis_multilayercloud_flag
integer :: modis_multilayercloud_ph_flag
integer :: modis_multilayercloud_ph_turned_off_flag
integer :: modis_cloud_phase_flag
integer :: caliop_cloud_phase_flag
integer :: caliop_multilayercloud_flag

integer :: tmp_i1
integer :: caliop_1km_nb_cloud_layers,caliop_5km_nb_cloud_layers,nb_cloud_layers
integer, dimension(10) :: caliop_1km_cloud_layer_phase, caliop_5km_cloud_layer_phase
real(pr), dimension(10) :: caliop_1km_cloud_layer_top, caliop_5km_cloud_layer_top
real(pr), dimension(10) :: caliop_1km_cloud_layer_base, caliop_5km_cloud_layer_base
integer, dimension(10) :: cal_Feature_Type_1km, cal_Feature_Type_5km
real(pr), dimension(10) :: caliop_1km_cloud_layer_optical_deth, caliop_5km_cloud_layer_optical_deth
integer, dimension(20) :: cloud_layer_phase
real(pr), dimension(20) :: cloud_layer_top
real(pr), dimension(20) :: cloud_layer_base
real(pr), dimension(20) :: cloud_layer_optical_depth
integer, dimension(20) :: cloud_layer_source
logical :: found_ice_cloud_layers,found_liq_cloud_layers
logical :: found_cloud_layers, found_aerosol_layers

real(pr) :: tmp1, tmp2

integer :: idx_d

real(pr) :: distance_max_found
integer :: idx_upper_layer_base
real(pr) :: tau_upper_layers
real(pr) :: tau_max_found
logical :: cst_cond
integer :: cloud_scene_type
integer :: cloud_phase_upper, cloud_phase_lower
logical :: multiple_cloud_phase_upper, multiple_cloud_phase_lower
integer, dimension(20) :: cloud_phase_profile
integer :: nb_cloud_phase_layer_found
logical, dimension(8) :: logical_cloud_scene_type

logical	upper_ice_cloud, upper_liq_cloud
logical	lower_ice_cloud, lower_liq_cloud
integer :: nb_upper_cloud_layer, nb_lower_cloud_layer
integer :: nb_ice_cloud_layer, nb_liq_cloud_layer

logical :: cloud_layer_ice = .False.
logical :: cloud_layer_liq = .False.

logical :: cloudsat_cloud_layer

integer :: cldclass_lidar_nb_cloud_layers

integer, dimension(10) :: cldclass_lidar_cloud_layer_phase
real(pr), dimension(10) :: cldclass_lidar_cloud_layer_top
real(pr), dimension(10) :: cldclass_lidar_cloud_layer_base
real(pr), dimension(10) :: cldclass_lidar_cloud_layer_optical_depth

integer, dimension(2) :: shape_caliop_v4_05kmCLay_Number_Layers_Found

integer :: s0_size, s1_size, s2_size, s3_size, s4_size, s5_size
integer :: s0_idx, s1_idx, s2_idx, s3_idx, s4_idx, s5_idx
integer :: s0_tot, s1_tot, s2_tot, s3_tot, s4_tot, s5_tot

real(pr), dimension(:,:), allocatable :: sampling_01, sampling_02, sampling_03 
real(pr), dimension(:,:), allocatable :: sampling_04, sampling_05, sampling_00

integer, dimension(2) :: array_shape

character(len=4) :: arg

logical :: file_exists

integer, dimension(:,:), allocatable :: gridded_data

integer, dimension(2) :: gridded_data_shape

real(pr) :: grid_resol

integer :: grid_i_max, grid_j_max

integer :: grid_i, grid_j, grid_bin

character*24 :: output_filename

integer :: modis_multilayer_cloud_c6_flag
integer :: modis_multilayer_cloud_c6_ph_turned_off_flag

!----------------------------------------------------------------------------------------!
! Settings

!year = 2008
!month = 07

call getarg(1, arg) 
read( arg, '(i4)' ) year

call getarg(2, arg) 
read( arg, '(i2)' ) month

call getarg(3, arg) 
read( arg, '(i2)' ) day

call getarg(4, arg) 
read( arg, '(i2)' ) hour

call getarg(5, arg) 
read( arg, '(i2)' ) minute

!write(6,*) '---- test ----',year, month, day

!do i = 1, iargc() 
!  call getarg(i, arg) 
!  write (6,*) arg 
!  read( arg, '(i4)' ) j ! Convert character to INTEGER
!  write (6,*) j 
!end do

output_file_prefix = '    _  _  _'

write(output_file_prefix(1:4),'(i4.4)') year
write(output_file_prefix(6:7),'(i2.2)') month
write(output_file_prefix(9:10),'(i2.2)') day

output_filename = 'gridded_data_    _  .txt'

write(output_filename(14:17),'(i4.4)') year
write(output_filename(19:20),'(i2.2)') month

julian_day = 0
do i = 1, month - 1
julian_day = julian_day + dpm_byear(i)
end do
julian_day = julian_day + day

if( trim(machine_id) == 'nasa_desktop' )then
  file_name_collocated_pixels = 'modis_caliop_cloudsat_collocated_indexes_'
  file_name_modis_linked = 'modis_linked_files_'
  path_to_modis_caliop_collocated_pixels = '/Users/marchant/Desktop/modis_caliop_cloudsat_collocated_pixel_files/'
  path_to_modis_linked_files = '/Users/marchant/Desktop/modis_based_match_files/'
else 
  if( trim(machine_id) == 'icare_server' )then
    file_name_collocated_pixels = 'modis_caliop_cloudsat_collocated_indexes_'
    file_name_modis_linked = 'modis_linked_files_'
    path_to_modis_caliop_collocated_pixels = '/home/benjamin.hg.marchant/data/modis_caliop_cloudsat_collocated_pixels_files/'
    path_to_modis_linked_files = '/home/benjamin.hg.marchant/data/modis_based_match_files/'
 else
   STOP
  end if
end if

file_common_path = '    /    _  _  /'

write(file_common_path(1:4),'(i4.4)') year
write(file_common_path(6:9),'(i4.4)') year
write(file_common_path(11:12),'(i2.2)') month

s0_size = 5000
s1_size = 5000
s2_size = 5000
s3_size = 5000
s4_size = 5000
s5_size = 5000

allocate( sampling_00(s0_size,9) )
allocate( sampling_01(s1_size,69) )
allocate( sampling_02(s2_size,69) )
allocate( sampling_03(s3_size,69) )
allocate( sampling_04(s4_size,114) )
allocate( sampling_05(s5_size,115) )

grid_resol = 1.0

grid_i_max = int( 180.0 / grid_resol )
grid_j_max = int( 360.0 / grid_resol ) 

allocate( gridded_data( grid_i_max * grid_j_max , 5 ) )

gridded_data_shape = shape(gridded_data)

gridded_data = 0

!----------------------------------------------------------------------------------------!
! Get previous Data

INQUIRE(file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_00.txt', EXIST=file_exists)

if(file_exists)then

	array_shape = shape(sampling_00)
	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_00.txt')
	read(1,*) s0_tot
	do i = 1, s0_size
	read(1,*) (sampling_00(i,j),j=1,array_shape(2))
	end do
	close(1)

	array_shape = shape(sampling_01)
	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_01.txt')
	read(1,*) s1_tot
	do i = 1, s1_size
	read(1,*) (sampling_01(i,j),j=1,array_shape(2))
	end do
	close(1)

	array_shape = shape(sampling_02)
	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_02.txt')
	read(1,*) s2_tot
	do i = 1, s2_size
	read(1,*) (sampling_02(i,j),j=1,array_shape(2))
	end do
	close(1)
	
	array_shape = shape(sampling_03)
	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_03.txt')
	read(1,*) s3_tot
	do i = 1, s3_size
	read(1,*) (sampling_03(i,j),j=1,array_shape(2))
	end do
	close(1)
	
	array_shape = shape(sampling_04)
	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_04.txt')
	read(1,*) s4_tot
	do i = 1, s4_size
	read(1,*) (sampling_04(i,j),j=1,array_shape(2))
	end do
	close(1)
	
	array_shape = shape(sampling_05)
	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_05.txt')
	read(1,*) s5_tot
	do i = 1, s5_size
	read(1,*) (sampling_05(i,j),j=1,array_shape(2))
	end do
	close(1)				

else

	write(6,*) "no previous data 1"

	sampling_00 = 0.0
	sampling_01 = 0.0
	sampling_02 = 0.0
	sampling_03 = 0.0
	sampling_04 = 0.0
	sampling_05 = 0.0

	s0_tot = 0
	s1_tot = 0
	s2_tot = 0
	s3_tot = 0
	s4_tot = 0
	s5_tot = 0	

end if

!write(6,*) s0_tot,s1_tot,s2_tot,s3_tot,s4_tot,s5_tot

INQUIRE(file="./outputs/modis_caliop_cloudsat_multilayer_clouds/"//output_filename, EXIST=file_exists)

if(file_exists)then

	open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_filename)
	read(1,*) 
	do i = 1, grid_i_max * grid_j_max
		read(1,*) (gridded_data(i,j),j=1,gridded_data_shape(2))
	end do
	close(1)

else

	write(6,*) "no previous data 2"

end if

!----------------------------------------------------------------------------------------!
! Iterate over day

!do day = 1, 1! 31

write(file_common_path(14:15),'(i2.2)') day

string_len = len( trim(path_to_modis_caliop_collocated_pixels) )
path_to_modis_caliop_collocated_pixels_01 = path_to_modis_caliop_collocated_pixels(1:string_len)//file_common_path

string_len = len( trim(path_to_modis_linked_files) )
path_to_modis_linked_files_01 = path_to_modis_linked_files(1:string_len)//file_common_path

!hour = 0
!minute = 0

!do mod_file = 1, 288

hhmm = '    '

write(hhmm(1:2),'(i2.2)') hour
write(hhmm(3:4),'(i2.2)') minute

string_len = len( trim(file_name_collocated_pixels) )
input_colocated_file_index = file_name_collocated_pixels(1:string_len)//hhmm

string_len = len( trim(file_name_modis_linked) )
input_colocated_file_name = file_name_modis_linked(1:string_len)//hhmm

string_len = len( trim(path_to_modis_caliop_collocated_pixels_01) )
input_colocated_file_index = path_to_modis_caliop_collocated_pixels_01(1:string_len)//input_colocated_file_index

string_len = len( trim(path_to_modis_linked_files_01) )
input_colocated_file_name = path_to_modis_linked_files_01(1:string_len)//input_colocated_file_name

found_colocated_data = .TRUE.
open(1,file=input_colocated_file_index,form="formatted",iostat=ierr,status="old")
if (ierr/=0) then
!write(*,*) "Impossible d ouvrir le fichier ",trim(input_colocated_file_index)
!stop
found_colocated_data = .FALSE.
endif
nb_colocated_data = 0
do while (ierr==0)
read(1,*,iostat=ierr)
if (ierr==0) nb_colocated_data = nb_colocated_data + 1
enddo
close(1)

missing_file = .FALSE.
open(1,file=input_colocated_file_name,form="formatted",iostat=ierr,status="old")
if (ierr/=0) then
!write(*,*) "Impossible d ouvrir le fichier ",trim(input_colocated_file_name)
!stop
missing_file = .TRUE.
endif
close(1)

if( found_colocated_data .and. .not. missing_file )then

nb_colocated_data = nb_colocated_data - 3

if( nb_colocated_data > 0 )then

allocate( modis_index(nb_colocated_data,2), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_index***"
allocate( caliop_index(nb_colocated_data), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory caliop_index***"
allocate( cloudsat_index(nb_colocated_data), STAT = AllocateStatus )
IF (AllocateStatus /= 0) STOP "*** Not enough memory cloudsat_index***"

open(1,file=input_colocated_file_index, status="old")
read(1,'(a)') 
read(1,'(a)') 
read(1,'(a)') 
do i = 1, nb_colocated_data 
read(1,*) modis_index(i,1), modis_index(i,2), caliop_index(i), cloudsat_index(i)
end do
close(1)

input_files_exist = .TRUE.

if( trim(machine_id) == 'nasa_desktop' )then
  open(1,file=input_colocated_file_name, status="old")
  read(1,'(a)') input_MYD021
  read(1,'(a)') input_MYD03
  read(1,'(a)') 
  read(1,'(a)') 
  read(1,'(a)') input_MYD06_c6
  read(1,'(a)') input_CAL_01km
  read(1,'(a)') input_CAL_05km
  read(1,'(a)') input_cloudsat_2b_geoprof
  read(1,'(a)') input_cldclass_lidar
  close(1)
  if( input_MYD03 == 'none' .OR. input_MYD06_c6 == 'none' )then
  	input_files_exist = .FALSE.
  end if
  if( input_CAL_01km == 'none' .OR. input_CAL_05km == 'none' )then
  	input_files_exist = .FALSE.
  end if
  if( input_cldclass_lidar == 'none' )then
  	input_files_exist = .FALSE.
  end if    
  input_MYD021 = '/Volumes/Data/MYD021km/'//file_common_path//input_MYD021
  input_MYD03 = '/Volumes/Data/Myd03/'//file_common_path//input_MYD03
  input_MYD06_c6 = '/Volumes/DATA_03/modis/Myd06_c6/'//file_common_path//input_MYD06_c6
  input_CAL_01km = '/Volumes/Data/Caliop/v4/01km_CLay/'//file_common_path//input_CAL_01km
  input_CAL_05km = '/Volumes/Data/Caliop/v4/05km_CLay/'//file_common_path//input_CAL_05km
  input_cloudsat_2b_geoprof = '/Volumes/Data/Cloudsat/'//file_common_path//input_cloudsat_2b_geoprof
  input_cldclass_lidar = '/Volumes/Data/2B-CLDCLASS-LIDAR/'//file_common_path//input_cldclass_lidar
else 
  if( trim(machine_id) == 'icare_server' )then
	open(1,file=input_colocated_file_name, status="old")
	read(1,'(a)') input_MYD021
	read(1,'(a)') input_MYD03
	read(1,'(a)') 
	read(1,'(a)') 
	read(1,'(a)') input_MYD06_c6
	read(1,'(a)') input_CAL_01km
	read(1,'(a)') input_CAL_05km
	read(1,'(a)') input_cloudsat_2b_geoprof
	read(1,'(a)') input_cldclass_lidar
	close(1)
	if( input_MYD03 == 'none' .OR. input_MYD06_c6 == 'none' )then
		input_files_exist = .FALSE.
	end if
	if( input_CAL_01km == 'none' .OR. input_CAL_05km == 'none' )then
		input_files_exist = .FALSE.
	end if
	if( input_cldclass_lidar == 'none' )then
		input_files_exist = .FALSE.
	end if     
    input_MYD021 = '/DATA/LIENS/MODIS/MYD021KM.006/'//file_common_path//input_MYD021
    input_MYD03 = '/DATA/LIENS/MODIS/MYD03/'//file_common_path//input_MYD03
    input_MYD06_c6 = '/DATA/LIENS/MODIS/MYD06_L2.006/'//file_common_path//input_MYD06_c6
    input_CAL_01km = '/DATA/LIENS/CALIOP/01kmCLay.v4.10/'//file_common_path//input_CAL_01km
    input_CAL_05km = '/DATA/LIENS/CALIOP/05kmCLay.v4.10/'//file_common_path//input_CAL_05km
    input_cloudsat_2b_geoprof = '/DATA/LIENS/CLOUDSAT/2B-GEOPROF/'//file_common_path//input_cloudsat_2b_geoprof
    input_cldclass_lidar = '/DATA/LIENS/CLOUDSAT/2B-CLDCLASS-LIDAR.v04.02/'//file_common_path//input_cldclass_lidar
 else
   STOP
  end if
end if

if(input_files_exist)then

CALL read_modis_myd021(input_MYD021)
CALL read_modis_c6_myd03(input_MYD03)
CALL read_modis_c6_myd06(input_MYD06_c6)
CALL read_caliop_v4_01kmCLay(input_CAL_01km)
CALL read_caliop_v4_05kmCLay(input_CAL_05km)
CALL read_cldclass_lidar(input_cldclass_lidar)

!----------------------------------------------------------------------------------------!
! MODIS CALIOP colocated data iteration

do i = 1, nb_colocated_data 

j = caliop_index(i)
k = modis_index(i,1)
l = modis_index(i,2)
cldsat_idx = cloudsat_index(i)

if( j > 0 .and. k > 0 .and. l > 0 .and. cldsat_idx > 0 ) then 

!write(6,*) modis_c6_myd03_latitude(k,l), modis_c6_myd03_longitude(k,l)
!write(6,*) cloudsat_latitude(cldsat_idx), cloudsat_longitude(cldsat_idx)
!write(6,*) '---------'

if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1) == 1 .and. & ! Cloud Mask determined 
&   ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),1,2) < 2 )then  ! Cloudy or Probably Cloudy Pixel  
cloud_mask_flag = 1
else
cloud_mask_flag = 2
end if

day_flag = 0
if( modis_c6_myd03_solarzenith(k,l) * 0.01 < 81.36  ) day_flag = 1 ! MYD06 definition of the day 

modis_cloud_phase_flag = 0
if(modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l) == 2) modis_cloud_phase_flag = 1
if(modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l) == 3) modis_cloud_phase_flag = 2

if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 0 ) surface_flag = 0 ! Ocean
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 3 ) surface_flag = 3 ! Land
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 2 ) surface_flag = 2 ! Desert
if( ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),6,2) == 1 ) surface_flag = 1 ! Coastal

shape_caliop_v4_05kmCLay_Number_Layers_Found = shape(caliop_v4_05kmCLay_Number_Layers_Found)

modis_multilayer_cloud_c6_flag = 1
if( modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l) > 1.0 ) modis_multilayer_cloud_c6_flag = 2

tmp1 = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
if( ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1) == 1 ) tmp1 = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l) - 3.0

modis_multilayer_cloud_c6_ph_turned_off_flag = 1
if( tmp1 > 1.0 ) modis_multilayer_cloud_c6_ph_turned_off_flag = 2

!----------------------------------------------------------------------------------------!
! Get CALIOP 1km & 5km cloud layers

caliop_1km_cloud_layer_phase = 0
caliop_5km_cloud_layer_phase = 0

cal_Feature_Type_1km = 0 
cal_Feature_Type_5km = 0 

caliop_1km_cloud_layer_top = 0
caliop_5km_cloud_layer_top = 0

caliop_1km_cloud_layer_base = 0.0
caliop_5km_cloud_layer_base = 0.0

caliop_1km_cloud_layer_optical_deth = 0.0
caliop_5km_cloud_layer_optical_deth = 0.0

caliop_1km_nb_cloud_layers = 0
if( caliop_v4_01kmCLay_Number_Layers_Found(1,j) >= 1 )then
	do n = 1, caliop_v4_01kmCLay_Number_Layers_Found(1,j)
	if( ibits(caliop_v4_01kmCLay_Feature_Classification_Flags(n,j),0,3) == 2 )then
		caliop_1km_nb_cloud_layers = caliop_1km_nb_cloud_layers + 1	
		caliop_1km_cloud_layer_phase(caliop_1km_nb_cloud_layers) = ibits(caliop_v4_01kmCLay_Feature_Classification_Flags(n,j),5,2)
		cal_Feature_Type_1km(caliop_1km_nb_cloud_layers) = ibits(caliop_v4_01kmCLay_Feature_Classification_Flags(n,j),0,3)
		caliop_1km_cloud_layer_top(caliop_1km_nb_cloud_layers) = caliop_v4_01kmCLay_Layer_Top_Altitude(n,j)
		caliop_1km_cloud_layer_base(caliop_1km_nb_cloud_layers) = caliop_v4_01kmCLay_Layer_Base_Altitude(n,j)
		caliop_1km_cloud_layer_optical_deth(caliop_1km_nb_cloud_layers) = -9999.0
		!write(6,*) caliop_v4_01kmCLay_Layer_Top_Altitude(n,j), caliop_v4_01kmCLay_Layer_Base_Altitude(n,j)
	end if
	end do
end if

m = int(j/5.0) + 1
caliop_5km_nb_cloud_layers = 0
if( m <= shape_caliop_v4_05kmCLay_Number_Layers_Found(2) )then
if( caliop_v4_05kmCLay_Number_Layers_Found(1,m) >= 1 )then
	do n = 1, caliop_v4_05kmCLay_Number_Layers_Found(1,m)
	if( ibits(caliop_v4_05kmCLay_Feature_Classification_Flags(n,m),0,3) == 2 )then
		caliop_5km_nb_cloud_layers = caliop_5km_nb_cloud_layers + 1	
		caliop_5km_cloud_layer_phase(caliop_5km_nb_cloud_layers) = ibits(caliop_v4_05kmCLay_Feature_Classification_Flags(n,m),5,2)
		cal_Feature_Type_5km(caliop_5km_nb_cloud_layers) = ibits(caliop_v4_05kmCLay_Feature_Classification_Flags(n,m),0,3)
		caliop_5km_cloud_layer_top(caliop_5km_nb_cloud_layers) = caliop_v4_05kmCLay_Layer_Top_Altitude(n,m)
		caliop_5km_cloud_layer_base(caliop_5km_nb_cloud_layers) = caliop_v4_05kmCLay_Layer_Base_Altitude(n,m)
		!write(6,*) 'caliop_v4_05kmCLay_Feature_Optical_Depth_532', caliop_v4_05kmCLay_Feature_Optical_Depth_532(n,m)
		caliop_5km_cloud_layer_optical_deth(caliop_5km_nb_cloud_layers) = caliop_v4_05kmCLay_Feature_Optical_Depth_532(n,m)
	end if
	end do
end if
end if

!----------------------------------------------------------------------------------------!
! Get 2b cldclass lidar cloud layers 

cldclass_lidar_cloud_layer_phase = 0.0
cldclass_lidar_cloud_layer_top = 0.0
cldclass_lidar_cloud_layer_base = 0.0
cldclass_lidar_cloud_layer_optical_depth = 0.0

cldclass_lidar_nb_cloud_layers = 0

do n = 1, 10
	if (cldclass_lidar_CloudPhase(n,cldsat_idx) /= 0)then
		cldclass_lidar_nb_cloud_layers = cldclass_lidar_nb_cloud_layers + 1
		cldclass_lidar_cloud_layer_phase(cldclass_lidar_nb_cloud_layers) = cldclass_lidar_CloudPhase(n,cldsat_idx)
		cldclass_lidar_cloud_layer_top(cldclass_lidar_nb_cloud_layers) = cldclass_lidar_CloudLayerTop(n,cldsat_idx)
		cldclass_lidar_cloud_layer_base(cldclass_lidar_nb_cloud_layers) = cldclass_lidar_CloudLayerBase(n,cldsat_idx)
	end if
end do

!----------------------------------------------------------------------------------------!
! Concatenate 2b cldclass lidar, caliop 1km & 5km cloud layer products (max 10 layers)

nb_cloud_layers = 0

cloud_layer_phase = 0
cloud_layer_top = 0.0
cloud_layer_base = 0.0
cloud_layer_optical_depth = -9999.0
cloud_layer_source = 0

nb_cloud_layers = 0
do n = 1, caliop_1km_nb_cloud_layers
	nb_cloud_layers = nb_cloud_layers + 1
	cloud_layer_phase(nb_cloud_layers) = caliop_1km_cloud_layer_phase(n)
	cloud_layer_top(nb_cloud_layers) = caliop_1km_cloud_layer_top(n)
	cloud_layer_base(nb_cloud_layers) = caliop_1km_cloud_layer_base(n)
	cloud_layer_source(nb_cloud_layers) = 1
end do
do n = 1, caliop_5km_nb_cloud_layers
	nb_cloud_layers = nb_cloud_layers + 1
	cloud_layer_phase(nb_cloud_layers) = caliop_5km_cloud_layer_phase(n)
	cloud_layer_top(nb_cloud_layers) = caliop_5km_cloud_layer_top(n)
	cloud_layer_base(nb_cloud_layers) = caliop_5km_cloud_layer_base(n)
	cloud_layer_optical_depth(nb_cloud_layers) = caliop_5km_cloud_layer_optical_deth(n)
	cloud_layer_source(nb_cloud_layers) = 2
end do

do n = 1, cldclass_lidar_nb_cloud_layers
	nb_cloud_layers = nb_cloud_layers + 1
	cloud_layer_phase(nb_cloud_layers) = cldclass_lidar_cloud_layer_phase(n)
	cloud_layer_top(nb_cloud_layers) = cldclass_lidar_cloud_layer_top(n)
	cloud_layer_base(nb_cloud_layers) = cldclass_lidar_cloud_layer_base(n)
	cloud_layer_optical_depth(nb_cloud_layers) = cldclass_lidar_cloud_layer_optical_depth(n)
	cloud_layer_source(nb_cloud_layers) = 3
end do

!----------------------------------------------------------------------------------------!
! Sort cloud layers by cloud top altitude

do n = 1, nb_cloud_layers - 1
	do m = n + 1, nb_cloud_layers
	if( cloud_layer_top(m) > cloud_layer_top(n) )then

	tmp1 = cloud_layer_top(n)
	tmp2 = cloud_layer_base(n)
	cloud_layer_top(n) = cloud_layer_top(m)
	cloud_layer_base(n) = cloud_layer_base(m)
	cloud_layer_top(m) = tmp1
	cloud_layer_base(m) = tmp2
		
	!----- cot_5km -----!

	tmp1 = cloud_layer_optical_depth(n)
	cloud_layer_optical_depth(n) = cloud_layer_optical_depth(m)
	cloud_layer_optical_depth(m) = tmp1

	!----- cloud phase -----!
	
	tmp1 = cloud_layer_phase(n)
	cloud_layer_phase(n) = cloud_layer_phase(m)
	cloud_layer_phase(m) = tmp1

	!----- source -----!

	tmp1 = cloud_layer_source(n)
	cloud_layer_source(n) = cloud_layer_source(m)
	cloud_layer_source(m) = tmp1
	
	end if
	end do 
end do

!----------------------------------------------------------------------------------------!
! Random Sample 0

if( day_flag == 1 )then

s0_tot = s0_tot + 1

if( s0_tot <= s0_size )then

s0_idx = s0_tot
	
sampling_00(s0_idx,1) = modis_c6_myd03_latitude(k,l)
sampling_00(s0_idx,2) = modis_c6_myd03_longitude(k,l)
sampling_00(s0_idx,3) = ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1)
sampling_00(s0_idx,4) = ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),1,2)
sampling_00(s0_idx,5) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_00(s0_idx,6) = nb_cloud_layers
sampling_00(s0_idx,7) = caliop_1km_nb_cloud_layers
sampling_00(s0_idx,8) = caliop_5km_nb_cloud_layers
sampling_00(s0_idx,9) = cldclass_lidar_nb_cloud_layers

else

s0_idx = int( rand(0) * ( s0_tot - 1 ) ) + 1

if( s0_idx <= s0_size )then

sampling_00(s0_idx,1) = modis_c6_myd03_latitude(k,l)
sampling_00(s0_idx,2) = modis_c6_myd03_longitude(k,l)
sampling_00(s0_idx,3) = ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1)
sampling_00(s0_idx,4) = ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),1,2)
sampling_00(s0_idx,5) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_00(s0_idx,6) = nb_cloud_layers
sampling_00(s0_idx,7) = caliop_1km_nb_cloud_layers
sampling_00(s0_idx,8) = caliop_5km_nb_cloud_layers
sampling_00(s0_idx,9) = cldclass_lidar_nb_cloud_layers

end if
end if

end if 

!----------------------------------------------------------------------------------------!
! Random Sample 1

if( day_flag == 1 )then
if( cloud_mask_flag == 1 )then
if( nb_cloud_layers > 0 )then

s1_tot = s1_tot + 1

if( s1_tot <= s1_size )then

s1_idx = s1_tot
	
sampling_01(s1_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_01(s1_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_01(s1_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_01(s1_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_01(s1_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_01(s1_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_01(s1_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_01(s1_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_01(s1_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_01(s1_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_01(s1_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_01(s1_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_01(s1_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_01(s1_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_01(s1_idx,15) = surface_flag	
sampling_01(s1_idx,16) = nb_cloud_layers
sampling_01(s1_idx,17) = caliop_1km_nb_cloud_layers
sampling_01(s1_idx,18) = caliop_5km_nb_cloud_layers
sampling_01(s1_idx,19) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 10
sampling_01(s1_idx,m) = cloud_layer_top(n) 
sampling_01(s1_idx,m+1) = cloud_layer_base(n)
sampling_01(s1_idx,m+2) = cloud_layer_optical_depth(n)
sampling_01(s1_idx,m+3) = cloud_layer_phase(n)
sampling_01(s1_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do
	
else

s1_idx = int( rand(0) * ( s1_tot - 1 ) ) + 1

if( s1_idx <= s1_size )then
	
sampling_01(s1_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_01(s1_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_01(s1_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_01(s1_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_01(s1_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_01(s1_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_01(s1_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_01(s1_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_01(s1_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_01(s1_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_01(s1_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_01(s1_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_01(s1_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_01(s1_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_01(s1_idx,15) = surface_flag	
sampling_01(s1_idx,16) = nb_cloud_layers
sampling_01(s1_idx,17) = caliop_1km_nb_cloud_layers
sampling_01(s1_idx,18) = caliop_5km_nb_cloud_layers
sampling_01(s1_idx,19) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 10
sampling_01(s1_idx,m) = cloud_layer_top(n) 
sampling_01(s1_idx,m+1) = cloud_layer_base(n)
sampling_01(s1_idx,m+2) = cloud_layer_optical_depth(n)
sampling_01(s1_idx,m+3) = cloud_layer_phase(n)
sampling_01(s1_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do

end if
end if

end if
end if
end if 

!----------------------------------------------------------------------------------------!
! Random Sample 2

if( day_flag == 1 )then
if( cloud_mask_flag == 1 )then
if( nb_cloud_layers == 2 )then

s2_tot = s2_tot + 1

if( s2_tot <= s2_size )then

s2_idx = s2_tot
	
sampling_02(s2_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_02(s2_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_02(s2_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_02(s2_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_02(s2_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_02(s2_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_02(s2_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_02(s2_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_02(s2_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_02(s2_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_02(s2_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_02(s2_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_02(s2_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_02(s2_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_02(s2_idx,15) = surface_flag	
sampling_02(s2_idx,16) = nb_cloud_layers
sampling_02(s2_idx,17) = caliop_1km_nb_cloud_layers
sampling_02(s2_idx,18) = caliop_5km_nb_cloud_layers
sampling_02(s2_idx,19) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 10
sampling_02(s2_idx,m) = cloud_layer_top(n) 
sampling_02(s2_idx,m+1) = cloud_layer_base(n)
sampling_02(s2_idx,m+2) = cloud_layer_optical_depth(n)
sampling_02(s2_idx,m+3) = cloud_layer_phase(n)
sampling_02(s2_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do
	
else

s2_idx = int( rand(0) * ( s2_tot - 1 ) ) + 1

if( s2_idx <= s2_size )then
	
sampling_02(s2_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_02(s2_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_02(s2_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_02(s2_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_02(s2_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_02(s2_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_02(s2_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_02(s2_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_02(s2_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_02(s2_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_02(s2_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_02(s2_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_02(s2_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_02(s2_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_02(s2_idx,15) = surface_flag	
sampling_02(s2_idx,16) = nb_cloud_layers
sampling_02(s2_idx,17) = caliop_1km_nb_cloud_layers
sampling_02(s2_idx,18) = caliop_5km_nb_cloud_layers
sampling_02(s2_idx,19) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 10
sampling_02(s2_idx,m) = cloud_layer_top(n) 
sampling_02(s2_idx,m+1) = cloud_layer_base(n)
sampling_02(s2_idx,m+2) = cloud_layer_optical_depth(n)
sampling_02(s2_idx,m+3) = cloud_layer_phase(n)
sampling_02(s2_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do

end if
end if

end if 
end if 
end if 

!----------------------------------------------------------------------------------------!
! Random Sample 3

if( day_flag == 1 )then
if( cloud_mask_flag == 1 )then
if( nb_cloud_layers == 2 )then

cloud_layer_ice = .False.
cloud_layer_liq = .False.

do n = 1, 2
if( cloud_layer_phase(n) == 1 .and. cloud_layer_source(n) /= 3) cloud_layer_ice = .True.
if( cloud_layer_phase(n) == 2 .and. cloud_layer_source(n) /= 3) cloud_layer_liq = .True.
if( cloud_layer_phase(n) == 1 .and. cloud_layer_source(n) == 3) cloud_layer_ice = .True.
if( cloud_layer_phase(n) == 3 .and. cloud_layer_source(n) == 3) cloud_layer_liq = .True.
end do

if( cloud_layer_ice .and. cloud_layer_liq )then

s3_tot = s3_tot + 1

if( s3_tot <= s3_size )then

s3_idx = s3_tot
	
sampling_03(s3_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_03(s3_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_03(s3_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_03(s3_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_03(s3_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_03(s3_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_03(s3_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_03(s3_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_03(s3_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_03(s3_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_03(s3_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_03(s3_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_03(s3_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_03(s3_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_03(s3_idx,15) = surface_flag	
sampling_03(s3_idx,16) = nb_cloud_layers
sampling_03(s3_idx,17) = caliop_1km_nb_cloud_layers
sampling_03(s3_idx,18) = caliop_5km_nb_cloud_layers
sampling_03(s3_idx,19) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 10
sampling_03(s3_idx,m) = cloud_layer_top(n) 
sampling_03(s3_idx,m+1) = cloud_layer_base(n)
sampling_03(s3_idx,m+2) = cloud_layer_optical_depth(n)
sampling_03(s3_idx,m+3) = cloud_layer_phase(n)
sampling_03(s3_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do
	
else

s3_idx = int( rand(0) * ( s3_tot - 1 ) ) + 1

if( s3_idx <= s3_size )then
	
sampling_03(s3_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_03(s3_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_03(s3_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_03(s3_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_03(s3_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_03(s3_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_03(s3_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_03(s3_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_03(s3_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_03(s3_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_03(s3_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_03(s3_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_03(s3_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_03(s3_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_03(s3_idx,15) = surface_flag	
sampling_03(s3_idx,16) = nb_cloud_layers
sampling_03(s3_idx,17) = caliop_1km_nb_cloud_layers
sampling_03(s3_idx,18) = caliop_5km_nb_cloud_layers
sampling_03(s3_idx,19) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 10
sampling_03(s3_idx,m) = cloud_layer_top(n) 
sampling_03(s3_idx,m+1) = cloud_layer_base(n)
sampling_03(s3_idx,m+2) = cloud_layer_optical_depth(n)
sampling_03(s3_idx,m+3) = cloud_layer_phase(n)
sampling_03(s3_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do

end if
end if

end if 
end if 
end if 
end if 

!----------------------------------------------------------------------------------------!
! Random Sample 4

if( day_flag == 1 )then
if( cloud_mask_flag == 1 )then
if( nb_cloud_layers > 0 )then

s4_tot = s4_tot + 1

if( s4_tot <= s4_size )then

s4_idx = s4_tot
	
sampling_04(s4_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_04(s4_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_04(s4_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_04(s4_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_04(s4_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_04(s4_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_04(s4_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_04(s4_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_04(s4_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_04(s4_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_04(s4_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_04(s4_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_04(s4_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_04(s4_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_04(s4_idx,15) = surface_flag	
sampling_04(s4_idx,16) = nb_cloud_layers
sampling_04(s4_idx,17) = caliop_1km_nb_cloud_layers
sampling_04(s4_idx,18) = caliop_5km_nb_cloud_layers
sampling_04(s4_idx,19) = cldclass_lidar_nb_cloud_layers

m = 19
do n = 1, 2
sampling_04(s4_idx,m+n) = EV_250_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_250_Aggr1km_RefSB(k,l,n) - EV_250_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 21
do n = 1, 5
sampling_04(s4_idx,m+n) = EV_500_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_500_Aggr1km_RefSB(k,l,n) - EV_500_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 26
do n = 1, 15
sampling_04(s4_idx,m+n) = EV_1KM_RefSB_reflectance_scales(n) * &
& ( EV_1KM_RefSB(k,l,n) - EV_1KM_RefSB_reflectance_offsets(n) )
end do

m = 41
do n = 1, 16
sampling_04(s4_idx,m+n) = EV_1KM_Emissive_radiance_scales(n) * &
& ( EV_1KM_Emissive(k,l,n) - EV_1KM_Emissive_radiance_offsets(n) )
end do

m = 57
do n = 1, 6
sampling_04(s4_idx,m+n) = modis_c6_myd06_Atm_Corr_Refl(n,k,l)
end do

m = 64
do n = 1, 10
sampling_04(s4_idx,m) = cloud_layer_top(n) 
sampling_04(s4_idx,m+1) = cloud_layer_base(n)
sampling_04(s4_idx,m+2) = cloud_layer_optical_depth(n)
sampling_04(s4_idx,m+3) = cloud_layer_phase(n)
sampling_04(s4_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do
	
else

s4_idx = int( rand(0) * ( s4_tot - 1 ) ) + 1

if( s4_idx <= s4_size )then
	
sampling_04(s4_idx,1) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_04(s4_idx,2) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_04(s4_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_04(s4_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_04(s4_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_04(s4_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_04(s4_idx,7) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_04(s4_idx,8) = modis_c6_myd03_latitude(k,l)
sampling_04(s4_idx,9) = modis_c6_myd03_longitude(k,l)
sampling_04(s4_idx,10) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_04(s4_idx,11) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_04(s4_idx,12) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_04(s4_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_04(s4_idx,14) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_04(s4_idx,15) = surface_flag	
sampling_04(s4_idx,16) = nb_cloud_layers
sampling_04(s4_idx,17) = caliop_1km_nb_cloud_layers
sampling_04(s4_idx,18) = caliop_5km_nb_cloud_layers
sampling_04(s4_idx,19) = cldclass_lidar_nb_cloud_layers

m = 19
do n = 1, 2
sampling_04(s4_idx,m+n) = EV_250_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_250_Aggr1km_RefSB(k,l,n) - EV_250_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 21
do n = 1, 5
sampling_04(s4_idx,m+n) = EV_500_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_500_Aggr1km_RefSB(k,l,n) - EV_500_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 26
do n = 1, 15
sampling_04(s4_idx,m+n) = EV_1KM_RefSB_reflectance_scales(n) * &
& ( EV_1KM_RefSB(k,l,n) - EV_1KM_RefSB_reflectance_offsets(n) )
end do

m = 41
do n = 1, 16
sampling_04(s4_idx,m+n) = EV_1KM_Emissive_radiance_scales(n) * &
& ( EV_1KM_Emissive(k,l,n) - EV_1KM_Emissive_radiance_offsets(n) )
end do

m = 57
do n = 1, 6
sampling_04(s4_idx,m+n) = modis_c6_myd06_Atm_Corr_Refl(n,k,l)
end do

m = 64
do n = 1, 10
sampling_04(s4_idx,m) = cloud_layer_top(n) 
sampling_04(s4_idx,m+1) = cloud_layer_base(n)
sampling_04(s4_idx,m+2) = cloud_layer_optical_depth(n)
sampling_04(s4_idx,m+3) = cloud_layer_phase(n)
sampling_04(s4_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do

end if
end if

end if ! caliop cloud flag
end if ! modis cloud flag
end if ! day flag

!----------------------------------------------------------------------------------------!
! Random Sample 5

if( day_flag == 1 )then

s5_tot = s5_tot + 1

if( s5_tot <= s5_size )then

s5_idx = s5_tot

sampling_05(s5_idx,1) = ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1)
sampling_05(s5_idx,2) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_05(s5_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_05(s5_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_05(s5_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_05(s5_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_05(s5_idx,7) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_05(s5_idx,8) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_05(s5_idx,9) = modis_c6_myd03_latitude(k,l)
sampling_05(s5_idx,10) = modis_c6_myd03_longitude(k,l)
sampling_05(s5_idx,11) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_05(s5_idx,12) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_05(s5_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_05(s5_idx,14) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_05(s5_idx,15) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_05(s5_idx,16) = surface_flag	
sampling_05(s5_idx,17) = nb_cloud_layers
sampling_05(s5_idx,18) = caliop_1km_nb_cloud_layers
sampling_05(s5_idx,19) = caliop_5km_nb_cloud_layers
sampling_05(s5_idx,20) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 2
sampling_05(s5_idx,m+n) = EV_250_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_250_Aggr1km_RefSB(k,l,n) - EV_250_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 22
do n = 1, 5
sampling_05(s5_idx,m+n) = EV_500_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_500_Aggr1km_RefSB(k,l,n) - EV_500_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 27
do n = 1, 15
sampling_05(s5_idx,m+n) = EV_1KM_RefSB_reflectance_scales(n) * &
& ( EV_1KM_RefSB(k,l,n) - EV_1KM_RefSB_reflectance_offsets(n) )
end do

m = 42
do n = 1, 16
sampling_05(s5_idx,m+n) = EV_1KM_Emissive_radiance_scales(n) * &
& ( EV_1KM_Emissive(k,l,n) - EV_1KM_Emissive_radiance_offsets(n) )
end do

m = 58
do n = 1, 6
sampling_05(s5_idx,m+n) = modis_c6_myd06_Atm_Corr_Refl(n,k,l)
end do

m = 65
do n = 1, 10
sampling_05(s5_idx,m) = cloud_layer_top(n) 
sampling_05(s5_idx,m+1) = cloud_layer_base(n)
sampling_05(s5_idx,m+2) = cloud_layer_optical_depth(n)
sampling_05(s5_idx,m+3) = cloud_layer_phase(n)
sampling_05(s5_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do
	
else

s5_idx = int( rand(0) * ( s5_tot - 1 ) ) + 1

if( s5_idx <= s5_size )then
	
sampling_05(s5_idx,1) = ibits(modis_c6_myd06_Cloud_Mask_1km(1,k,l),0,1)
sampling_05(s5_idx,2) = modis_c6_myd06_Cloud_Multi_Layer_Flag(k,l)
sampling_05(s5_idx,3) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),0,1)
sampling_05(s5_idx,4) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),1,1)
sampling_05(s5_idx,5) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),2,1)
sampling_05(s5_idx,6) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),3,1)
sampling_05(s5_idx,7) = ibits(modis_c6_myd06_Quality_Assurance_1km(6,k,l),4,1)
sampling_05(s5_idx,8) = modis_c6_myd06_Cloud_Phase_Optical_Properties(k,l)
sampling_05(s5_idx,9) = modis_c6_myd03_latitude(k,l)
sampling_05(s5_idx,10) = modis_c6_myd03_longitude(k,l)
sampling_05(s5_idx,11) = modis_c6_myd06_Cloud_Top_Height_1km(k,l)
sampling_05(s5_idx,12) = modis_c6_myd06_Cloud_Effective_Radius(k,l)
sampling_05(s5_idx,13) = modis_c6_myd06_Cloud_Effective_Radius_16(k,l)
sampling_05(s5_idx,14) = modis_c6_myd06_Cloud_Effective_Radius_37(k,l)
sampling_05(s5_idx,15) = modis_c6_myd06_Cloud_Optical_Thickness(k,l)
sampling_05(s5_idx,16) = surface_flag	
sampling_05(s5_idx,17) = nb_cloud_layers
sampling_05(s5_idx,18) = caliop_1km_nb_cloud_layers
sampling_05(s5_idx,19) = caliop_5km_nb_cloud_layers
sampling_05(s5_idx,20) = cldclass_lidar_nb_cloud_layers

m = 20
do n = 1, 2
sampling_05(s5_idx,m+n) = EV_250_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_250_Aggr1km_RefSB(k,l,n) - EV_250_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 22
do n = 1, 5
sampling_05(s5_idx,m+n) = EV_500_Aggr1km_RefSB_reflectance_scales(n) * &
& ( EV_500_Aggr1km_RefSB(k,l,n) - EV_500_Aggr1km_RefSB_reflectance_offsets(n) )
end do

m = 27
do n = 1, 15
sampling_05(s5_idx,m+n) = EV_1KM_RefSB_reflectance_scales(n) * &
& ( EV_1KM_RefSB(k,l,n) - EV_1KM_RefSB_reflectance_offsets(n) )
end do

m = 42
do n = 1, 16
sampling_05(s5_idx,m+n) = EV_1KM_Emissive_radiance_scales(n) * &
& ( EV_1KM_Emissive(k,l,n) - EV_1KM_Emissive_radiance_offsets(n) )
end do

m = 58
do n = 1, 6
sampling_05(s5_idx,m+n) = modis_c6_myd06_Atm_Corr_Refl(n,k,l)
end do

m = 65
do n = 1, 10
sampling_05(s5_idx,m) = cloud_layer_top(n) 
sampling_05(s5_idx,m+1) = cloud_layer_base(n)
sampling_05(s5_idx,m+2) = cloud_layer_optical_depth(n)
sampling_05(s5_idx,m+3) = cloud_layer_phase(n)
sampling_05(s5_idx,m+4) = cloud_layer_source(n)
m = m + 5	
end do

end if
end if

end if ! day flag

!----------------------------------------------------------------------------------------!
! Gridded Data

if( day_flag == 1 )then

grid_i = int( (modis_c6_myd03_latitude(k,l) + 90) / grid_resol ) + 1
grid_j = int( (modis_c6_myd03_longitude(k,l) + 180) / grid_resol ) + 1

if( grid_i > grid_i_max) grid_i = grid_i_max
if( grid_j > grid_j_max) grid_j = grid_j_max

grid_bin = grid_j + (grid_i-1) * grid_j_max

if( grid_bin > 0 .AND. grid_bin <= grid_i_max * grid_j_max )then

gridded_data(grid_bin,1) = gridded_data(grid_bin,1)  + 1

if(cloud_mask_flag == 1)then
gridded_data(grid_bin,2) = gridded_data(grid_bin,2)  + 1
end if

if(modis_multilayer_cloud_c6_flag == 2)then
gridded_data(grid_bin,3) = gridded_data(grid_bin,3)  + 1
end if

if(modis_multilayer_cloud_c6_ph_turned_off_flag == 2)then
gridded_data(grid_bin,4) = gridded_data(grid_bin,4)  + 1
end if

if(nb_cloud_layers > 1)then
gridded_data(grid_bin,5) = gridded_data(grid_bin,5)  + 1
end if

end if ! grid bin 

end if ! day flag

!----------------------------------------------------------------------------------------!

end if

end do ! Colocated Data

deallocate(EV_250_Aggr1km_RefSB, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (EV_250_Aggr1km_RefSB) ***"
deallocate(EV_500_Aggr1km_RefSB, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (EV_500_Aggr1km_RefSB) ***"
deallocate(EV_1KM_RefSB, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (EV_1KM_RefSB) ***"
deallocate(EV_1KM_Emissive, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (EV_1KM_Emissive) ***"

deallocate(modis_c6_myd03_latitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_latitude) ***"
deallocate(modis_c6_myd03_longitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_longitude) ***"
deallocate(modis_c6_myd03_sensorzenith, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_sensorzenith) ***"
deallocate(modis_c6_myd03_sensorazimuth, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_sensorazimuth) ***"
deallocate(modis_c6_myd03_solarzenith, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_solarzenith) ***"
deallocate(modis_c6_myd03_solarazimuth, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd03_solarazimuth) ***"

deallocate(modis_c6_myd06_Cloud_Phase_Infrared_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Phase_Infrared_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Phase_Optical_Properties, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Phase_Optical_Properties) baseline***"
deallocate(modis_c6_myd06_Cloud_Mask_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Mask_1km) baseline***"
deallocate(modis_c6_myd06_Quality_Assurance_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Quality_Assurance_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Optical_Thickness, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Optical_Thickness) baseline***"
deallocate(modis_c6_myd06_Cloud_Optical_Thickness_PCL, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Optical_Thickness_PCL) baseline***"
deallocate(modis_c6_myd06_cloud_top_temperature_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_cloud_top_temperature_1km) baseline***"
deallocate(modis_c6_myd06_Cloud_Top_Temperature, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Top_Temperature) baseline***"
deallocate(modis_c6_myd06_Cloud_Multi_Layer_Flag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Multi_Layer_Flag) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius_16, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius_16) baseline***"
deallocate(modis_c6_myd06_Cloud_Effective_Radius_37, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Effective_Radius_37) baseline***"
deallocate(modis_c6_myd06_Cloud_Top_Height_1km, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Cloud_Top_Height_1km) baseline***"
deallocate(modis_c6_myd06_Atm_Corr_Refl, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (modis_c6_myd06_Atm_Corr_Refl) baseline***"

deallocate(caliop_v4_01kmCLay_latitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_latitude)***"
deallocate(caliop_v4_01kmCLay_longitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_longitude)***"
deallocate(caliop_v4_01kmCLay_Opacity_Flag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Opacity_Flag)***"
deallocate(caliop_v4_01kmCLay_Number_Layers_Found, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Number_Layers_Found)***"
deallocate(caliop_v4_01kmCLay_Feature_Classification_Flags, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Feature_Classification_Flags)***"
deallocate(caliop_v4_01kmCLay_Layer_Top_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Layer_Top_Altitude)***"
deallocate(caliop_v4_01kmCLay_Layer_Base_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_01kmCLay_Layer_Base_Altitude)***"

deallocate(caliop_v4_05kmCLay_Number_Layers_Found, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Number_Layers_Found)***"
deallocate(caliop_v4_05kmCLay_Feature_Classification_Flags, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Feature_Classification_Flags)***"
deallocate(caliop_v4_05kmCLay_Feature_Optical_Depth_532, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Feature_Optical_Depth_532)***"
deallocate(caliop_v4_05kmCLay_Layer_Top_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Layer_Top_Altitude)***"
deallocate(caliop_v4_05kmCLay_Layer_Base_Altitude, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (caliop_v4_05kmCLay_Layer_Base_Altitude)***"

deallocate(cldclass_lidar_Height, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_Height)***"
deallocate(cldclass_lidar_CloudLayerBase, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudLayerBase)***"
deallocate(cldclass_lidar_LayerBaseFlag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_LayerBaseFlag)***"
deallocate(cldclass_lidar_CloudLayerTop, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudLayerTop)***"
deallocate(cldclass_lidar_LayerTopFlag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_LayerTopFlag)***"
deallocate(cldclass_lidar_CloudFraction, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudFraction)***"
deallocate(cldclass_lidar_CloudPhase, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudPhase)***"
deallocate(cldclass_lidar_CloudPhaseConfidenceLevel, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudPhaseConfidenceLevel)***"
deallocate(cldclass_lidar_CloudLayerType, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudLayerType)***"
deallocate(cldclass_lidar_CloudTypeQuality, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_CloudTypeQuality)***"
deallocate(cldclass_lidar_PrecipitationFlag, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_PrecipitationFlag)***"
deallocate(cldclass_lidar_Phase_log, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_Phase_log)***"
deallocate(cldclass_lidar_Water_layer_top, STAT = DeAllocateStatus)
if (DeAllocateStatus /= 0) STOP "*** Not enough memory (cldclass_lidar_Water_layer_top)***"

end if ! if input files exist

deallocate(modis_index, STAT = DeAllocateStatus)
deallocate(caliop_index, STAT = DeAllocateStatus)
deallocate(cloudsat_index, STAT = DeAllocateStatus)

end if !----- if nb_colocated_data > 0
end if !----- if found colocated file

!write(6,*) 'hour, minute: ', hour, minute

!minute = minute + 5
!if ( minute > 55 ) then
!	hour = hour + 1
!	minute = 0
!end if

!end do ! MODIS Granule Iteration
!write(6,*) 'day: ', day

!end do ! Day Iteration

!----------------------------------------------------------------------------------------!
! Random Sampling

array_shape = shape(sampling_00)
open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_00.txt')
write(1,*) s0_tot, s0_size
do i = 1, s0_size
write(1,*) (sampling_00(i,j),j=1,array_shape(2))
end do
close(1)

array_shape = shape(sampling_01)
open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_01.txt')
write(1,*) s1_tot, s1_size
do i = 1, s1_size
write(1,*) (sampling_01(i,j),j=1,array_shape(2))
end do
close(1)

array_shape = shape(sampling_02)
open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_02.txt')
write(1,*) s2_tot, s2_size
do i = 1, s2_size
write(1,*) (sampling_02(i,j),j=1,array_shape(2))
end do
close(1)

array_shape = shape(sampling_03)
open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_03.txt')
write(1,*) s3_tot, s3_size
do i = 1, s3_size
write(1,*) (sampling_03(i,j),j=1,array_shape(2))
end do
close(1)

array_shape = shape(sampling_04)
open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_04.txt')
write(1,*) s4_tot, s4_size
do i = 1, s4_size
write(1,*) (sampling_04(i,j),j=1,array_shape(2))
end do
close(1)

array_shape = shape(sampling_05)
open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_file_prefix//'sample_05.txt')
write(1,*) s5_tot, s5_size
do i = 1, s5_size
write(1,*) (sampling_05(i,j),j=1,array_shape(2))
end do
close(1)
!----------------------------------------------------------------------------------------!

open(1,file='./outputs/modis_caliop_cloudsat_multilayer_clouds/'//output_filename)
write(1,*) grid_i_max, grid_j_max
do i = 1, grid_i_max * grid_j_max
write(1,*) (gridded_data(i,j),j=1,gridded_data_shape(2))
end do
close(1)

!----------------------------------------------------------------------------------------!

end program create_modis_caliop_cloudsat_daily_samples

Create a fortran 90 module to read MODIS MYD06 C6 data

					
MODULE get_modis_c6_myd06_data

implicit none

!----------------------------------------------------------------------------------------!
! Dynamic Allocation of Arrays

integer*1, dimension(:,:,:), allocatable :: modis_c6_myd06_Cloud_Mask_1km 
integer*1, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Phase_Optical_Properties 
integer*1, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Phase_Infrared_1km 
integer*1, dimension(:,:,:), allocatable :: modis_c6_myd06_Quality_Assurance_1km 
integer*1, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Multi_Layer_Flag 
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Top_Temperature 
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_cloud_top_temperature_1km 
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Effective_Radius 
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Effective_Radius_16 
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Effective_Radius_37
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Optical_Thickness
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Optical_Thickness_PCL
integer*2, dimension(:,:), allocatable :: modis_c6_myd06_Cloud_Top_Height_1km 
integer*2, dimension(:,:,:), allocatable :: modis_c6_myd06_Atm_Corr_Refl

!----------------------------------------------------------------------------------------!
! User Inputs (user_selected_sds)

integer :: nb_user_selected_sds
parameter ( nb_user_selected_sds = 14 )

character*300 selected_sds( nb_user_selected_sds )      

data selected_sds  / 'Cloud_Phase_Infrared_1km', 'Cloud_Phase_Optical_Properties', &
                   & 'Cloud_Mask_1km', 'Quality_Assurance_1km', &
                   & 'Cloud_Optical_Thickness', 'Cloud_Optical_Thickness_PCL', &
                   & 'cloud_top_temperature_1km', 'Cloud_Top_Temperature', &
                   & 'Cloud_Multi_Layer_Flag', 'Cloud_Effective_Radius', &
                   & 'Cloud_Effective_Radius_16', 'Cloud_Effective_Radius_37', &
                   & 'cloud_top_height_1km', 'Atm_Corr_Refl' /

!----------------------------------------------------------------------------------------!
! Subroutine to read selected SDSs

CONTAINS

subroutine read_modis_c6_myd06(modis_c6_myd06_file_name)

use read_write_hdf_file_variables

implicit none

integer :: i,j,k,l,m,n
integer :: AllocateStatus, DeAllocateStatus
integer :: ierr ! etat de reussite (si ==0) ou echec (si /= 0)
character*300 :: modis_c6_myd06_file_name

integer sfstart, sfselect, sfrdata, sfendacc, sfend , sffattr , sfrattr, var_id, attr_id, err_code
integer status, n_attrs
integer sffinfo , sfginfo

sd_id = sfstart(modis_c6_myd06_file_name, DFACC_READ)
status = sffinfo( sd_id , n_datasets , n_file_attrs )

do index = 0, n_datasets - 1

	sds_id = sfselect(sd_id, index)
	status = sfginfo( sds_id , name , rank , dim_sizes , data_type , n_attrs )
    sds_name = name(1:len(name))

if ( ANY( selected_sds==sds_name ) ) then

	do i = 1 , rank
	start(i) = 0
	edges(i) = dim_sizes(i)
	stride(i) = 1
	end do

	if( sds_name == "Cloud_Phase_Infrared_1km" )then
	allocate( modis_c6_myd06_Cloud_Phase_Infrared_1km(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Phase_Infrared_1km***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Phase_Infrared_1km )
	end if

	if( sds_name == 'Cloud_Phase_Optical_Properties' )then
	allocate(   modis_c6_myd06_Cloud_Phase_Optical_Properties(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus modis_c6_myd06_Cloud_Phase_Optical_Properties)***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Phase_Optical_Properties )
	end if

	if( sds_name == 'Cloud_Mask_1km' )then
	allocate(   modis_c6_myd06_Cloud_Mask_1km(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Mask_1km***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Mask_1km )
	end if

	if( sds_name == 'Quality_Assurance_1km' )then
	allocate(   modis_c6_myd06_Quality_Assurance_1km(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Quality_Assurance_1km***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Quality_Assurance_1km )
	end if

	if( sds_name == 'Cloud_Optical_Thickness' )then
	allocate( modis_c6_myd06_Cloud_Optical_Thickness(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Optical_Thickness***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Optical_Thickness )
	end if
	
	if( sds_name == 'Cloud_Optical_Thickness_PCL' )then
	allocate( modis_c6_myd06_Cloud_Optical_Thickness_PCL(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Optical_Thickness_PCL***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Optical_Thickness_PCL )
	end if	

	if( sds_name == 'cloud_top_temperature_1km' )then
	allocate( modis_c6_myd06_cloud_top_temperature_1km(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_cloud_top_temperature_1km***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_cloud_top_temperature_1km )
	end if

	if( sds_name == 'Cloud_Top_Temperature' )then
	allocate( modis_c6_myd06_Cloud_Top_Temperature(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Top_Temperature***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Top_Temperature )
	end if

	if( sds_name == 'Cloud_Multi_Layer_Flag' )then
	allocate( modis_c6_myd06_Cloud_Multi_Layer_Flag(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Multi_Layer_Flag***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Multi_Layer_Flag )
	end if	

	if( sds_name == 'Cloud_Effective_Radius' )then
	allocate( modis_c6_myd06_Cloud_Effective_Radius(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Effective_Radius***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Effective_Radius )
	end if	

	if( sds_name == 'Cloud_Effective_Radius_16' )then
	allocate( modis_c6_myd06_Cloud_Effective_Radius_16(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Effective_Radius_16***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Effective_Radius_16 )
	end if	

	if( sds_name == 'Cloud_Effective_Radius_37' )then
	allocate( modis_c6_myd06_Cloud_Effective_Radius_37(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Effective_Radius_37***"
	status = sfrdata( sds_id, start, stride, edges, modis_c6_myd06_Cloud_Effective_Radius_37 )
	end if	

	if( sds_name == 'cloud_top_height_1km' )then
	allocate(  modis_c6_myd06_Cloud_Top_Height_1km(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Cloud_Top_Height_1km***"
	status = sfrdata( sds_id, start, stride, edges,  modis_c6_myd06_Cloud_Top_Height_1km )
	end if	

	if( sds_name == 'Atm_Corr_Refl' )then
	allocate(  modis_c6_myd06_Atm_Corr_Refl(dim_sizes(1),dim_sizes(2),dim_sizes(3)), STAT = AllocateStatus )
	IF (AllocateStatus /= 0) STOP "*** Not enough memory modis_c6_myd06_Atm_Corr_Refl***"
	status = sfrdata( sds_id, start, stride, edges,  modis_c6_myd06_Atm_Corr_Refl )
	end if	
	
end if

end do

!----------------------------------------------------------------------------------------!

status = sfendacc(sds_id)
status = sfend(sd_id)

!----------------------------------------------------------------------------------------!

end subroutine read_modis_c6_myd06

!----------------------------------------------------------------------------------------!

END MODULE get_modis_c6_myd06_data