function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history

	dark = cal_pack.uv_channel.dark

	dit = header.dit
	nbin = header.nbin 
	ndit1 = header.ndit1
	ndit2 = header.ndit2
	obj_cnt = header.obj_cnt
	obt_beg = header.obt_beg
	masking = header.masking
	tsensor = header.tsensor
	
	; WARN - temporary patch to handle local l1 fits files
	if tsensor eq 0. then tsensor = -25.

	flag_use_extrapolated = 1
	flag_use_notfullset = 1
	flag_normalize_dark = 0

	obt_available = !null

	; check if one or more dark matrices with same dit, ndit1, ndit2, and tsensor exist
	for i = 0, n_elements(dark) - 1 do begin
		if $
			(dark[i].dit eq dit) and $
			(abs(dark[i].tsensor - tsensor) lt 5) and $
			(dark[i].ndit1 eq ndit1) and (dark[i].ndit2 eq ndit2) and $
			flag_use_notfullset and flag_use_extrapolated  then $
			obt_available = [obt_available, dark[i].obt_beg]
	endfor

	; if not, exit with error code
	if ~ isa(obt_available) then begin
		journal, 'Error 12: UV applicable dark file not found.'
		journal
		exit, status = 12
	endif

	; if yes, chose the temporally closest
	delta_obt = min(abs(obt_beg - obt_available), j)
	i = where(dark.obt_beg eq obt_available[j])

	transient = where(dark[i].filename.cnt_1 ne '', n)
	
	; select the dark matrix corresponding to the correct transient phase
	if obj_cnt le n then $
		dark_file = dark[i].filename.cnt_1[obj_cnt - 1] else $
		dark_file = dark[i].filename.cnt_2

	; apply normalization if required
	if flag_normalize_dark and masking.contains('disabled', /fold) then $
		dark_image = correzione_uvda(data, cal_pack.path + dark_file) else $
		dark_image = float(readfits(cal_pack.path + dark_file, /silent))

	dark_nobin = dark_image

	; rebin the dark and take into account the exposure time correction of l1 images
	dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin
	dark_image = dark_image * ndit1 * ndit2

	mask = where(data eq 0., count)
	data = data - dark_image
	
	; set values in masked areas equal to zero
	if count gt 0 then data[mask] = 0.

	; read the dark error matrix
	dark_error_file = dark[i].error
	if file_test(dark_error_file) then $
		dark_error = float(readfits(cal_pack.path + dark_error_file, /silent)) else $
		dark_error = fltarr(1024, 1024)

	; rebin it and take into account the exposure time correction
	dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin)
	
	; calculate the total uncertainty
	if isa(error) then error += (data/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2

	jump:
	s = where(data le 0., count)
	quality_matrix[s] = 0
	
	if ~ isa(history) then history = !null
	history = [history, 'Dark-current correction: ', '  ' + dark_file]

	journal, 'Dark-current correction:'
	journal, '  dit = ' + string(dit, format = '(I0)') + ' ms'
	journal, '  ndit1 = ' + string(ndit1, format = '(I0)')
	journal, '  ndit2 = ' + string(ndit2, format = '(I0)')
	journal, '  obj_cnt = ' + string(obj_cnt, format = '(I0)')
	journal, '  tsensor = ' + string(tsensor, format = '(F0)') + ' degC'
	journal, '  dark file = ' + dark_file
	journal, '  extrapolated dark = ' + dark[i].extrapol.tolower()

	return, data
end