Notebook Source Codes

Share your source codes

Source codes: 118

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')

a = Latitude[:]

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

#----------------------------------------------------------------------------------------#
# 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)

plt.title('cldclass-lidar')

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)

plt.title('cldclass-lidar')

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)

plt.title('cldclass-lidar')

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)

plt.title('cldclass-lidar')

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

Python script modis cldclass-lidar analysis (main)

					
#!/usr/bin/env python

import numpy as np

import modis_cldclass_lidar_analysis_001 as analysis_001

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

dict_settings = {}

dict_settings['input_path'] = ''

dict_settings['output_path'] = './media/images/modis_cldclass_lidar_analysis_main_py/'
dict_settings['output_file_name_prefix'] = 'modis_cldclass_lidar_analyisis_'

dict_settings['input_file_name'] = '2008_01_sample_01_modis_caliop_cldclass_lidar.txt'

data = np.loadtxt(
dict_settings['input_path'] + dict_settings['input_file_name'], skiprows=1)

fig_count = 0

#----------------------------------------------------------------------------------------#
# confusion matrices 

dict_settings['fig_title'] = 'modis c6 vs cldclass-lidar \n score = '
dict_settings['fig_ylabel'] = 'cldclass-lidar'
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['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 \n score = '
dict_settings['fig_ylabel'] = 'caliop 1km'
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['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 \n score = '
dict_settings['fig_ylabel'] = 'caliop 5km'
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['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 \n score = '
dict_settings['fig_ylabel'] = 'caliop 1km'
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['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 \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)

#----------------------------------------------------------------------------------------#
# 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.cloud_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.cloud_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.cloud_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.cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

Python script modis cldclass-lidar analysis 001

					
#!/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

#----------------------------------------------------------------------------------------#
# 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((2,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
		
			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 / sum
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = ['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], ['modis C6 \n multilayer cloud', 'modis C6 \n monolayer cloud'],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 (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

#----------------------------------------------------------------------------------------#
# Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def 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']]
	
	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)]
		
	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('------')
		
			cloud_layer_caliop_count_bar = [0 for i in range(200)]

			'''
			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
			'''

			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
				
	
	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()

	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

Create modis caliop cloudsat matched files for a month with python 3

					
#!/usr/bin/env python

import os
import glob
import numpy as np
import calendar
import re

#----------------------------------------------------------------------------------------#
# User inputs

year = 2008
month = 7

machine = 'nasa' # icare

create_database_directory = False

#----------------------------------------------------------------------------------------#
# local machine (paths)

if machine == 'nasa':

	path_to_myd021km = '/Volumes/Data/MYD021KM/'
	path_to_myd03 = '/Volumes/Data/Myd03/'
	path_to_myd35 = '/Volumes/DATA_03/modis/Myd35/'
	path_to_myd06_c5 = '/Volumes/DATA_03/modis/Myd06_c5/'
	path_to_myd06_c6 = '/Volumes/DATA_03/modis/Myd06_c6/'
	path_to_01km_clay = '/Volumes/Data/Caliop/v4/01km_CLay/'
	path_to_05km_clay = '/Volumes/Data/Caliop/v4/05km_CLay/'
	path_to_cloudsat = '/Volumes/Data/Cloudsat/'
	path_to_2b_cldclass_lidar = '/Volumes/Data/2B-CLDCLASS-LIDAR/'

	dest_path = '/Users/marchant/Desktop/modis_based_match_files/'

#----------------------------------------------------------------------------------------#
# icare server (paths)

if machine == 'icare':

	path_to_myd021km = '/DATA/LIENS/MODIS/MYD021KM.006/'
	path_to_myd03 = '/DATA/LIENS/MODIS/MYD03/'
	path_to_myd35 = '/DATA/LIENS/MODIS/MYD35_L2.006/'
	path_to_myd06_c5 = '/DATA/LIENS/MODIS/MYD06_L2.005/'
	path_to_myd06_c6 = '/DATA/LIENS/MODIS/MYD06_L2.006/'
	path_to_01km_clay = '/DATA/LIENS/CALIOP/01kmCLay.v4.10/'
	path_to_05km_clay = '/DATA/LIENS/CALIOP/05kmCLay.v4.10/'
	path_to_cloudsat = '/DATA/LIENS/CLOUDSAT/2B-GEOPROF/'
	path_to_2b_cldclass_lidar = '/DATA/LIENS/CLOUDSAT/2B-CLDCLASS-LIDAR.v04.02/'

	dest_path = '/home/benjamin.hg.marchant/data/modis_based_match_files/'

#----------------------------------------------------------------------------------------#
# Create directory (Optional)

if create_database_directory:

	os.system( 'mkdir ' + dest_path + str(year) )

	for day in np.arange(calendar.monthrange(year,month)[1]):
		d = './' + str(year) + '/' + str(year) + '_' + "%02d" % month + \
			 '_' + "%02d" % (day+1) + '/'
		os.system('mkdir ' + dest_path + d)

#----------------------------------------------------------------------------------------#
# create file

day_of_year_0 = 0
for i in range(month-1):
	day_of_year_0 = day_of_year_0 + calendar.monthrange(year,i+1)[1]

for day in np.arange(calendar.monthrange(year,month)[1]):
#for day in [0]: #np.arange(calendar.monthrange(year,month)[1]):

	day_of_year = day_of_year_0 + day + 1

	#------------------------------------------------------------------------------------# 
	# create paths

	d = str(year) + '/' + str(year) + '_' + "%02d" % month + '_' + "%02d" % (day+1) + '/'

	path_to_myd021km_day = path_to_myd021km + d
	path_to_myd03_day = path_to_myd03 + d
	path_to_myd35_day = path_to_myd35 + d
	path_to_myd06_c5_day = path_to_myd06_c5 + d
	path_to_myd06_c6_day = path_to_myd06_c6 + d			
	path_to_01km_clay_day = path_to_01km_clay + d
	path_to_05km_clay_day = path_to_05km_clay + d
	path_to_cloudsat_day = path_to_cloudsat + d
	path_to_2b_cldclass_lidar_day = path_to_2b_cldclass_lidar + d

	#------------------------------------------------------------------------------------# 
	# get all file names

	list_of_myd021km_files = glob.glob(path_to_myd021km_day + '/*hdf' )
	list_of_myd03_files = glob.glob(path_to_myd03_day + '/*hdf' )
	list_of_myd35_files = glob.glob(path_to_myd35_day + '/*hdf' )
	list_of_myd06_c5_files = glob.glob(path_to_myd06_c5_day + '/*hdf' )
	list_of_myd06_c6_files = glob.glob(path_to_myd06_c6_day + '/*hdf' )
	list_of_01km_clay_files = glob.glob(path_to_01km_clay_day + '/*hdf' )
	list_of_05km_clay_files = glob.glob(path_to_05km_clay_day + '/*hdf' )	
	list_of_cloudsat_files = glob.glob(path_to_cloudsat_day + '/*hdf' ) 
	list_of_2b_cldclass_lidar_files = glob.glob(path_to_2b_cldclass_lidar_day + '/*hdf' )  

	#------------------------------------------------------------------------------------# 
	# create file

	myd021km_name = 'none'		
	myd03_name = 'none'	
	myd35_name = 'none'					
	myd06_c5_name = 'none'				
	myd06_c6_name = 'none'				
	cal_1km_file = 'none'	
	cal_5km_file = 'none'	
	cloudsat_2B_geoprof_file = 'none'	
	cloudsat_2B_cldclass_lidar_file = 'none'	

	for modis_granule_name in list_of_myd03_files:
		
		myd03_name = modis_granule_name.split('/')[-1]
		
		hm = myd03_name.split('.')[2]

		modis_granule_hour = int(hm[0:2])
		modis_granule_minute = int(hm[2:4])

		#--------------------------------------------------------------------------------#
		# Caliop 01km_clay

		min = 200

		file_exists = False

		for caliop_01km_granule_name in list_of_01km_clay_files:
			
			caliop_granule_name = caliop_01km_granule_name.split('/')[-1]
			
			caliop_re = re.compile(str(year) + '-' + "%02d" % month + '-' + \
			                       "%02d" % (day+1) + "T(.*)-(.*)-(.*).hdf")

			obj = caliop_re.search(caliop_granule_name)

			caliop_granule_hour = int(obj.group(1))
			caliop_granule_minute = int(obj.group(2))
			
			diff = int(modis_granule_hour) * 60 + int(modis_granule_minute) \
			        - ( int(caliop_granule_hour) * 60 + int(caliop_granule_minute) )

			if diff < min and diff > 0 :
				min = diff
				cal_1km_file = caliop_granule_name
				file_exists = True

		#--------------------------------------------------------------------------------#
		# Caliop 05km_clay
			
		min = 200

		file_exists = False

		for caliop_05km_granule_name in list_of_05km_clay_files:
			
			caliop_granule_name = caliop_05km_granule_name.split('/')[-1]
			
			caliop_re = re.compile(str(year) + '-' + "%02d" % month + '-' + \
			                       "%02d" % (day+1) + "T(.*)-(.*)-(.*).hdf")

			obj = caliop_re.search(caliop_granule_name)

			caliop_granule_hour = int(obj.group(1))
			caliop_granule_minute = int(obj.group(2))

			diff = int(modis_granule_hour) * 60 + int(modis_granule_minute) \
			        - ( int(caliop_granule_hour) * 60 + int(caliop_granule_minute) )

			if diff < min and diff > 0 :
				min = diff
				cal_5km_file = caliop_granule_name
				file_exists = True

		#--------------------------------------------------------------------------------#
		# MYD021KM
		
		file_name_prefix = 'MYD021KM.A'+ str(year) + \
		                   "%03d" % day_of_year + \
		                   '.' + "%02d" % modis_granule_hour + \
		                   "%02d" % modis_granule_minute
		found_file = 0
		for myd021km_granule_name in list_of_myd021km_files:
			offset = len(path_to_myd021km_day)
			if file_name_prefix == myd021km_granule_name[offset:offset+22]:
				found_file = 1
				myd021km_name = myd021km_granule_name[offset:]
		if found_file == 0: myd021km_name = 'NoFile'

		#--------------------------------------------------------------------------------#
		# MYD35

		file_name_prefix = 'MYD35_L2.A'+ str(year) + \
		                   "%03d" % day_of_year + \
		                   '.' + "%02d" % modis_granule_hour + \
		                   "%02d" % modis_granule_minute
		found_file = 0
		for myd35_granule_name in list_of_myd35_files:
			offset = len(path_to_myd35_day)
			if file_name_prefix == myd35_granule_name[offset:offset+22]:
				found_file = 1
				myd35_name = myd35_granule_name[offset:]
		if found_file == 0: myd35_name = 'NoFile'

		#--------------------------------------------------------------------------------#
		# MYD06_c5

		file_name_prefix = 'MYD06_L2.A'+ str(year) + \
		                   "%03d" % day_of_year + \
		                   '.' + "%02d" % modis_granule_hour + \
		                   "%02d" % modis_granule_minute
		found_file = 0		
		for myd06_c5_granule_name in list_of_myd06_c5_files:
			offset = len(path_to_myd06_c5_day)
			if file_name_prefix == myd06_c5_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_c5_name = myd06_c5_granule_name[offset:]
		if found_file == 0: myd06_c5_name = 'NoFile'		

		#--------------------------------------------------------------------------------#
		# MYD06_c6

		file_name_prefix = 'MYD06_L2.A'+ str(year) + \
		                   "%03d" % day_of_year + \
		                   '.' + "%02d" % modis_granule_hour + \
		                   "%02d" % modis_granule_minute
		found_file = 0		
		for myd06_c6_granule_name in list_of_myd06_c6_files:
			offset = len(path_to_myd06_c6_day)
			if file_name_prefix == myd06_c6_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_c6_name = myd06_c6_granule_name[offset:]
		if found_file == 0: myd06_c6_name = 'NoFile'	

		#--------------------------------------------------------------------------------#
		# cloudsat (GEOPROFF)

		find_cloudsat_match_file = False
		min_minute_diff = 99
		for cloudsat_file in list_of_cloudsat_files:
			
			cloudsat_file_name = cloudsat_file.split('/')[-1]
			cloudsat_granule_hour = int(cloudsat_file_name[7:9])
			cloudsat_granule_minute = int(cloudsat_file_name[9:11])
			cloudsat_granule_second = int(cloudsat_file_name[11:13] )
			
			minute_diff = modis_granule_hour * 60 + modis_granule_minute - \
					( cloudsat_granule_hour * 60 + cloudsat_granule_minute )
			
			if (minute_diff < min_minute_diff ) and ( minute_diff > 0 ):
					min_minute_diff = minute_diff
					cloudsat_2B_geoprof_file = cloudsat_file_name
					find_cloudsat_match_file = True

		#--------------------------------------------------------------------------------#
		# cloudsat (2B CLDCLASS-LIDAR)

		find_2b_cldclass_lidar_match_file = False
		min_minute_diff = 99
		for cloudsat_file in list_of_2b_cldclass_lidar_files:

			cloudsat_file_name = cloudsat_file.split('/')[-1]
			cloudsat_granule_hour = int(cloudsat_file_name[7:9])
			cloudsat_granule_minute = int(cloudsat_file_name[9:11])
			cloudsat_granule_second = int(cloudsat_file_name[11:13] )

			minute_diff = modis_granule_hour * 60 + modis_granule_minute - \
					( cloudsat_granule_hour * 60 + cloudsat_granule_minute )

			if (minute_diff < min_minute_diff ) and ( minute_diff > 0 ):
					min_minute_diff = minute_diff
					cloudsat_2B_cldclass_lidar_file = cloudsat_file_name
					find_2b_cldclass_lidar_match_file = True

		#--------------------------------------------------------------------------------#
		# save result

		filename = 'modis_linked_files_' + \
		           "%02d" % modis_granule_hour + "%02d" % modis_granule_minute
		
		f = open(filename, 'w')

		f.write(myd021km_name + ' \n')			
		f.write(myd03_name + ' \n')	
		f.write(myd35_name + ' \n')				
		f.write(myd06_c5_name + ' \n')			
		f.write(myd06_c6_name + ' \n')			
		f.write(cal_1km_file + ' \n')
		f.write(cal_5km_file + ' \n')
		f.write(cloudsat_2B_geoprof_file + ' \n')
		f.write(cloudsat_2B_cldclass_lidar_file + ' \n')
		
		f.closed
		
		command = 'mv ' + filename + ' ' + dest_path + d
		
		os.system(command)

Fortran 90 module to read CloudSat 2B CLDCLASS products

					
MODULE get_cldclass_lidar_data

implicit none

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

integer*2, dimension(:,:), allocatable :: cldclass_lidar_Height
real(kind=4), dimension(:,:), allocatable :: cldclass_lidar_CloudLayerBase
integer*1, dimension(:,:), allocatable :: cldclass_lidar_LayerBaseFlag
real(kind=4), dimension(:,:), allocatable :: cldclass_lidar_CloudLayerTop
integer*1, dimension(:,:), allocatable :: cldclass_lidar_LayerTopFlag
real(kind=4), dimension(:,:), allocatable :: cldclass_lidar_CloudFraction
integer*1, dimension(:,:), allocatable :: cldclass_lidar_CloudPhase
integer*1, dimension(:,:), allocatable :: cldclass_lidar_CloudPhaseConfidenceLevel
integer*1, dimension(:,:), allocatable :: cldclass_lidar_CloudLayerType
real(kind=4), dimension(:,:), allocatable :: cldclass_lidar_CloudTypeQuality
integer*1, dimension(:,:), allocatable :: cldclass_lidar_PrecipitationFlag
integer*1, dimension(:,:), allocatable :: cldclass_lidar_Phase_log
real(kind=4), dimension(:,:), allocatable :: cldclass_lidar_Water_layer_top

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

integer :: nb_user_selected_sds
parameter ( nb_user_selected_sds = 13 )

character*300 selected_sds( nb_user_selected_sds )      

data selected_sds  / 'Height', &
& 'CloudLayerBase', &
& 'LayerBaseFlag', &
& 'CloudLayerTop', &
& 'LayerTopFlag', &
& 'CloudFraction', &
& 'CloudPhase', &
& 'CloudPhaseConfidenceLevel', &
& 'CloudLayerType', &
& 'CloudTypeQuality', &
& 'PrecipitationFlag', &
& 'Phase_log', &
& 'Water_layer_top' /

!----------------------------------------------------------------------------------------!
! Subroutine read_cldclass_lidar

CONTAINS

subroutine read_cldclass_lidar(cldclass_lidar_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 :: cldclass_lidar_file_name

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

character*300 vdata_name 
integer file_id, vdata_id, vdata_ref
integer n_records
character*80 fieldname_lists 
integer FULL_INTERLACE
parameter (FULL_INTERLACE = 0)
integer interlace_mode, vdata_size

!----------------------------------------------------------------------------------------!
! SDS

sd_id = sfstart(cldclass_lidar_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 == 'Height' )then
	allocate( cldclass_lidar_Height(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_Height)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_Height )
	end if	
	
	if( sds_name == 'CloudLayerBase' )then
	allocate( cldclass_lidar_CloudLayerBase(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudLayerBase)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudLayerBase )
	end if	
	
	if( sds_name == 'LayerBaseFlag' )then
	allocate( cldclass_lidar_LayerBaseFlag(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_LayerBaseFlag)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_LayerBaseFlag )
	end if	
	
	if( sds_name == 'CloudLayerTop' )then
	allocate( cldclass_lidar_CloudLayerTop(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudLayerTop)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudLayerTop )
	end if	
	
	if( sds_name == 'LayerTopFlag' )then
	allocate( cldclass_lidar_LayerTopFlag(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_LayerTopFlag)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_LayerTopFlag )
	end if	
	
	if( sds_name == 'CloudFraction' )then
	allocate( cldclass_lidar_CloudFraction(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudFraction)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudFraction )
	end if	
	
	if( sds_name == 'CloudPhase' )then
	allocate( cldclass_lidar_CloudPhase(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudPhase)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudPhase )
	end if	
	
	if( sds_name == 'CloudPhaseConfidenceLevel' )then
	allocate( cldclass_lidar_CloudPhaseConfidenceLevel(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudPhaseConfidenceLevel)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudPhaseConfidenceLevel )
	end if	
	
	if( sds_name == 'CloudLayerType' )then
	allocate( cldclass_lidar_CloudLayerType(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudLayerType)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudLayerType )
	end if	
	
	if( sds_name == 'CloudTypeQuality' )then
	allocate( cldclass_lidar_CloudTypeQuality(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_CloudTypeQuality)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_CloudTypeQuality )
	end if	
	
	if( sds_name == 'PrecipitationFlag' )then
	allocate( cldclass_lidar_PrecipitationFlag(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_PrecipitationFlag)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_PrecipitationFlag )
	end if	
	
	if( sds_name == 'Phase_log' )then
	allocate( cldclass_lidar_Phase_log(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_Phase_log)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_Phase_log )
	end if	
	
	if( sds_name == 'Water_layer_top' )then
	allocate( cldclass_lidar_Water_layer_top(dim_sizes(1),dim_sizes(2)), STAT = AllocateStatus)
	IF (AllocateStatus /= 0) STOP "*** Not enough memory (AllocateStatus cldclass_lidar_Water_layer_top)***"
	status = sfrdata( sds_id, start, stride, edges, cldclass_lidar_Water_layer_top )
	end if	

end if

end do

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

end subroutine read_cldclass_lidar

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

END MODULE get_cldclass_lidar_data

Read 2B-CLDCLASS-LIDAR SDS names and attributes using python 3 and pyhdf

					
#!/usr/bin/env python

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

import pprint

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

file_path = './products/2B-CLDCLASS-LIDAR/'
file_name = '2008183012329_11573_CS_2B-CLDCLASS-LIDAR_GRANULE_P_R04_E02.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files (VD data) 

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

print(vs.vdatainfo() )

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

a = Latitude[:]

print('VD data')
print(len(a),Latitude[0])

print( Latitude.attrinfo() )

print(min(a),max(a))

for i in range(20000,20020):
	print(a[i])

Latitude.detach() # "close" the vdata
Longitude.detach() # "close" the vdata
vs.end() # terminate the vdata interface
f.close() 

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

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

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

#----------------------------------------------------------------------------------------#
# print sds names

datasets_dic = file.datasets()

#for idx,sds in enumerate(datasets_dic.keys()):
#	print idx,sds 

sds_dic = {}
for key, value in datasets_dic.items():
	#print( key, value, value[3] )
	sds_dic[value[3]] = key

pprint.pprint( sds_dic )

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

datasets_dic = file.datasets()

for key, value in datasets_dic.items():
	#print( key, value, value[3] )
	print('---------------------------')
	print('------------>', key)
	sds_obj = file.select(key) # select sds
	data = sds_obj.get() 
	sds_info = sds_obj.info()
	print( 'sds info' )
	print( sds_info )
	print( sds_info[0], sds_info[1] )
	print( 'sds attributes' )
	pprint.pprint( sds_obj.attributes() )

About


To add a source code to your notebook just click on the link below. Sharing scientific source code allows to understand and reproduce your results:

Add a code

Search


Source codes can be retrieved using the following form: