Error using netcdflib NC_EINVALCOORD

Hi,

I'm trying to get some wind direction timelines for different heights. The dataset includes wind data for 37 pressure levels. And I'm getting the following error:

h =

     1

Error using netcdflib
The NetCDF library encountered an error during execution of 'getVarsShort' function - 'Index exceeds dimension bound (NC_EINVALCOORDS)'.

Error in netcdf.getVar (line 140)
    data = netcdflib(funcstr,ncid,varid,varargin{:});

Error in internal.matlab.imagesci.nc/read (line 635)
                data  = netcdf.getVar(gid, varid, ...

Error in ncread (line 66)
vardata = ncObj.read(varName, varargin{:});

Error in DirectionTimelines_from_ERAInt6hrly_SJ (line 44)
        uwnd_all_temp = ncread(Wu(f).name,'u', [(int16((VolcLong/0.75)+1)) (int16(241-((90+VolcLat)/0.75))) h 1],[1 1 1 (length(ncread(Wu(f).name,
        'time')))],[1 1 1 1]);      %Some years have a leap year and thus 1464 time periods, but leave at 1460 so doesn't need adapting for leap years.

Error in run (line 91)
evalin('caller', strcat(script, ';'));

The code is the following:

cd ~/Documents/MATLAB/Wind
Wu = dir('u*.nc');
Wv = dir('v*.nc');
YearRange = '2009-2018';
NumYears = 10;  
Levels = length(ncread(Wu(1).name, 'level'));
for f = 1:length(Wu) 
    if (f==1)
        date(1:16072) = ncread(Wu(f).name, 'time');    
    else
        date((length(date)+1):((length(date)+1)+length(ncread(Wu(f).name, 'time'))-1)) = ncread(Wu(f).name, 'time');
    end
end
Time = length(date);

%Change for each volcano run...
Volcano = '5S'; VolcLat = -5; VolcLong = 95;
mkdir(strcat('../../Output-Wind/',Volcano,'_WindInfo_2019-2018'));

if (VolcLong < 0)   %i.e. in Western Hem
   VolcLong = 360 + VolcLong; 
end

uMatrix = zeros(Time, 37);  %Matrix of n rows (1460 time records per year for 6 hourly daily data) by 37 (heights)

%To convert mb or hPa (unit for pressure level data) to altitude, try: www.csgnetwork.com/pressurealtcalc.html
for h = 1:37        %Number of levels
    h               %So can track progress when running
    count = 0;
    for f = 1:length(Wu)
        uwnd_all_temp = ncread(Wu(f).name,'u', [(int16((VolcLong/0.75)+1)) (int16(241-((90+VolcLat)/0.75))) h 1],[1 1 1 (length(ncread(Wu(f).name, 'time')))],[1 1 1 1]);      %Some years have a leap year and thus 1464 time periods, but leave at 1460 so doesn't need adapting for leap years.
        uwnd_all_temp = squeeze(uwnd_all_temp);
        uMatrix((count+1):(count+length(uwnd_all_temp)),h) = uwnd_all_temp;  
        count = count+length(uwnd_all_temp);
    end
end

vMatrix = zeros(Time, 37);  %Matrix of rows (no. years * no. days * no. records/day) by 37 cols (heights)

for h = 1:37        %Number of levels
    h
    count = 0;
    for f = 1:length(Wv)
        vwnd_all_temp = ncread(Wv(f).name,'v', [(int16((VolcLong/0.75)+1)) (int16(241-((90+VolcLat)/0.75))) h 1],[1 1 1 (length(ncread(Wu(f).name, 'time')))],[1 1 1 1]);      %Some years have a leap year and thus 1464 time periods, but leave at 1460 so doesn't need adapting for leap years.
        vwnd_all_temp = squeeze(vwnd_all_temp);
        vMatrix((count+1):(count+length(vwnd_all_temp)),h) = vwnd_all_temp;  
        count = count+length(vwnd_all_temp);
    end
end

Dmatrix = zeros(Time,Levels);       % matrix of direction (rows) at each height (cols)
Imatrix = zeros(Time,Levels);       % matrix of intensity (rows) at each height (cols)

AltitudeLabels = {'32.5 km' '31 km' '29.5 km' '28 km' '27 km' '26 km' '23.5 km' '21.5 km' '19.5 km' '17.5 km' '16 km' '14.5 km' '13.5 km' '12.5 km' '12 km' '11 km' '10.5 km' '9 km' '8 km' '7 km' '6.5 km' '5.5 km' '5 km' '4.2 km' '3.6 km' '3 km' '2.5 km' '2.2 km' '1.9 km' '1.7 km' '1.5 km' '1.2 km' '1 km' '800 m' '500 m' '300 m' '100 m'};

%To calculate and save direction and intensity per record with height:
for h = 1:Levels  
    Dmatrix(:, h) = mod(90-atan2d(vMatrix(:,h),uMatrix(:,h)), 360); %calculation from u and vmatrix... for Direction.
    Imatrix(:, h) = sqrt((uMatrix(:,h).^2)+(vMatrix(:,h).^2)); %calculation from u and vmatrix... for Speed.
end
save(strcat('../../Output-Wind-Info/',Volcano,'_WindInfo_2009-2018/Dir-Int_Matrix.mat'),'Dmatrix','Imatrix');

%From: http://www.csgnetwork.com/pressurealtcalc.html
Altitude = {32435,30760,29675,28180,27115,25905,23315,21630,19315,17660,15790,14555,13505,12585,11770,11030,10360,9160,8115,7180,6340,5570,4865,4205,3590,3010,2465,2205,1950,1700,1455,1220,990,760,540,325,110};
AltAsNumber = cell2mat(Altitude);
AltAsNumber = reshape(AltAsNumber,37,1);

Has anyone an idea how to solve the problem?

Thanks in advance

Steffen