function [out]=zint4(in,plev);

%calculates the vertical density average via pressure integration at each grid point and time
%input is a 4d matrix (level,time,lat,lon)
%plev is the levels corresponding to the first dimension of input matrix

[nlev,ntime, NN, N] = size(in);

if ((length(plev) == nlev) == 0)
    error('vector of pressure values must have same length as 1st matrix dimension')
end

sgn=(plev(2)-plev(1))/(abs(plev(2)-plev(1))); %are the levels in ascending or descinding order

for n=1:nlev-1
    if  ((plev(n+1)-plev(n)) > 0) == (.5 - .5*sgn) %do the levels jump up and down
        error('levels must be monotonic')
    end
end


if sgn == -1
    plevo=plev;
    plev=sort(plev);
    display('levels are in descending order')
else
    plevo=plev;
    display('levels are in ascending order')
end
    
    
display(['data matrix has ', num2str(nlev) , ' levels ', num2str(ntime),' time iterations and a horizontal mesh of ', num2str(NN),' latitudes and ', num2str(N), ' longitudes'])
display('most of which does not matter anyways as long as the number of levels is right')

out=zeros(ntime,NN,N);
if sgn == 1
    for nt=1:ntime
        for nn=1:NN
            for n=1:N
                tot=0;
                num=max((find(isnan(in(:,nt,nn,n)) == 0))); %number of data points for integration
                mx = max(plev(find(isnan(in(:,nt,nn,n)) == 0)));
                mn=min(plev);
                out(nt,nn,n) = trapz(plev(1:num), in(1:num,nt,nn,n))/(mx-mn);
            end
        end
    end


else
    for nt=1:ntime
        for nn=1:NN
            for n=1:N
                tot=0;
                num=min((find(isnan(in(:,nt,nn,n)) == 0))); %number of data points for integration
                mx = max(plev(find(isnan(in(:,nt,nn,n)) == 0)));
                mn=min(plev);
                out(nt,nn,n) = trapz(plevo(nlev:-1:num),in(nlev:-1:num,nt,nn,n) )/(mx-mn);
            end
        end
    end
end