DFT for signals of different length and sampling rate
Fft normalization is often an overlooked topic and many articles dealing with FFT just say to divide the fft ouput by N (divide by either N or N/2 depending on the actual
FFT algorithm used). 
Real FFT  usually discards the upper part of the spectrum (due to FFT redundancy) so the normalization coefficient here is N/2, rather than N.


  So, my question is: what normalization coefficient for PSD should I
  use in my case?


I prefer the approach taken by Heinzel et al, and according to them, the normalized/scaled PSD  can be computed as

PSD = ( 2 * (yi * yi) ) / ( fs * S2 ); // i = 0…N/2  (1)

S2 = sum(wi * wi); //i = 0...N (2)

(yi is the complex result of the fft at bin i and wi is the value of the window function at point i )


The sum S2 is used instead of N in order to better estimate the scaling
factor, given a particular window (rectangular, Hann etc). Using no windows at all corresponds to a rectangular window and in the case of a rectangular window S2 is equal to N.

Heinzel et al also define a metric called ENBW (equivalent noise bandwidth), which is at the heart of the here defined scaling technique and for those interested, complete derivations (including 1 and 2) can be found here https://holometer.fnal.gov/GH_FFT.pdf

Note also that for a rectangular window, ENBW is equal to fsampling / N, which may explain why the results started to make sense for you once this normalization factor was applied.