dropsonde_raw_gap_filler.py

from __future__ import print_function, division
import numpy as np
import glob
import copy
import os
import pdb
import warnings
from general_importer import * 		# General Importer reporting for duty.
import datetime
import matplotlib
import matplotlib.pyplot as plt



def run_dropsonde_gap_filler(path_raw_sondes, data_out_path_halo,
	dropsonde_dataset='raw', path_BAH_data=''):

	"""
	Parameters
	----------
	path_raw_sondes : str
		Path of raw dropsonde data.
	data_out_path_halo : str
		Path of repaired dropsonde data (gaps filled).
	dropsonde_dataset : str, optional
		Dropsonde dataset specification. 'raw' is default.
	path_BAH_data : str, optional
		Path of BAHAMAS data in unified netCDF files.
	"""

	def run_dropsonde_gap_filler_raw(path_raw_sondes, data_out_path_halo,
		path_BAH_data=''):

		"""
		Parameters
		----------
		path_raw_sondes : str
			Path of raw dropsonde data.
		data_out_path_halo : str
			Path of repaired dropsonde data (gaps filled).
		path_BAH_data : str, optional
			Path of BAHAMAS data in unified netCDF files.
		"""


		############################################################################################
		# FUNCTIONS

		def fill_gaps(old_var):
		# old variable gets linearly interpolated for each sonde launch. The function is ignoring nan values at the surface and above
		# the launch altitude.

			new_var = copy.deepcopy(old_var)
			# create flag variable indicating if an entry of old_var has been changed: if = 0: not interpol.
			interp_flag = np.zeros(old_var.shape)


			# identify regions of nan values in the middle of the drop. Extrapolation will be handled in another function.
			# identify the highest non-nan entry so we can cut the values above that highest entry:
			# identify the lowest non-nan entry for similar reasons:
			non_nan_idx = [idx for idx, x in np.ndenumerate(old_var) if (not np.isnan(x))]
			non_nan_idx = np.where(~np.isnan(old_var))[0]
			limits = np.array([non_nan_idx[0], non_nan_idx[-1]])

			temp_var = copy.deepcopy(old_var)
			temp_var = temp_var[limits[0]:limits[1]+1]		# will be the variable where the gaps are filled

			interp_flag_temp = np.zeros(temp_var.shape)

			# identify mid-drop-nan-values: need values after and before the nan:
			nan_idx = np.argwhere(np.isnan(temp_var))
			interp_flag_temp[nan_idx] = 1

			if nan_idx.size == 0:
				return new_var, interp_flag

			else: # correct nan values: find the hole size via subtraction of subsequent indices
				hole_size = np.zeros((len(nan_idx)+1,)).astype(int)
				# hole_size = np.zeros(nan_idx.shape)			# old version
				k = 0		# index to address a hole ('hole number')

				for m in range(0, len(temp_var)-1):

					if not np.isnan(temp_var[m+1] - temp_var[m]):
						hole_size[k] = 0

					elif np.isnan(temp_var[m+1] - temp_var[m]): # k shall only be incremented if an END of a hole has been identified:
						if len(nan_idx) == 1: 	# must be handled seperately in case that merely one nan value exists in temp_var
							hole_size[k] = 1
							break

						else:
							if (not np.isnan(temp_var[m+1])) & (np.isnan(temp_var[m])): # END of a hole
								k = k + 1
								continue
							hole_size[k] = hole_size[k] + 1		# k won't be incremented until m finds another non-nan value

					else:
						print("\n Something unexpected happened when trying to find the nan values in the middle of the dropsonde launch... Contact 'a.walbroel@uni-koeln.de'. \n")

				# holes have been identified: edit the FIRST hole (editing depends on the size of the hole...)
				c = 0 		# dummy variable needed for the right jumps in hole_size and nan_idx. c is used to address nan_idx and therefore new_var...

				# meanwhile 'd' just runs through the array hole_size:
				for d in range(0, len(hole_size)):
					for L in range(0, hole_size[d]):		# range(0, 1): L = 0
						temp_var[nan_idx[c] + L] = temp_var[nan_idx[c] - 1] + (L + 1)*(temp_var[int(nan_idx[c] + hole_size[d])] - temp_var[nan_idx[c]-1]) / (hole_size[d] + 1)

					c = c + int(hole_size[d])
					if c > len(hole_size)-1:
						break


			# overwrite the possibly holey section:
			new_var[limits[0]:limits[1]+1] = temp_var
			# update interp_flag
			interp_flag[limits[0]:limits[1]+1] = interp_flag_temp

			return new_var, interp_flag


		def std_extrapol_BAH(old_dict, ill_keys, bah_filename, old_ipflag_dict=dict()):
		# Will extrapolate some atmospheric variables to the ceiling of the dropsonde; old_ipflag will be updated.
		# Needs the old variable, the interpolation flag (should've been generated by fill_gaps()), the key and height levels as INPUT

			new_dict = old_dict
			n_alt = len(new_dict['Z'])

			new_ipflag_dict = old_ipflag_dict


			# Need BAHAMAS information to set the new ceiling:
			# Import altitude and time data from BAHAMAS:
			bah_keys = ['time', 'altitude', 'ta', 'p', 'rh']
			bah_dict = import_BAHAMAS_unified(bah_filename[0], bah_keys)
			bah_dict['time'] = np.rint(bah_dict['time']).astype(float) # must be done to avoid small fractions of seconds

			# to get the obs_height: average BAHAMAS altitude over +/- 10 seconds around launch_time:
			# find time index of the sonde launches:
			timestamp = new_dict['launch_time']
			bah_launch_idx = np.asarray([np.argwhere(bah_dict['time'] == timestamp)]).flatten()		# had some dimensions too many -> flattened
			drop_alt = np.floor(np.asarray([np.mean(bah_dict['altitude'][i-10:i+10]) for i in bah_launch_idx])/100)*100
			obs_height = np.max(np.unique(drop_alt))		# some values are repeated ... omit them and find the max value: this will be used for top of extrapolation!

			# BAHAMAS temperature, pressure, relative humidity at launch time:
			bah_T = np.nanmean(bah_dict['ta'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])
			bah_P = np.nanmean(bah_dict['p'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])
			bah_RH = np.nanmean(bah_dict['rh'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])
			bah_alt = np.nanmean(bah_dict['altitude'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])


			ceiling = obs_height	# last entry of altitude
			if ceiling > 15000:
				print("Ceiling appears to be > 15000 m. Aborted extrapolation to dropsonde ceiling because the tropopause may intervene.\n")
				return new_dict, new_ipflag_dict


			# any value above obs_height will be deleted: So if e.g. Z has got values above obs_height, delete them:
			# find the first index that overshoots obs_height:
			with warnings.catch_warnings():
				warnings.simplefilter("ignore")
				overshoot = np.argwhere(new_dict['Z'] >= obs_height)
			if len(overshoot) > 0:
				overshoot = overshoot[0][0] + 1

				for key in new_dict.keys():
					if key in ['trajectory', 'fillValues', 'ipflag']:	# skip these ones ... it s not interesting anyway
						continue

					if new_dict[key].ndim > 0:		# otherwise: error when using len()

						if len(new_dict[key]) == n_alt:
							new_dict[key] = new_dict[key][:overshoot]	# limit variable to obs_height
							if key in ill_keys or key == 'Z':
								new_ipflag_dict[key] = new_ipflag_dict[key][:overshoot]


			# at the end of 'Z' there may still be nans -> so we don't know to which altitude meteorological variables belong to in this region:
			# therefore: delete it and replace by extrapolation:
			n_alt = len(new_dict['Z'])
			last_nonnan_alt = np.argwhere(~np.isnan(new_dict['Z']))[-1][0]
			if ceiling - new_dict['Z'][last_nonnan_alt] > 1000:
				print("WARNING: highest GPS altitude measurement is at least 1000 m below the aircraft. Extrapolation may be erroneous.\n")

			for key in new_dict.keys():
				if key in ['trajectory', 'fillValues', 'ipflag']:	# skip these ones ... it s not interesting anyway
					continue

				if new_dict[key].ndim > 0:		# otherwise: error when using len()

					if len(new_dict[key]) == n_alt:
						new_dict[key] = new_dict[key][:last_nonnan_alt+1]	# limit variable to obs_height
						if key in ill_keys or key == 'Z':
							new_ipflag_dict[key] = new_ipflag_dict[key][:last_nonnan_alt+1]


			# extend the old height grid up to the ceiling if the distance is greater than 10 meters:
			alt = new_dict['Z']
			n_alt = len(alt)
			alt = np.append(alt, np.arange(alt[np.argwhere(~np.isnan(alt))[-1]]+10, ceiling+11, 10))
			n_alt_new = len(alt)

			# update the altitude variable in the dictionary: & update ipflag for gpsalt:
			new_dict['Z'] = alt
			new_ipflag_dict['Z'] = np.append(new_ipflag_dict['Z'], np.ones((n_alt_new - n_alt,)))


			np.random.seed(42)		# needed for the added noise
			launch_time = datetime.datetime.utcfromtimestamp(new_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S") # for printing

			for key in ill_keys:

				new_var = new_dict[key]
				# must be expanded to the new height grid (up to the new ceiling)
				new_var = np.append(new_var, np.full((n_alt_new - n_alt,), np.nan), axis=0) # append nans at the top of the profile

				if not new_ipflag_dict: # in case fill_gaps(...) wasn't called before this one, it's assumed that nothing has been interpolated yet.
					new_ipflag_dict[key] = np.zeros(new_var.shape)

				else: # new_ipflag also has to be extended to the new hgt grid:
					new_ipflag_dict[key] = np.append(new_ipflag_dict[key], np.zeros((n_alt_new - n_alt,)), axis=0)



				if key == 'T':
					# If BAHAMAS Temperature measurement is available use it as target in case only the top 15 % of measurements
					# are missing. Otherwise:
					# Temperature: If dropsondes with measurements (ipflag = 0) from that day exist, estimate their average T gradient.
					# If the extrapolated dropsonde temperature then deviates from BAH T by more than 5 K, use the ICAO std atmosphere
					# as T gradient:
					# Standard atmosphere (shifted accordingly to avoid a jump between the last known value and the extrapolation).
					# ICAO standard atmosphere taken from:
					# https://www.dwd.de/DE/service/lexikon/begriffe/S/Standardatmosphaere_pdf.pdf?__blob=publicationFile&v=3
					ICAO_standard_T = 288.15 - 0.0065*alt
					noise_strength = 0/2


					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.6*ceiling:
						print("Insufficient amount of measurements for temperature extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no temperature measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var			# then just overwrite the dictionary entry with the nonedited (but extended) variable
						continue

					if alt[idx] < 0.85*ceiling: # then use BAHAMAS temperature as extrapolation target:
						new_var[idx+1:] = (new_var[idx] + (bah_T - new_var[idx]) / (bah_alt - alt[idx]) * (alt[idx+1:] - alt[idx]) +
							np.random.normal(0.0, noise_strength, n_alt_new-idx-1))

					else:
						# Or use mean T gradient of highest 20 measurements and continue with this gradient:
						# compute mean T gradient of highest 20 measurements:
						mean_T_grad = np.mean(np.asarray([(new_var[idx-19:idx+1] - new_var[idx-20:idx]) / (alt[idx-19:idx+1] - alt[idx-20:idx])]))
						extra_T = 288.15 + mean_T_grad*alt

						new_var[idx+1:] = extra_T[idx+1:] - (extra_T[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)

						if np.abs(new_var[-1] - bah_T) > 5:    # then the deviation from BAHAMAS T is too great and we use the ICAO std atmosphere
							new_var[idx+1:] = ICAO_standard_T[idx+1:] - (ICAO_standard_T[idx] - new_var[idx])


					new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'P':
					# Pressure: use hydrostatic eq. with scale height H = R <T> / g0, R = 287 J kg^-1 K^-1, g0 = 9.81 m s^-2, using the vertican mean temperature <T>:
					# p(z) = p0 exp( -z / H) (Holton, p.21); + noise (Hock and Franklin 1999)
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < ceiling/3:
						print("Insufficient amount of measurements for pressure extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no pressure measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					# MAKE SURE THAT mean TEMPERATURE CAPTURES THE ACTUAL MEAN TEMPERATURE!!
					if np.count_nonzero(~np.isnan(new_dict['T'][:])) / float(n_alt_new) <= 0.75: 		# in this case you may expect that a mean temperature would yield
																									# a bad representation of the true scale height. 0.75 was chosen arbitrarily.
						T_icao_0 = 12*np.cos(4*np.pi*np.nanmean(new_dict['lat'][:])/360) + 288.15	# strongly simplified meridional surface temperature structure
						H = 287 * np.mean(288.15 - 0.0065*alt) / 9.81		# using the ICAO standard atmosphere to compute the mean temperature

						print("Warning: Because insufficient non-nan temperature values were given for launch " +
							launch_time + ", '" + str(np.mean(288.15 - 0.0065*alt)) +
							" K' was assumed to be the mean temperature for hydrostatic pressure calculation. Can possibly be avoided if the temperature is extrapolated before the pressure.\n")

					else:
						H = 287 * np.nanmean(new_dict['T'][:]) / 9.81		# scale height

						# find index of lowest non nan pressure measurement:
						l_idx = np.argwhere(~np.isnan(new_var[:]))[0]
						p_ref = new_var[l_idx]		# in Pa
						alt_ref = alt[l_idx]
						p_hydrostat = p_ref * np.exp(-(alt - alt_ref) / H)		# in Pa

						new_var[idx+1:] = p_hydrostat[idx+1:] - (p_hydrostat[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'u_wind' or key == 'v_wind':
					# Wind: idea: fill nan values with the mean wind gradient of the highest 20 (non-nan)measurents. It will only be extrapolated if the the last non-nan entry
					# is higher than 0.80*ceiling:
					# other idea: just keep the same wind value
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.8*ceiling:
						print("Insufficient amount of measurements for wind extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no wind measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						extra_speed_length = 20		# amount of indices used for wind speed gradient calculation

						# # # for k in range(idx, n_alt_new):
							# # # new_var[n,k] = new_var[n,idx] + (k-idx)*(new_var[n,idx] - new_var[n,idx-extra_speed_length]) / (extra_speed_length + (k-idx))

						# alternative: just use the latest value for higher altitudes:
						new_var[idx+1:] = new_var[idx]
						new_var[idx+1:] = new_var[idx+1:] + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'RH':
					# Relative humidity (RH): Idea: fill nan values with the mean RH of the highest 10 measurements but only if the highest measurement exceeds or is equal to 0.90*ceiling!
					# A greater range probably doesn't make sense due to the high variability of relative humidity.
					# Other idea: Marek's suggestion: use rh = 0 for greater altitudes. I have decided not to use rh = 0 % because 0 % relative humidity will hardly ever be measured
					# Other idea: Linearly interpolate to the BAHAMAS value
					noise_strength = 0/2		# percent

					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.65*ceiling:
						print("Insufficient amount of measurements for relative humidity extrapolation at the top of the dropsonde grid (" + launch_time +
						"). There are no rel. hum. measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						# # # new_var[n,idx+1:] = np.mean(new_var[n,idx-9:idx+1]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_var[idx+1:] = (new_var[idx] + (bah_RH - new_var[idx]) / (bah_alt - alt[idx]) * (alt[idx+1:] - alt[idx]) +
							np.random.normal(0.0, noise_strength, n_alt_new - idx-1))
						# # # new_var[idx+1:] = 1.5 + np.random.normal(0.0, noise_strength, n_alt_new - idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag
						new_var[np.argwhere(new_var[:] < 0)] = 0.0



				# # # # # # More variables may be added here, if desired.


				new_dict[key] = new_var

			return new_dict, new_ipflag_dict, obs_height


		def std_extrapol(old_dict, ill_keys, old_ipflag_dict=dict()):
		# Will extrapolate some atmospheric variables to the ceiling of the dropsonde; old_ipflag will be updated.
		# Needs the old variable, the interpolation flag (should've been generated by fill_gaps()), the key and height levels as INPUT

			new_dict = old_dict
			n_alt = len(new_dict['Z'])

			new_ipflag_dict = old_ipflag_dict

			# to get the obs_height: find highest ... (?)
			if np.isnan(new_dict['reference_alt']):
				# select the highest non nan index of T or P. 
				highest_nonnan_Z = np.argwhere(~np.isnan(new_dict['Z']))[-1]
				obs_height = np.floor(new_dict['Z'][highest_nonnan_Z[0]]/100)*100
			else:	# use the reference_alt
				obs_height = (np.floor(new_dict['reference_alt']/100)*100)[0]
			

			ceiling = obs_height	# last entry of altitude
			if ceiling > 15000:
				print("Ceiling appears to be > 15000 m. Aborted extrapolation to dropsonde ceiling because the tropopause may intervene.\n")
				return new_dict, new_ipflag_dict


			# any value above obs_height will be deleted: So if e.g. Z has got values above obs_height, delete them:
			# find the first index that overshoots obs_height:
			with warnings.catch_warnings():
				warnings.simplefilter("ignore")
				overshoot = np.argwhere(new_dict['Z'] >= obs_height)
			if len(overshoot) > 0:
				overshoot = overshoot[0][0] + 1

				for key in new_dict.keys():
					if key in ['trajectory', 'fillValues', 'ipflag']:	# skip these ones ... it s not interesting anyway
						continue

					if new_dict[key].ndim > 0:		# otherwise: error when using len()

						if len(new_dict[key]) == n_alt:
							new_dict[key] = new_dict[key][:overshoot]	# limit variable to obs_height
							if key in ill_keys or key == 'Z':
								new_ipflag_dict[key] = new_ipflag_dict[key][:overshoot]


			# at the end of 'Z' there may still be nans -> so we don't know to which altitude meteorological variables belong to in this region:
			# therefore: delete it and replace by extrapolation:
			n_alt = len(new_dict['Z'])
			last_nonnan_alt = np.argwhere(~np.isnan(new_dict['Z']))[-1][0]
			if ceiling - new_dict['Z'][last_nonnan_alt] > 1000:
				print("WARNING: highest GPS altitude measurement is at least 1000 m below the aircraft. Extrapolation may be erroneous.\n")

			for key in new_dict.keys():
				if key in ['trajectory', 'fillValues', 'ipflag']:	# skip these ones ... it s not interesting anyway
					continue

				if new_dict[key].ndim > 0:		# otherwise: error when using len()

					if len(new_dict[key]) == n_alt:
						new_dict[key] = new_dict[key][:last_nonnan_alt+1]	# limit variable to obs_height
						if key in ill_keys or key == 'Z':
							new_ipflag_dict[key] = new_ipflag_dict[key][:last_nonnan_alt+1]


			# extend the old height grid up to the ceiling if the distance is greater than 10 meters:
			alt = new_dict['Z']
			n_alt = len(alt)
			alt = np.append(alt, np.arange(alt[np.argwhere(~np.isnan(alt))[-1]]+10, ceiling+11, 10))
			n_alt_new = len(alt)

			# update the altitude variable in the dictionary: & update ipflag for gpsalt:
			new_dict['Z'] = alt
			new_ipflag_dict['Z'] = np.append(new_ipflag_dict['Z'], np.ones((n_alt_new - n_alt,)))


			np.random.seed(42)		# needed for the added noise
			launch_time = datetime.datetime.utcfromtimestamp(new_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S") # for printing

			for key in ill_keys:

				new_var = new_dict[key]
				# must be expanded to the new height grid (up to the new ceiling)
				new_var = np.append(new_var, np.full((n_alt_new - n_alt,), np.nan), axis=0) # append nans at the top of the profile

				if not new_ipflag_dict: # in case fill_gaps(...) wasn't called before this one, it's assumed that nothing has been interpolated yet.
					new_ipflag_dict[key] = np.zeros(new_var.shape)

				else: # new_ipflag also has to be extended to the new hgt grid:
					new_ipflag_dict[key] = np.append(new_ipflag_dict[key], np.zeros((n_alt_new - n_alt,)), axis=0)



				if key == 'T':
					# Temperature: If dropsondes with measurements (ipflag = 0) from that day exist, estimate their average T gradient.
					# If it clearly deviates from the standard atmospheric T gradient, then use a modified ICAO_standard atmosphere which
					# has an adapted T gradient. Otherwise:
					# Assume standard atmosphere (shifted accordingly to avoid a jump between the last known value and the extrapolation;
					# + noise (whose strength is according to instrument noise (from Hock and Franklin 1999)). ICAO standard atmosphere taken from:
					# https://www.dwd.de/DE/service/lexikon/begriffe/S/Standardatmosphaere_pdf.pdf?__blob=publicationFile&v=3
					ICAO_standard_T = 288.15 - 0.0065*alt
					noise_strength = 0/2


					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.6*ceiling:
						print("Insufficient amount of measurements for temperature extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no temperature measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var			# then just overwrite the dictionary entry with the nonedited (but extended) variable
						continue

					if alt[idx] < 0.85*ceiling: # then use standard atmosphere (ICAO):
						new_var[idx+1:] = ICAO_standard_T[idx+1:] + (new_var[idx] - ICAO_standard_T[idx])

					else:
						# Or use mean T gradient of highest 20 measurements and continue with this gradient:
						# compute mean T gradient of highest 20 measurements:
						mean_T_grad = np.mean(np.asarray([(new_var[idx-19:idx+1] - new_var[idx-20:idx]) / (alt[idx-19:idx+1] - alt[idx-20:idx])]))
						T_continued = 288.15 + mean_T_grad*alt

						new_var[idx+1:] = T_continued[idx+1:] - (T_continued[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)


					new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'P':
					# Pressure: use hydrostatic eq. with scale height H = R <T> / g0, R = 287 J kg^-1 K^-1, g0 = 9.81 m s^-2, using the vertican mean temperature <T>:
					# p(z) = p0 exp( -z / H) (Holton, p.21); + noise (Hock and Franklin 1999)
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < ceiling/3:
						print("Insufficient amount of measurements for pressure extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no pressure measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					# MAKE SURE THAT mean TEMPERATURE CAPTURES THE ACTUAL MEAN TEMPERATURE!!
					if np.count_nonzero(~np.isnan(new_dict['T'][:])) / float(n_alt_new) <= 0.75: 		# in this case you may expect that a mean temperature would yield
																									# a bad representation of the true scale height. 0.75 was chosen arbitrarily.
						T_icao_0 = 12*np.cos(4*np.pi*np.nanmean(new_dict['lat'][:])/360) + 288.15	# strongly simplified meridional surface temperature structure
						H = 287 * np.mean(288.15 - 0.0065*alt) / 9.81		# using the ICAO standard atmosphere to compute the mean temperature

						print("Warning: Because insufficient non-nan temperature values were given for launch " +
							launch_time + ", '" + str(np.mean(288.15 - 0.0065*alt)) +
							" K' was assumed to be the mean temperature for hydrostatic pressure calculation. Can possibly be avoided if the temperature is extrapolated before the pressure.\n")

					else:
						H = 287 * np.nanmean(new_dict['T'][:]) / 9.81		# scale height

						# find index of lowest non nan pressure measurement:
						l_idx = np.argwhere(~np.isnan(new_var[:]))[0]
						p_ref = new_var[l_idx]		# in Pa
						alt_ref = alt[l_idx]
						p_hydrostat = p_ref * np.exp(-(alt - alt_ref) / H)		# in Pa

						new_var[idx+1:] = p_hydrostat[idx+1:] - (p_hydrostat[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'u_wind' or key == 'v_wind':
					# Wind: idea: fill nan values with the mean wind gradient of the highest 20 (non-nan)measurents. It will only be extrapolated if the the last non-nan entry
					# is higher than 0.80*ceiling:
					# other idea: just keep the same wind value
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.8*ceiling:
						print("Insufficient amount of measurements for wind extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no wind measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						# # # extra_speed_length = 20		# amount of indices used for wind speed gradient calculation

						# # # for k in range(idx, n_alt_new):
							# # # new_var[n,k] = new_var[n,idx] + (k-idx)*(new_var[n,idx] - new_var[n,idx-extra_speed_length]) / (extra_speed_length + (k-idx))

						# alternative: just use the latest value for higher altitudes:
						new_var[idx+1:] = new_var[idx]
						new_var[idx+1:] = new_var[idx+1:] + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'RH':
					# Relative humidity (RH): Idea: fill nan values with the mean RH of the highest 10 measurements but only if the highest measurement exceeds or is equal to 0.90*ceiling!
					# A greater range probably doesn't make sense due to the high variability of relative humidity.
					# Other idea: Marek's suggestion: use RH approx 0 for greater altitudes. I have decided not to use rh = 0 % because 0 % relative humidity will hardly ever be measured
					# Other idea: Linearly interpolate to the BAHAMAS value
					noise_strength = 0/2		# percent

					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.65*ceiling:
						print("Insufficient amount of measurements for relative humidity extrapolation at the top of the dropsonde grid (" + launch_time +
						"). There are no rel. hum. measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						# # # new_var[idx+1:] = np.mean(new_var[idx-9:idx+1]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_var[idx+1:] = 1.5 + np.random.normal(0.0, noise_strength, n_alt_new - idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag
						new_var[np.argwhere(new_var[:] < 0)] = 0.0



				# # # # # # More variables may be added here, if desired.


				new_dict[key] = new_var

			return new_dict, new_ipflag_dict, obs_height


		def regridding(new_dict, obs_height, ill_keys, resolution=10):
			'''
				Regridding variables specified in ill_keys to a uniform grid 
				(from the surface up to obs_height) with a user-defined 
				resolution (in meters, default=10).
			'''

			new_alt = np.arange(0, obs_height+1, 10)

			for key in ill_keys:
				new_dict[key] = np.interp(new_alt, new_dict['Z'], new_dict[key])

			new_dict['Z'] = new_alt

			return new_dict


		def repair_surface(old_dict, ill_keys, old_ipflag_dict=dict()):
		# Filling nan values at the surface if the gap to the surface isn't too large (e.g. measurements below 150 m must exist (roughly 10-15 seconds before splash).
			new_dict = old_dict
			alt = old_dict['Z']
			n_alt = len(alt)
			new_ipflag_dict = old_ipflag_dict
			launch_time = datetime.datetime.utcfromtimestamp(new_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S")

			lim = 200		# if there are no measurements below this altitude then the extrapolation at the surface won't be performed

			if ill_keys == ['Z']:
				threshold_list = [ill_keys, [200], ['m']]
			else:
				threshold_list = [ill_keys, [5.0, 4000.0, 50.0, 0.1, 0.1, 5.0, 5.0, 1.0],
					['K', 'hPa', '%', 'deg', 'deg', 'm/s', 'm/s', 'm/s']]		# used to check if surface value deviates siginificantly from lowest measurement

			for key in ill_keys:

				new_var = new_dict[key]

				if not new_ipflag_dict: # in case fill_gaps(...) wasn't called before this one, it's assumed that nothing has been interpolated yet.
					new_ipflag_dict[key] = np.zeros(new_var.shape)


				# find the first non-nan entry
				idx = np.argwhere(~np.isnan(new_var[:]))[0][0]

				if alt[idx] < lim:
					sfc_gap = np.arange(0,idx)

					if len(sfc_gap) == 0:
						continue
					else:

						# create mean gradient of the variable of 10 measurements above the lowest measurement, or, if grid is too coarse, take 200-400 m average:
						if alt[idx+10] > 2*lim: # default: if alt[idx+10] > 400: # take lowest measurement to 400 m mean gradient
							# find index closest to 400 m:
							idx2 = np.argmin(np.abs(alt - 2*lim))

						else: # take mean grad. of 10 measurem. above lowest measurement:
							idx2 = idx+10

						mean_grad = np.mean([new_var[j+1] - new_var[j] for j in range(idx,idx2)])	# theoretically, should never be nan because fill_gaps
																										# should've fixed the holes between the first and last measurement
						for j in sfc_gap:
							new_var[idx-j-1] = new_var[idx] - mean_grad*(j+1)


						# check if sfc value not too far off the lowest measurement:
						if key == 'RH':
							if np.any(new_var[sfc_gap] < 0):
								new_var[sfc_gap] = 0
								print("Caution, '" + key + "' surface repair resulted in negative values. Manually set the missing values at the ground to 0 for launch "
									+ launch_time + ".\n")

							elif np.any(new_var[sfc_gap] > 100):
								new_var[sfc_gap] = 100
								print("Caution, '" + key + "' surface repair resulted in >100 %. Manually set the missing values at the ground to 100 for launch "
									+ launch_time + ".\n")

						threshold = threshold_list[1][threshold_list[0].index(key)]
						si_unit = threshold_list[2][threshold_list[0].index(key)]
						if np.abs(new_var[0] - new_var[idx]) > threshold:
							print("Caution, '" + key + "' surface value deviates more than " + str(threshold) + " " + si_unit + " from the lowest measurement (launch "
								+ launch_time + ").\n")


						new_ipflag_dict[key][sfc_gap] = 1

				else:
					print("No measurements below " + str(lim) + " m. Extrapolation of '" + key + "', launch " + launch_time +
						" would eventually lead to wrong assumptions at the surface. Therefore aborted.\n")
					continue

			return new_dict, new_ipflag_dict


		def mark_outliers(sonde_dict, ill_keys): # mark outliers: outliers defined when exceeding certain thresholds

			new_dict = sonde_dict

			# thresholds are defined by change of meteorol. variable with altitude: e.g. delta p / delta z
			thresholds = [0.065, 40, 2.5, 1, 1]	# in [K/m, Pa/m, %/m, ms^-1/m, ms^-1/m]
			dz = new_dict['Z'][1:] - new_dict['Z'][:-1]	# delta z

			for key in ill_keys:
				if key == 'lat' or key == 'lon':
					continue

				met_threshold = thresholds[ill_keys.index(key)]		# threshold for key

				d_met = new_dict[key][1:] - new_dict[key][:-1]		# change of meteorological variable 'key'

				with warnings.catch_warnings():
					warnings.simplefilter("ignore")
					exceed_idx = np.argwhere(np.abs(d_met / dz) >= met_threshold)
				new_dict[key][exceed_idx] = float('nan')

			return new_dict


		def plot_met_profile(sonde_dict, ill_keys, plot_path, plot_filename_base): # plots T profile and saves it in 'plot_path'

			units = ['K', 'Pa', '%']

			for key in ill_keys:	# plot each meteorological variable that has been modified:
				# Plotting after extrapolation:
				font_size = 14
				fig = plt.figure(figsize=(6,9))
				a1 = plt.axes()

				launch_date = datetime.datetime.utcfromtimestamp(sonde_dict['launch_time']).strftime("%Y%m%d_%H%M%S")
				a1.plot(sonde_dict[key], sonde_dict['Z'], linewidth=1.2, color=(0,0,0))

				titletext = r"Dropsonde " + key + " profile from EUREC4A campaign: " + launch_date
				plt.title(titletext, fontsize=font_size, wrap=True)
				a1.set_xlabel(key + " [" + units[ill_keys.index(key)] + "]", fontsize=font_size)
				a1.set_ylabel(r"Height [m]", fontsize=font_size)
				a1.grid(True, axis='x', which='both')
				a1.grid(True, axis='y', which='major')
				a1.set_ylim(bottom=0, top=sonde_dict['Z'][-1])

				if key == 'tdry':
					a1.set_xlim(left=240, right=305)
				elif key == 'pres':
					a1.set_xlim(left=10000, right=105000)
				elif key == 'rh':
					a1.set_xlim(left=0, right=100)

				plt.savefig(plot_path + plot_filename_base + "_" + key + ".png") #, dpi=250, bbox_inches='tight'
				plt.close()


		def saveExpolSondeAsNC(sonde_dict, out_filename): # saves the sonde_dict as an nc file named out_filename into data_out_path

			new_nc = nc.Dataset(out_filename, "w", format="NETCDF4")

			# create dimensions:
			n_alt = len(sonde_dict['Z'])
			new_nc.createDimension("alt", n_alt)

			# convert


			# Global attributes:
			new_nc.description = """EUREC4A campaign RAW dropsondes. Extrapolated when enough measurements were given. More information required? Contact author (listed below).
				'alt_old' represents the dimension of height levels of the uninterpolated / extrapolated variables. alt represents the height levels of the extrapolated variables."""
			new_nc.history = "Created: " + datetime.datetime.utcnow().strftime("%Y-%m-%d %H:%M:%S")
			new_nc.author = "Andreas Walbroel (Mail: a.walbroel@uni-koeln.de)"


			# create variables, set attributes and save values into the created variables:
			# attributes = [[units], [description]]
			attributes = [['m', 'm s^-1', 'degree East', 'seconds since 1970-01-01 00:00:00 UTC', 'degree East', '%', 'm s^-1',
			'Pa', 'Pa', 'degree North', 'm', '%', 'K', 'seconds since 1970-01-01 00:00:00 UTC',
			'K', 'degree North'],
			['gps reported altitude above mean sea level', 'meridional wind (northwards is > 0)', 'reference longitude',
			'reference time', 'longitude', 'reference relative humidity', 'zonal wind (eastwards is > 0)', 'air pressure',
			'reference pressure', 'latitude', 'reference altitude (aircraft height)', ' relative humidity', 'reference temperature',
			'sonde launch time', 'temperature', 'reference latitude']]

			sonde_keys = ['Z', 'v_wind', 'reference_lon', 'reference_time', 'lon', 'reference_rh', 'u_wind', 'P', 'reference_pres',
				'lat', 'reference_alt', 'RH', 'reference_tdry', 'launch_time', 'T', 'reference_lat']

			for at_idx, key in enumerate(sonde_keys):

				# find out about this variable's dimension:
				if key == 'launch_time':
					var = np.asarray([sonde_dict[key]])

				else:
					var = sonde_dict[key]


				if var.ndim == 1:
					if len(var) == n_alt:
						nc_var = new_nc.createVariable(key, "f8", ("alt"))

					elif len(var) == 1:
						nc_var = new_nc.createVariable(key, "f8")

					else:
						print("""Hmmm.... something went wrong. Actually all variables should either have the dimension of the
							altitude or length = 1. But """ + key + " has got another shape. Happy debugging!")

					nc_var[:] = var


				nc_var.units = attributes[0][at_idx]
				nc_var.description  = attributes[1][at_idx]

			new_nc.close()
			print("Output: ", out_filename)



		############################################################################################



		# Dropsonde quality control: if there are small gaps of measurements between launch altitude and surface, they will be filled via interpolation.
		# Additionally, the sondes will be extrapolated to a certain altitude.

		# Reading in dropsonde files:
		HALO_sondes_NC = sorted(glob.glob(path_raw_sondes + "*PQC.nc"))
		if path_BAH_data:
			BAH_files_NC = sorted(glob.glob(path_BAH_data + "bahamas*.nc"))
		else:
			BAH_files_NC = []


		tot_failure_warnings = 0			# counts the amount of times a critical variable has got no measurements
		tot_sonde_stuck = 0
		failed_sondes = []
		stuck_sondes = []

		print('Found %d sondes.' % len(HALO_sondes_NC))

		for sonde_nc in HALO_sondes_NC:

			sonde_dict = readrawNCraw(sonde_nc)
			launch_date = datetime.datetime.utcfromtimestamp(sonde_dict['launch_time']).strftime("%Y-%m-%d")	# time delta required
			dropsonde_date = (datetime.datetime.strptime(launch_date, "%Y-%m-%d")).strftime("%Y%m%d")	# date displayed in the filename ... comfy way to find the right BAHAMAS data for std_extrapol
			print("########## Day: " + datetime.datetime.utcfromtimestamp(sonde_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S") + " ##########\n")
			print("Input: ", sonde_nc)


			# add another condition that checks if e.g. nearly no measurements exist at all (for T, P and RH):
			if np.any([np.count_nonzero(~np.isnan(sonde_dict['T'])) < 0.1*len(sonde_dict['T']),
				np.count_nonzero(~np.isnan(sonde_dict['P'])) < 0.1*len(sonde_dict['P']),
				np.count_nonzero(~np.isnan(sonde_dict['RH'])) < 0.1*len(sonde_dict['RH']),
				np.count_nonzero(~np.isnan(sonde_dict['lat'])) < 0.05*len(sonde_dict['lat']),
				np.count_nonzero(~np.isnan(sonde_dict['lon'])) < 0.05*len(sonde_dict['lon']),
				np.count_nonzero(~np.isnan(sonde_dict['u_wind'])) < 0.05*len(sonde_dict['u_wind']),
				np.count_nonzero(~np.isnan(sonde_dict['v_wind'])) < 0.05*len(sonde_dict['v_wind'])]):
				tot_failure_warnings = tot_failure_warnings + 1
				failed_sondes.append(sonde_nc)
				print("One PAMTRA-critical variable measurement failed. Skipping this dropsonde.\n")
				continue

			# add jet another condition that checks if the sonde got stuck mid air (gps alt. values don't really decrease with time):
			if not np.any(sonde_dict['Z'][~np.isnan(sonde_dict['Z'])] < 1500):	# then I assume that the whole launch was doomed
				print("Sonde got stuck in mid air. 'gpsalt' doesn't seem to include any values < 1500 m.\n")
				stuck_sondes.append(sonde_nc)
				tot_sonde_stuck = tot_sonde_stuck + 1
				continue


			# it's known that the altitude axis (gpsalt) does have some broken values (e.g. sudden jumps or exceeding the aircraft altitude):
			# but before we deal with those, the 'normal' gaps in gpsalt shall be fixed


			# subsequent variables will be cured from holey nan value disease...:
			# pressure, temperature, relhum, wind (u & v & w), lat, lon.
			# Other variables, of which you can expect a linear interpolated over gaps to be applicable, may be added.
			ill_keys = ['T', 'P', 'RH', 'lat', 'lon', 'u_wind', 'v_wind']


			sonde_ipflag = dict()		# will contain the interpolation flags for interpolated nan values in the middle of the drop

			# before filliing the gaps, the altitude axis must be fixed:
			sonde_dict['Z'], sonde_ipflag['Z'] = fill_gaps(sonde_dict['Z'])
			for key in ill_keys:
				sonde_dict[key], sonde_ipflag[key] = fill_gaps(sonde_dict[key])		# altitude must be passed to check for dimensions of the to-be-cured variable...

			sonde_dict['ipflag'] = sonde_ipflag



			# the raw dropsonde files show an altitude variable with increases from [0 to -1] in general: but probably due to gps tracking,
			# the altitude decreases at the "top": this must be filtered out:
			# find highest index where altitude[idx+1] - altitude[idx] > 0:
			# CATCH WARNINGS: This is a warning that appears because nans still lie within sonde_dict['Z'] and a comparison if nan
			# is greater than 0 is not possible
			with warnings.catch_warnings():
				warnings.simplefilter("ignore")
				altitude_stop = np.argwhere(sonde_dict['Z'][1:] - sonde_dict['Z'][0:-1] > 0)[-1][0] + 2
				# +2 because it's used as indexing [... : altitude_stop] => +1 and because the array size had been reduced by 1 during argwhere(...)


			# if the lowest non nan value of the altitude coordinate is < 0: cut the rest off:
			# find lowest non nan value of altitude:
			lowest = np.argwhere(~np.isnan(sonde_dict['Z'][:]))[0][0]

			# then cut each variable at this altitude index:
			ndata = len(sonde_dict['Z'])

			if sonde_dict['Z'][lowest] < 0:
				for key in sonde_dict.keys():
					if key in ['trajectory', 'fillValues', 'ipflag']:	# skip these ones ... it s not interesting anyway
						continue

					if sonde_dict[key].ndim > 0:		# otherwise: error when using len()

						if len(sonde_dict[key]) == ndata:
							sonde_dict[key] = sonde_dict[key][lowest:altitude_stop]
							if key in ill_keys or key == 'Z':
								sonde_dict['ipflag'][key] = sonde_dict['ipflag'][key][lowest:altitude_stop]

			else:
				for key in sonde_dict.keys():
					if key in ['trajectory', 'fillValues', 'ipflag']:	# skip these ones ... it s not interesting anyway
						continue

					if sonde_dict[key].ndim > 0:		# otherwise: error when using len()

						if len(sonde_dict[key]) == ndata:
							sonde_dict[key] = sonde_dict[key][:altitude_stop]
							if key in ill_keys or key == 'Z':
								sonde_dict['ipflag'][key] = sonde_dict['ipflag'][key][:altitude_stop]


			# the altitude index may be a bit broken... needs to be fixed. mark them as nan and let it run through fill_gaps again:
			dz = sonde_dict['Z'][1:] - sonde_dict['Z'][:-1]		# dz[i] = z[i+1] - z[i]

			for k in range(len(dz)):
				if (dz[k] < 0) or (dz[k] > 15):		# this filters too big jumps in the altitude coordinate
					sonde_dict['Z'][k+1] = float('nan')

			sonde_dict['Z'], sonde_ipflag['Z'] = fill_gaps(sonde_dict['Z'])


			# perform surface repair for altitude coordinate:
			sonde_dict, sonde_dict['ipflag'] = repair_surface(sonde_dict, ['Z'], sonde_dict['ipflag'])
			# for some reasons, 'alt' is sort of unused but an assigned key in the dictionary. So ... we can also just set it to gpsalt:
			sonde_dict['alt'] = sonde_dict['Z']

			# now we still need to handle the nan values at the surface: If there are no non-nan values in the lowest 5 % of the variable
			# -> don't interpolate because the assumption would eventually lead to senseless surface values:
			sonde_dict, sonde_dict['ipflag'] = repair_surface(sonde_dict, ill_keys, sonde_dict['ipflag'])


			# Extrapolating the ill_keys to the ceiling of the dropsondes (e.g. below aircraft altitude):
			# CAUTION: it is expected that the dropsondes start BELOW THE TROPOPAUSE!
			# It is advisable to do the TEMPERATURE EXTRAPOLATION FIRST because it improves the pressure extrapolation.
			# Need bahamas file for extrapolation limit:
			if BAH_files_NC:
				bah_filename = [bah_file for idx, bah_file in enumerate(BAH_files_NC) if dropsonde_date in bah_file]
				sonde_dict, sonde_dict['ipflag'], obs_height = std_extrapol_BAH(sonde_dict, ill_keys, bah_filename, sonde_dict['ipflag'])

			else:
				sonde_dict, sonde_dict['ipflag'], obs_height = std_extrapol(sonde_dict, ill_keys, sonde_dict['ipflag'])

			# Regridding to a uniform vertical grid with a user-specified resolution:
			sonde_dict = regridding(sonde_dict, obs_height, ill_keys, 10)


			# find outliers and mark them (as nan): afterwards fill them again
			sonde_dict = mark_outliers(sonde_dict, ['T', 'P', 'RH', 'u_wind', 'v_wind'])
			for key in ill_keys:
				sonde_dict[key], sonde_ipflag[key] = fill_gaps(sonde_dict[key])


			# Save the extrapolated sonde dictionary to a new nc file:
			data_out_path_halo_dir = os.path.dirname(data_out_path_halo)
			if not os.path.exists(data_out_path_halo_dir):
				os.makedirs(data_out_path_halo_dir)
			out_filename = os.path.basename(sonde_nc)		# removes the path in the string HALO_sondes_NC[m] so that the filename remains
			out_filename = os.path.join(data_out_path_halo, out_filename[0:-3] + "_RAW_v01.nc")

			saveExpolSondeAsNC(sonde_dict, out_filename)


		# # # # # # save the lists containing the broken and stuck files:
		# # # # # np.save(data_out_path_halo + "tot_failure_warnings.npy", np.asarray(tot_failure_warnings))
		# # # # # np.save(data_out_path_halo + "tot_sonde_stuck.npy", np.asarray(tot_sonde_stuck))


###############################################################################################
###############################################################################################
###############################################################################################
###############################################################################################


	def run_dropsonde_gap_filler_joanne(path_raw_sondes, data_out_path_halo,
		path_BAH_data=''):

		"""
		Parameters
		----------
		path_raw_sondes : str
			Path of raw dropsonde data.
		data_out_path_halo : str
			Path of repaired dropsonde data (gaps filled).
		path_BAH_data : str, optional
			Path of BAHAMAS data in unified netCDF files.
		"""


		############################################################################################
		# FUNCTIONS

		def fill_gaps(old_var):
		# old variable gets linearly interpolated for each sonde launch. The function is ignoring nan values at the surface and above
		# the launch altitude.

			new_var = copy.deepcopy(old_var)
			# create flag variable indicating if an entry of old_var has been changed: if = 0: not interpol.
			interp_flag = np.zeros(old_var.shape)


			# identify regions of nan values in the middle of the drop. Extrapolation will be handled in another function.
			# identify the highest non-nan entry so we can cut the values above that highest entry:
			# identify the lowest non-nan entry for similar reasons:
			non_nan_idx = [idx for idx, x in np.ndenumerate(old_var) if (not np.isnan(x))]
			non_nan_idx = np.where(~np.isnan(old_var))[0]
			limits = np.array([non_nan_idx[0], non_nan_idx[-1]])

			temp_var = copy.deepcopy(old_var)
			temp_var = temp_var[limits[0]:limits[1]+1]		# will be the variable where the gaps are filled

			interp_flag_temp = np.zeros(temp_var.shape)

			# identify mid-drop-nan-values: need values after and before the nan:
			nan_idx = np.argwhere(np.isnan(temp_var))
			interp_flag_temp[nan_idx] = 1

			if nan_idx.size == 0:
				return new_var, interp_flag

			else: # correct nan values: find the hole size via subtraction of subsequent indices
				hole_size = np.zeros((len(nan_idx)+1,)).astype(int)
				# hole_size = np.zeros(nan_idx.shape)			# old version
				k = 0		# index to address a hole ('hole number')

				for m in range(0, len(temp_var)-1):

					if not np.isnan(temp_var[m+1] - temp_var[m]):
						hole_size[k] = 0

					elif np.isnan(temp_var[m+1] - temp_var[m]): # k shall only be incremented if an END of a hole has been identified:
						if len(nan_idx) == 1: 	# must be handled seperately in case that merely one nan value exists in temp_var
							hole_size[k] = 1
							break

						else:
							if (not np.isnan(temp_var[m+1])) & (np.isnan(temp_var[m])): # END of a hole
								k = k + 1
								continue
							hole_size[k] = hole_size[k] + 1		# k won't be incremented until m finds another non-nan value

					else:
						print("\n Something unexpected happened when trying to find the nan values in the middle of the dropsonde launch... Contact 'a.walbroel@uni-koeln.de'. \n")

				# holes have been identified: edit the FIRST hole (editing depends on the size of the hole...)
				c = 0 		# dummy variable needed for the right jumps in hole_size and nan_idx. c is used to address nan_idx and therefore new_var...

				# meanwhile 'd' just runs through the array hole_size:
				for d in range(0, len(hole_size)):
					for L in range(0, hole_size[d]):		# range(0, 1): L = 0
						temp_var[nan_idx[c] + L] = temp_var[nan_idx[c] - 1] + (L + 1)*(temp_var[int(nan_idx[c] + hole_size[d])] - temp_var[nan_idx[c]-1]) / (hole_size[d] + 1)

					c = c + int(hole_size[d])
					if c > len(hole_size)-1:
						break


			# overwrite the possibly holey section:
			new_var[limits[0]:limits[1]+1] = temp_var
			# update interp_flag
			interp_flag[limits[0]:limits[1]+1] = interp_flag_temp

			return new_var, interp_flag


		def std_extrapol_BAH(old_dict, ill_keys, bah_filename, old_ipflag_dict=dict()):
		# Will extrapolate some atmospheric variables to the ceiling of the dropsonde; old_ipflag will be updated.
		# Needs the old variable, the interpolation flag (should've been generated by fill_gaps()), the key and height levels as INPUT

			new_dict = old_dict
			n_alt = len(new_dict['Z'])

			new_ipflag_dict = old_ipflag_dict


			# Need BAHAMAS information to set the new ceiling:
			# Import altitude and time data from BAHAMAS:
			bah_keys = ['time', 'altitude', 'ta', 'p', 'rh']
			bah_dict = import_BAHAMAS_unified(bah_filename[0], bah_keys)
			bah_dict['time'] = np.rint(bah_dict['time']).astype(float) # must be done to avoid small fractions of seconds

			# to get the obs_height: average BAHAMAS altitude over +/- 10 seconds around launch_time:
			# find time index of the sonde launches:
			timestamp = new_dict['launch_time']
			bah_launch_idx = np.asarray([np.argwhere(bah_dict['time'] == timestamp)]).flatten()		# had some dimensions too many -> flattened
			drop_alt = np.floor(np.asarray([np.mean(bah_dict['altitude'][i-10:i+10]) for i in bah_launch_idx])/100)*100
			obs_height = np.max(np.unique(drop_alt))		# some values are repeated ... omit them and find the max value: this will be used for top of extrapolation!

			# BAHAMAS temperature, pressure, relative humidity at launch time:
			bah_T = np.nanmean(bah_dict['ta'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])
			bah_P = np.nanmean(bah_dict['p'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])
			bah_RH = np.nanmean(bah_dict['rh'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])
			bah_alt = np.nanmean(bah_dict['altitude'][bah_launch_idx[0]-10:bah_launch_idx[0]+10])


			ceiling = obs_height	# last entry of altitude
			if ceiling > 15000:
				print("Ceiling appears to be > 15000 m. Aborted extrapolation to dropsonde ceiling because the tropopause may intervene.\n")
				return new_dict, new_ipflag_dict


			# any value above obs_height will be deleted: So if e.g. Z has got values above obs_height, delete them:
			# find the first index that overshoots obs_height:
			with warnings.catch_warnings():
				warnings.simplefilter("ignore")
				overshoot = np.argwhere(new_dict['Z'] >= obs_height)
			ill_keys.append('Z')
			if len(overshoot) > 0:
				overshoot = overshoot[0][0] + 1

				for key in ill_keys:

					if new_dict[key].ndim > 0:		# otherwise: error when using len()

						if len(new_dict[key]) == n_alt:
							new_dict[key] = new_dict[key][:overshoot]	# limit variable to obs_height
							if key in ill_keys or key == 'Z':
								new_ipflag_dict[key] = new_ipflag_dict[key][:overshoot]


			# at the end of 'Z' there may still be nans -> so we don't know to which altitude meteorological variables belong to in this region:
			# therefore: delete it and replace by extrapolation:
			n_alt = len(new_dict['Z'])
			last_nonnan_alt = np.argwhere(~np.isnan(new_dict['Z']))[-1][0]
			if ceiling - new_dict['Z'][last_nonnan_alt] > 1000:
				print("WARNING: highest GPS altitude measurement is at least 1000 m below the aircraft. Extrapolation may be erroneous.\n")

			for key in ill_keys:

				if new_dict[key].ndim > 0:		# otherwise: error when using len()

					if len(new_dict[key]) == n_alt:
						new_dict[key] = new_dict[key][:last_nonnan_alt+1]	# limit variable to obs_height
						if key in ill_keys or key == 'Z':
							new_ipflag_dict[key] = new_ipflag_dict[key][:last_nonnan_alt+1]


			# extend the old height grid up to the ceiling if the distance is greater than 10 meters:
			alt = new_dict['Z']
			n_alt = len(alt)
			alt = np.append(alt, np.arange(alt[np.argwhere(~np.isnan(alt))[-1]]+10, ceiling+1, 10))
			n_alt_new = len(alt)

			# update the altitude variable in the dictionary: & update ipflag for gpsalt:
			new_dict['Z'] = alt
			new_ipflag_dict['Z'] = np.append(new_ipflag_dict['Z'], np.ones((n_alt_new - n_alt,)))


			np.random.seed(42)		# needed for the added noise
			launch_time = datetime.datetime.utcfromtimestamp(new_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S") # for printing

			ill_keys.pop()
			for key in ill_keys: # .pop() removes 'Z' from ill_keys

				new_var = new_dict[key]
				# must be expanded to the new height grid (up to the new ceiling)
				new_var = np.append(new_var, np.full((n_alt_new - n_alt,), np.nan), axis=0) # append nans at the top of the profile

				if not new_ipflag_dict: # in case fill_gaps(...) wasn't called before this one, it's assumed that nothing has been interpolated yet.
					new_ipflag_dict[key] = np.zeros(new_var.shape)

				else: # new_ipflag also has to be extended to the new hgt grid:
					new_ipflag_dict[key] = np.append(new_ipflag_dict[key], np.zeros((n_alt_new - n_alt,)), axis=0)



				if key == 'T':
					# If BAHAMAS Temperature measurement is available use it as target in case only the top 15 % of measurements
					# are missing. Otherwise:
					# Temperature: If dropsondes with measurements (ipflag = 0) from that day exist, estimate their average T gradient.
					# If the extrapolated dropsonde temperature then deviates from BAH T by more than 5 K, use the ICAO std atmosphere
					# as T gradient:
					# Standard atmosphere (shifted accordingly to avoid a jump between the last known value and the extrapolation).
					# ICAO standard atmosphere taken from:
					# https://www.dwd.de/DE/service/lexikon/begriffe/S/Standardatmosphaere_pdf.pdf?__blob=publicationFile&v=3
					ICAO_standard_T = 288.15 - 0.0065*alt
					noise_strength = 0/2


					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.6*ceiling:
						print("Insufficient amount of measurements for temperature extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no temperature measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var			# then just overwrite the dictionary entry with the nonedited (but extended) variable
						continue

					if alt[idx] < 0.85*ceiling: # then use BAHAMAS temperature as extrapolation target:
						new_var[idx+1:] = (new_var[idx] + (bah_T - new_var[idx]) / (bah_alt - alt[idx]) * (alt[idx+1:] - alt[idx]) +
							np.random.normal(0.0, noise_strength, n_alt_new-idx-1))

					else:
						# Or use mean T gradient of highest 20 measurements and continue with this gradient:
						# compute mean T gradient of highest 20 measurements:
						mean_T_grad = np.mean(np.asarray([(new_var[idx-19:idx+1] - new_var[idx-20:idx]) / (alt[idx-19:idx+1] - alt[idx-20:idx])]))
						extra_T = 288.15 + mean_T_grad*alt

						new_var[idx+1:] = extra_T[idx+1:] - (extra_T[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)

						if np.abs(new_var[-1] - bah_T) > 5:    # then the deviation from BAHAMAS T is too great and we use the ICAO std atmosphere
							new_var[idx+1:] = ICAO_standard_T[idx+1:] - (ICAO_standard_T[idx] - new_var[idx])


					new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'P':
					# Pressure: use hydrostatic eq. with scale height H = R <T> / g0, R = 287 J kg^-1 K^-1, g0 = 9.81 m s^-2, using the vertican mean temperature <T>:
					# p(z) = p0 exp( -z / H) (Holton, p.21); + noise (Hock and Franklin 1999)
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < ceiling/3:
						print("Insufficient amount of measurements for pressure extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no pressure measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					# MAKE SURE THAT mean TEMPERATURE CAPTURES THE ACTUAL MEAN TEMPERATURE!!
					if np.count_nonzero(~np.isnan(new_dict['T'][:])) / float(n_alt_new) <= 0.75: 		# in this case you may expect that a mean temperature would yield
																									# a bad representation of the true scale height. 0.75 was chosen arbitrarily.
						T_icao_0 = 12*np.cos(4*np.pi*np.nanmean(new_dict['lat'][:])/360) + 288.15	# strongly simplified meridional surface temperature structure
						H = 287 * np.mean(288.15 - 0.0065*alt) / 9.81		# using the ICAO standard atmosphere to compute the mean temperature

						print("Warning: Because insufficient non-nan temperature values were given for launch " +
							launch_time + ", '" + str(np.mean(288.15 - 0.0065*alt)) +
							" K' was assumed to be the mean temperature for hydrostatic pressure calculation. Can possibly be avoided if the temperature is extrapolated before the pressure.\n")

					else:
						H = 287 * np.nanmean(new_dict['T'][:]) / 9.81		# scale height

						# find index of lowest non nan pressure measurement:
						l_idx = np.argwhere(~np.isnan(new_var[:]))[0]
						p_ref = new_var[l_idx]		# in Pa
						alt_ref = alt[l_idx]
						p_hydrostat = p_ref * np.exp(-(alt - alt_ref) / H)		# in Pa

						new_var[idx+1:] = p_hydrostat[idx+1:] - (p_hydrostat[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'u_wind' or key == 'v_wind':
					# Wind: idea: fill nan values with the mean wind gradient of the highest 20 (non-nan)measurents. It will only be extrapolated if the the last non-nan entry
					# is higher than 0.80*ceiling:
					# other idea: just keep the same wind value
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.8*ceiling:
						print("Insufficient amount of measurements for wind extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no wind measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						extra_speed_length = 20		# amount of indices used for wind speed gradient calculation

						# # # for k in range(idx, n_alt_new):
							# # # new_var[n,k] = new_var[n,idx] + (k-idx)*(new_var[n,idx] - new_var[n,idx-extra_speed_length]) / (extra_speed_length + (k-idx))

						# alternative: just use the latest value for higher altitudes:
						new_var[idx+1:] = new_var[idx]
						new_var[idx+1:] = new_var[idx+1:] + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'RH':
					# Relative humidity (RH): Idea: fill nan values with the mean RH of the highest 10 measurements but only if the highest measurement exceeds or is equal to 0.90*ceiling!
					# A greater range probably doesn't make sense due to the high variability of relative humidity.
					# Other idea: Marek's suggestion: use rh = 0 for greater altitudes. I have decided not to use rh = 0 % because 0 % relative humidity will hardly ever be measured
					# Other idea: Linearly interpolate to the BAHAMAS value
					noise_strength = 0/2		# percent

					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.65*ceiling:
						print("Insufficient amount of measurements for relative humidity extrapolation at the top of the dropsonde grid (" + launch_time +
						"). There are no rel. hum. measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						# # # new_var[n,idx+1:] = np.mean(new_var[n,idx-9:idx+1]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_var[idx+1:] = (new_var[idx] + (bah_RH - new_var[idx]) / (bah_alt - alt[idx]) * (alt[idx+1:] - alt[idx]) +
							np.random.normal(0.0, noise_strength, n_alt_new - idx-1))
						# # # new_var[idx+1:] = 1.5 + np.random.normal(0.0, noise_strength, n_alt_new - idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag
						new_var[np.argwhere(new_var[:] < 0)] = 0.0



				# # # # # # More variables may be added here, if desired.


				new_dict[key] = new_var

			return new_dict, new_ipflag_dict, obs_height


		def std_extrapol(old_dict, ill_keys, old_ipflag_dict=dict()):
		# Will extrapolate some atmospheric variables to the ceiling of the dropsonde; old_ipflag will be updated.
		# Needs the old variable, the interpolation flag (should've been generated by fill_gaps()), the key and height levels as INPUT

			new_dict = old_dict
			n_alt = len(new_dict['Z'])

			new_ipflag_dict = old_ipflag_dict


			# to get the obs_height: find highest ... (?)
			if np.isnan(new_dict['flight_height']):
				# select the highest non nan index of Z. 
				highest_nonnan_Z = np.argwhere(~np.isnan(new_dict['Z']))[-1]
				obs_height = np.floor(new_dict['Z'][highest_nonnan_Z[0]]/100)*100
			else:	# use the flight_height
				obs_height = (np.floor(new_dict['flight_height']/100)*100)
			

			ceiling = obs_height	# last entry of altitude
			if ceiling > 15000:
				print("Ceiling appears to be > 15000 m. Aborted extrapolation to dropsonde ceiling because the tropopause may intervene.\n")
				return new_dict, new_ipflag_dict


			# any value above obs_height will be deleted: So if e.g. Z has got values above obs_height, delete them:
			# find the first index that overshoots obs_height:
			with warnings.catch_warnings():
				warnings.simplefilter("ignore")
				overshoot = np.argwhere(new_dict['Z'] >= obs_height)
			ill_keys.append('Z')
			if len(overshoot) > 0:
				overshoot = overshoot[0][0] + 1

				for key in ill_keys:

					if new_dict[key].ndim > 0:		# otherwise: error when using len()

						if len(new_dict[key]) == n_alt:
							new_dict[key] = new_dict[key][:overshoot]	# limit variable to obs_height
							if key in ill_keys or key == 'Z':
								new_ipflag_dict[key] = new_ipflag_dict[key][:overshoot]


			# at the end of 'Z' there may still be nans -> so we don't know to which altitude meteorological variables belong to in this region:
			# therefore: delete it and replace by extrapolation:
			n_alt = len(new_dict['Z'])
			last_nonnan_alt = np.argwhere(~np.isnan(new_dict['Z']))[-1][0]
			if ceiling - new_dict['Z'][last_nonnan_alt] > 1000:
				print("WARNING: highest GPS altitude measurement is at least 1000 m below the aircraft. Extrapolation may be erroneous.\n")

			for key in ill_keys:

				if new_dict[key].ndim > 0:		# otherwise: error when using len()

					if len(new_dict[key]) == n_alt:
						new_dict[key] = new_dict[key][:last_nonnan_alt+1]	# limit variable to obs_height
						if key in ill_keys or key == 'Z':
							new_ipflag_dict[key] = new_ipflag_dict[key][:last_nonnan_alt+1]


			# extend the old height grid up to the ceiling if the distance is greater than 10 meters:
			alt = new_dict['Z']
			n_alt = len(alt)
			alt = np.append(alt, np.arange(alt[np.argwhere(~np.isnan(alt))[-1]]+10, ceiling+1, 10))
			n_alt_new = len(alt)

			# update the altitude variable in the dictionary: & update ipflag for gpsalt:
			new_dict['Z'] = alt
			new_ipflag_dict['Z'] = np.append(new_ipflag_dict['Z'], np.ones((n_alt_new - n_alt,)))


			np.random.seed(42)		# needed for the added noise
			launch_time = datetime.datetime.utcfromtimestamp(new_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S") # for printing

			ill_keys.pop()
			for key in ill_keys: # .pop() removes 'Z' from ill_keys

				new_var = new_dict[key]
				# must be expanded to the new height grid (up to the new ceiling)
				new_var = np.append(new_var, np.full((n_alt_new - n_alt,), np.nan), axis=0) # append nans at the top of the profile

				if not new_ipflag_dict: # in case fill_gaps(...) wasn't called before this one, it's assumed that nothing has been interpolated yet.
					new_ipflag_dict[key] = np.zeros(new_var.shape)

				else: # new_ipflag also has to be extended to the new hgt grid:
					new_ipflag_dict[key] = np.append(new_ipflag_dict[key], np.zeros((n_alt_new - n_alt,)), axis=0)


				if key == 'T':
					# Temperature: If dropsondes with measurements (ipflag = 0) from that day exist, estimate their average T gradient.
					# If it clearly deviates from the standard atmospheric T gradient, then use a modified ICAO_standard atmosphere which
					# has an adapted T gradient. Otherwise:
					# Assume standard atmosphere (shifted accordingly to avoid a jump between the last known value and the extrapolation;
					# + noise (whose strength is according to instrument noise (from Hock and Franklin 1999)). ICAO standard atmosphere taken from:
					# https://www.dwd.de/DE/service/lexikon/begriffe/S/Standardatmosphaere_pdf.pdf?__blob=publicationFile&v=3
					ICAO_standard_T = 288.15 - 0.0065*alt
					noise_strength = 0/2


					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.6*ceiling:
						print("Insufficient amount of measurements for temperature extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no temperature measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var			# then just overwrite the dictionary entry with the nonedited (but extended) variable
						continue

					if alt[idx] < 0.85*ceiling: # then use standard atmosphere (ICAO):
						new_var[idx+1:] = ICAO_standard_T[idx+1:] + (new_var[idx] - ICAO_standard_T[idx])

					else:
						# Or use mean T gradient of highest 20 measurements and continue with this gradient:
						# compute mean T gradient of highest 20 measurements:
						mean_T_grad = np.mean(np.asarray([(new_var[idx-19:idx+1] - new_var[idx-20:idx]) / (alt[idx-19:idx+1] - alt[idx-20:idx])]))
						T_continued = 288.15 + mean_T_grad*alt

						new_var[idx+1:] = T_continued[idx+1:] - (T_continued[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)


					new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'P':
					# Pressure: use hydrostatic eq. with scale height H = R <T> / g0, R = 287 J kg^-1 K^-1, g0 = 9.81 m s^-2, using the vertican mean temperature <T>:
					# p(z) = p0 exp( -z / H) (Holton, p.21); + noise (Hock and Franklin 1999)
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < ceiling/3:
						print("Insufficient amount of measurements for pressure extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no pressure measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					# MAKE SURE THAT mean TEMPERATURE CAPTURES THE ACTUAL MEAN TEMPERATURE!!
					if np.count_nonzero(~np.isnan(new_dict['T'][:])) / float(n_alt_new) <= 0.75: 		# in this case you may expect that a mean temperature would yield
																									# a bad representation of the true scale height. 0.75 was chosen arbitrarily.
						T_icao_0 = 12*np.cos(4*np.pi*np.nanmean(new_dict['lat'][:])/360) + 288.15	# strongly simplified meridional surface temperature structure
						H = 287 * np.mean(288.15 - 0.0065*alt) / 9.81		# using the ICAO standard atmosphere to compute the mean temperature

						print("Warning: Because insufficient non-nan temperature values were given for launch " +
							launch_time + ", '" + str(np.mean(288.15 - 0.0065*alt)) +
							" K' was assumed to be the mean temperature for hydrostatic pressure calculation. Can possibly be avoided if the temperature is extrapolated before the pressure.\n")

					else:
						H = 287 * np.nanmean(new_dict['T'][:]) / 9.81		# scale height

						# find index of lowest non nan pressure measurement:
						l_idx = np.argwhere(~np.isnan(new_var[:]))[0]
						p_ref = new_var[l_idx]		# in Pa
						alt_ref = alt[l_idx]
						p_hydrostat = p_ref * np.exp(-(alt - alt_ref) / H)		# in Pa

						new_var[idx+1:] = p_hydrostat[idx+1:] - (p_hydrostat[idx] - new_var[idx]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'u_wind' or key == 'v_wind':
					# Wind: idea: fill nan values with the mean wind gradient of the highest 20 (non-nan)measurents. It will only be extrapolated if the the last non-nan entry
					# is higher than 0.80*ceiling:
					# other idea: just keep the same wind value
					noise_strength = 0/2

					# find highest non nan value if it lies below the ceiling:
					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.8*ceiling:
						print("Insufficient amount of measurements for wind extrapolation at the top of the dropsonde grid (" + launch_time +
							"). There are no wind measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						# # # extra_speed_length = 20		# amount of indices used for wind speed gradient calculation

						# # # for k in range(idx, n_alt_new):
							# # # new_var[n,k] = new_var[n,idx] + (k-idx)*(new_var[n,idx] - new_var[n,idx-extra_speed_length]) / (extra_speed_length + (k-idx))

						# alternative: just use the latest value for higher altitudes:
						new_var[idx+1:] = new_var[idx]
						new_var[idx+1:] = new_var[idx+1:] + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag


				elif key == 'RH':
					# Relative humidity (RH): Idea: fill nan values with the mean RH of the highest 10 measurements but only if the highest measurement exceeds or is equal to 0.90*ceiling!
					# A greater range probably doesn't make sense due to the high variability of relative humidity.
					# Other idea: Marek's suggestion: use RH approx 0 for greater altitudes. I have decided not to use rh = 0 % because 0 % relative humidity will hardly ever be measured
					# Other idea: Linearly interpolate to the BAHAMAS value
					noise_strength = 0/2		# percent

					idx = np.argwhere(~np.isnan(new_var)).flatten()[-1]

					if alt[idx] < 0.65*ceiling:
						print("Insufficient amount of measurements for relative humidity extrapolation at the top of the dropsonde grid (" + launch_time +
						"). There are no rel. hum. measurements above " + str(alt[idx]) + " m.\n")
						new_dict[key] = new_var
						continue

					else:
						# # # new_var[idx+1:] = np.mean(new_var[idx-9:idx+1]) + np.random.normal(0.0, noise_strength, n_alt_new-idx-1)
						new_var[idx+1:] = 1.5 + np.random.normal(0.0, noise_strength, n_alt_new - idx-1)
						new_ipflag_dict[key][idx+1:] = 1		# setting the interpol flag
						new_var[np.argwhere(new_var[:] < 0)] = 0.0



				# # # # # # More variables may be added here, if desired.


				new_dict[key] = new_var

			return new_dict, new_ipflag_dict, obs_height


		def regridding(new_dict, obs_height, ill_keys, resolution=10):
			'''
				Regridding variables specified in ill_keys to a uniform grid 
				(from the surface up to obs_height) with a user-defined 
				resolution (in meters, default=10).
			'''

			new_alt = np.arange(0, obs_height+1, 10)

			for key in ill_keys:
				new_dict[key] = np.interp(new_alt, new_dict['Z'], new_dict[key])

			new_dict['Z'] = new_alt

			return new_dict


		def repair_surface(old_dict, ill_keys, old_ipflag_dict=dict()):
		# Filling nan values at the surface if the gap to the surface isn't too large (e.g. measurements below 150 m must exist (roughly 10-15 seconds before splash).
			new_dict = old_dict
			alt = old_dict['Z']
			n_alt = len(alt)
			new_ipflag_dict = old_ipflag_dict
			launch_time = datetime.datetime.utcfromtimestamp(new_dict['launch_time']).strftime("%Y-%m-%d %H:%M:%S")

			lim = 200		# if there are no measurements below this altitude then the extrapolation at the surface won't be performed

			if ill_keys == ['Z']:
				threshold_list = [ill_keys, [200], ['m']]
			else:
				threshold_list = [ill_keys, [5.0, 4000.0, 50.0, 0.1, 0.1, 5.0, 5.0, 1.0],
					['K', 'hPa', '%', 'deg', 'deg', 'm/s', 'm/s', 'm/s']]		# used to check if surface value deviates siginificantly from lowest measurement

			for key in ill_keys:

				new_var = new_dict[key]

				if not new_ipflag_dict: # in case fill_gaps(...) wasn't called before this one, it's assumed that nothing has been interpolated yet.
					new_ipflag_dict[key] = np.zeros(new_var.shape)


				# find the first non-nan entry
				idx = np.argwhere(~np.isnan(new_var[:]))[0][0]

				if alt[idx] < lim:
					sfc_gap = np.arange(0,idx)

					if len(sfc_gap) == 0:
						continue
					else:

						# create mean gradient of the variable of 10 measurements above the lowest measurement, or, if grid is too coarse, take 200-400 m average:
						if alt[idx+10] > 2*lim: # default: if alt[idx+10] > 400: # take lowest measurement to 400 m mean gradient
							# find index closest to 400 m:
							idx2 = np.argmin(np.abs(alt - 2*lim))

						else: # take mean grad. of 10 measurem. above lowest measurement:
							idx2 = idx+10

						mean_grad = np.mean([new_var[j+1] - new_var[j] for j in range(idx,idx2)])	# theoretically, should never be nan because fill_gaps
																										# should've fixed the holes between the first and last measurement
						for j in sfc_gap:
							new_var[idx-j-1] = new_var[idx] - mean_grad*(j+1)


						# check if sfc value not too far off the lowest measurement:
						if key == 'RH':
							if np.any(new_var[sfc_gap] < 0):
								new_var[sfc_gap] = 0
								print("Caution, '" + key + "' surface repair resulted in negative values. Manually set the missing values at the ground to 0 for launch "
									+ launch_time + ".\n")

							elif np.any(new_var[sfc_gap] > 100):
								new_var[sfc_gap] = 100
								print("Caution, '" + key + "' surface repair resulted in >100 %. Manually set the missing values at the ground to 100 for launch "
									+ launch_time + ".\n")

						threshold = threshold_list[1][threshold_list[0].index(key)]
						si_unit = threshold_list[2][threshold_list[0].index(key)]
						if np.abs(new_var[0] - new_var[idx]) > threshold:
							print("Caution, '" + key + "' surface value deviates more than " + str(threshold) + " " + si_unit + " from the lowest measurement (launch "
								+ launch_time + ").\n")


						new_ipflag_dict[key][sfc_gap] = 1

				else:
					print("No measurements below " + str(lim) + " m. Extrapolation of '" + key + "', launch " + launch_time +
						" would eventually lead to wrong assumptions at the surface. Therefore aborted.\n")
					continue

			return new_dict, new_ipflag_dict


		def mark_outliers(sonde_dict, ill_keys): # mark outliers: outliers defined when exceeding certain thresholds

			new_dict = sonde_dict

			# thresholds are defined by change of meteorol. variable with altitude: e.g. delta p / delta z
			thresholds = [0.065, 40, 2.5, 1, 1]	# in [K/m, Pa/m, %/m, ms^-1/m, ms^-1/m]
			dz = new_dict['Z'][1:] - new_dict['Z'][:-1]	# delta z

			for key in ill_keys:
				if key == 'lat' or key == 'lon':
					continue

				met_threshold = thresholds[ill_keys.index(key)]		# threshold for key

				d_met = new_dict[key][1:] - new_dict[key][:-1]		# change of meteorological variable 'key'

				with warnings.catch_warnings():
					warnings.simplefilter("ignore")
					exceed_idx = np.argwhere(np.abs(d_met / dz) >= met_threshold)
				new_dict[key][exceed_idx] = float('nan')

			return new_dict


		def plot_met_profile(sonde_dict, ill_keys, plot_path, plot_filename_base): # plots T profile and saves it in 'plot_path'

			units = ['K', 'Pa', '%']

			for key in ill_keys:	# plot each meteorological variable that has been modified:
				# Plotting after extrapolation:
				font_size = 14
				fig = plt.figure(figsize=(6,9))
				a1 = plt.axes()

				launch_date = datetime.datetime.utcfromtimestamp(sonde_dict['launch_time']).strftime("%Y%m%d_%H%M%S")
				a1.plot(sonde_dict[key], sonde_dict['Z'], linewidth=1.2, color=(0,0,0))

				titletext = r"Dropsonde " + key + " profile from EUREC4A campaign: " + launch_date
				plt.title(titletext, fontsize=font_size, wrap=True)
				a1.set_xlabel(key + " [" + units[ill_keys.index(key)] + "]", fontsize=font_size)
				a1.set_ylabel(r"Height [m]", fontsize=font_size)
				a1.grid(True, axis='x', which='both')
				a1.grid(True, axis='y', which='major')
				a1.set_ylim(bottom=0, top=sonde_dict['Z'][-1])

				if key == 'tdry':
					a1.set_xlim(left=240, right=305)
				elif key == 'pres':
					a1.set_xlim(left=10000, right=105000)
				elif key == 'rh':
					a1.set_xlim(left=0, right=100)

				plt.savefig(plot_path + plot_filename_base + "_" + key + ".png") #, dpi=250, bbox_inches='tight'
				plt.close()


		def saveExpolSondeAsNC(sonde_dict, out_filename): # saves the sonde_dict as an nc file named out_filename into data_out_path

			new_nc = nc.Dataset(out_filename, "w", format="NETCDF4")

			# create dimensions:
			n_alt = len(sonde_dict['Z'])
			new_nc.createDimension("alt", n_alt)

			# convert


			# Global attributes:
			new_nc.description = """EUREC4A campaign JOANNE Level 3 dropsondes. Extrapolated when enough measurements were given. More information required? Contact author (listed below).
				'alt_old' represents the dimension of height levels of the uninterpolated / extrapolated variables. alt represents the height levels of the extrapolated variables."""
			new_nc.history = "Created: " + datetime.datetime.utcnow().strftime("%Y-%m-%d %H:%M:%S")
			new_nc.author = "Andreas Walbroel (Mail: a.walbroel@uni-koeln.de)"


			# create variables, set attributes and save values into the created variables:
			# attributes = [[units], [description]]
			attributes = [
							['m', 'degree North', 'degree East', 'degree North', 'degree East', 'seconds since 1970-01-01 00:00:00 UTC',
								'Pa', 'K', 'm', 'm s^-1', 'm s^-1', '%'],
							['Flight altitude on dropsonde launch', 'Latitude on dropsonde launch', 'Longitude on dropsonde_launch',
								'Latitude', 'Longitude', 'Sonde launch time', 'Air pressure', 'Air temperature', 'Height',
								'Zonal wind speed (eastwards is > 0)', 'Meridional wind speed (northwards is > 0)', 'Relative humidity']
						]

			sonde_keys = ['reference_alt', 'reference_lat', 'reference_lon', 'lat', 'lon', 'launch_time', 'P', 'T', 
							'Z', 'u_wind', 'v_wind', 'RH']

			sonde_dict['reference_lon'] = sonde_dict['flight_lon']
			sonde_dict['reference_lat'] = sonde_dict['flight_lat']
			sonde_dict['reference_alt'] = sonde_dict['flight_height']
			sonde_dict['reference_time'] = sonde_dict['launch_time']

			for at_idx, key in enumerate(sonde_keys):

				# find out about this variable's dimension:
				if key in ['launch_time', 'reference_alt', 'reference_lat', 'reference_lon']:
					var = np.asarray([sonde_dict[key]])

				else:
					var = sonde_dict[key]


				if var.ndim == 1:
					if len(var) == n_alt:
						nc_var = new_nc.createVariable(key, "f8", ("alt"))

					elif len(var) == 1:
						nc_var = new_nc.createVariable(key, "f8")

					else:
						print("""Hmmm.... something went wrong. Actually all variables should either have the dimension of the
							altitude or length = 1. But """ + key + " has got another shape. Happy debugging!")

					nc_var[:] = var


				nc_var.units = attributes[0][at_idx]
				nc_var.description  = attributes[1][at_idx]

			new_nc.close()
			print("Output: ", out_filename)



		############################################################################################



		# Dropsonde quality control: if there are small gaps of measurements between launch altitude and surface, they will be filled via interpolation.
		# Additionally, the sondes will be extrapolated to a certain altitude.

		# Reading in dropsonde files:
		HALO_sondes_NC = sorted(glob.glob(path_raw_sondes + "*.nc"))
		if path_BAH_data:
			BAH_files_NC = sorted(glob.glob(path_BAH_data + "bahamas*.nc"))
		else:
			BAH_files_NC = []


		tot_failure_warnings = 0			# counts the amount of times a critical variable has got no measurements
		tot_sonde_stuck = 0
		failed_sondes = []
		stuck_sondes = []

		print('Found %d file.' % len(HALO_sondes_NC))

		for sonde_nc in HALO_sondes_NC:

			sonde_dict = readNCrawJOANNE3(sonde_nc)
			print("Input: ", sonde_nc)
			
			# Start iteration over the sondes in the sonde_nc file:
			n_sondes_in_file = len(sonde_dict['launch_time'])

			for sn in range(n_sondes_in_file):
				if 'HALO' in sonde_dict['platform'][sn]:

					# Only pass one dropsonde profile ... not all of them at once:
					sonde_dict_o = dict()
					for key in sonde_dict.keys():
						if key == 'fillValues':
							continue

						if sonde_dict[key].shape[0] == n_sondes_in_file and sonde_dict[key].ndim == 2:
							sonde_dict_o[key] = sonde_dict[key][sn,:]
						elif sonde_dict[key].shape[0] == n_sondes_in_file and sonde_dict[key].ndim == 1:
							sonde_dict_o[key] = sonde_dict[key][sn]
						elif sonde_dict[key].shape[0] == len(sonde_dict['Z']) and sonde_dict[key].ndim == 1:
							sonde_dict_o[key] = sonde_dict[key]
						else:
							raise IndexError("This index shape was rather unexpected. Something's wrong here.")

					launch_date = datetime.datetime.utcfromtimestamp(sonde_dict_o['launch_time']).strftime("%Y-%m-%d %H:%M:%S")	# time delta required
					dropsonde_date = (datetime.datetime.strptime(launch_date, "%Y-%m-%d %H:%M:%S")).strftime("%Y%m%d")	# date displayed in the filename ... comfy way to
																														# find the right BAHAMAS data for std_extrapol
					launch_date_for_filename = (datetime.datetime.strptime(launch_date, "%Y-%m-%d %H:%M:%S")).strftime("%Y%m%d_%H%M%ST")

					if launch_date_for_filename == '20200131_175736T':
						print("Skipp sonde on ", launch_date_for_filename, " as the sonde has a dry bias.")
						continue

					print("########## Day: " + launch_date + " ##########\n")

					# add another condition that checks if e.g. nearly no measurements exist at all (for T, P and RH):
					if np.any([np.count_nonzero(~np.isnan(sonde_dict_o['T'])) < 0.1*len(sonde_dict_o['T']),
						np.count_nonzero(~np.isnan(sonde_dict_o['P'])) < 0.1*len(sonde_dict_o['P']),
						np.count_nonzero(~np.isnan(sonde_dict_o['RH'])) < 0.1*len(sonde_dict_o['RH']),
						np.count_nonzero(~np.isnan(sonde_dict_o['lat'])) < 0.05*len(sonde_dict_o['lat']),
						np.count_nonzero(~np.isnan(sonde_dict_o['lon'])) < 0.05*len(sonde_dict_o['lon']),
						np.count_nonzero(~np.isnan(sonde_dict_o['u_wind'])) < 0.05*len(sonde_dict_o['u_wind']),
						np.count_nonzero(~np.isnan(sonde_dict_o['v_wind'])) < 0.05*len(sonde_dict_o['v_wind'])]):
						tot_failure_warnings = tot_failure_warnings + 1
						failed_sondes.append(sn)
						print("One PAMTRA-critical variable measurement failed. Skipping this dropsonde.\n")
						continue

					# add jet another condition that checks if the sonde got stuck mid air (gps alt. values don't really decrease with time):
					if not np.any(sonde_dict_o['Z'][~np.isnan(sonde_dict['Z'])] < 1500):	# then I assume that the whole launch was doomed
						print("Sonde got stuck in mid air. 'Z' doesn't seem to include any values < 1500 m.\n")
						stuck_sondes.append(sn)
						tot_sonde_stuck = tot_sonde_stuck + 1
						continue


					# it's known that the altitude axis (gpsalt) does have some broken values (e.g. sudden jumps or exceeding the aircraft altitude):
					# but before we deal with those, the 'normal' gaps in gpsalt shall be fixed


					# subsequent variables will be cured from holey nan value disease...:
					# pressure, temperature, relhum, wind (u & v & w), lat, lon.
					# Other variables, of which you can expect a linear interpolated over gaps to be applicable, may be added.
					ill_keys = ['T', 'P', 'RH', 'lat', 'lon', 'u_wind', 'v_wind']


					sonde_ipflag = dict()		# will contain the interpolation flags for interpolated nan values in the middle of the drop

					# before filliing the gaps, the altitude axis must be fixed:
					sonde_dict_o['Z'], sonde_ipflag['Z'] = fill_gaps(sonde_dict_o['Z'])
					for key in ill_keys:
						sonde_ipflag[key] = np.full_like(sonde_dict_o[key], 0)
						sonde_dict_o[key], sonde_ipflag[key] = fill_gaps(sonde_dict_o[key])		# altitude must be passed to check for dimensions of the to-be-cured variable...

					sonde_dict_o['ipflag'] = sonde_ipflag



					# now we still need to handle the nan values at the surface: perform surface repair for the atmospheric 
					# parameters If there are no non-nan values in the lowest 5 % of the variable --> don't interpolate 
					# because the assumption would eventually lead to senseless surface values:
					sonde_dict_o, sonde_dict_o['ipflag'] = repair_surface(sonde_dict_o, ill_keys, sonde_dict_o['ipflag'])



					# Extrapolating the ill_keys to the ceiling of the dropsondes (e.g. below aircraft altitude):
					# CAUTION: it is expected that the dropsondes start BELOW THE TROPOPAUSE!
					# It is advisable to do the TEMPERATURE EXTRAPOLATION FIRST because it improves the pressure extrapolation.
					# Need bahamas file for extrapolation limit:
					if BAH_files_NC:
						bah_filename = [bah_file for idx, bah_file in enumerate(BAH_files_NC) if dropsonde_date in bah_file]
						sonde_dict_o, sonde_dict_o['ipflag'], obs_height = std_extrapol_BAH(sonde_dict_o, ill_keys, bah_filename, sonde_dict_o['ipflag'])

					else:
						sonde_dict_o, sonde_dict_o['ipflag'], obs_height = std_extrapol(sonde_dict_o, ill_keys, sonde_dict_o['ipflag'])

					# Regridding to a uniform vertical grid with a user-specified resolution:
					sonde_dict_o = regridding(sonde_dict_o, obs_height, ill_keys, 10)


					# find outliers and mark them (as nan): afterwards fill them again
					sonde_dict_o = mark_outliers(sonde_dict_o, ['T', 'P', 'RH', 'u_wind', 'v_wind'])
					for key in ill_keys:
						sonde_dict_o[key], sonde_ipflag[key] = fill_gaps(sonde_dict_o[key])


					# Save the extrapolated sonde dictionary to a new nc file:
					data_out_path_halo_dir = os.path.dirname(data_out_path_halo)
					if not os.path.exists(data_out_path_halo_dir):
						os.makedirs(data_out_path_halo_dir)
					out_filename = os.path.basename(sonde_nc)		# removes the path in the string HALO_sondes_NC[m] so that the filename remains
					out_filename = os.path.join(data_out_path_halo, out_filename[0:-3] + "_" + launch_date_for_filename + "_v01.nc")

					saveExpolSondeAsNC(sonde_dict_o, out_filename)


				# # # # # # save the lists containing the broken and stuck files:
				# # # # # np.save(data_out_path_halo + "tot_failure_warnings.npy", np.asarray(tot_failure_warnings))
				# # # # # np.save(data_out_path_halo + "tot_sonde_stuck.npy", np.asarray(tot_sonde_stuck))



	if dropsonde_dataset == "raw":
		run_dropsonde_gap_filler_raw(path_raw_sondes, data_out_path_halo, path_BAH_data)
	elif dropsonde_dataset == 'joanne_level_3':
		run_dropsonde_gap_filler_joanne(path_raw_sondes, data_out_path_halo, path_BAH_data)
	else:
		raise ValueError("""The variable 'dropsonde_dataset' must be a string. 
			Only the following options are permitted: 'raw', 'joanne_level_3'""")