function [upsi,vpsi,levels,lats,lons,masterpsi,mastervort,mastertime]=non_div_wind(syear,eyear)

%This script returns the non-divergent wind using the ERA-interim relative
%vorticity and lat/lon grids as an input

close all
%clear all

%MAP STUFF
load conus;
S=shaperead('landareas','UseGeoCoords',true);

%add a bunch of paths to make sure netcdf handling is easier
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/'); %we call some scripts from here
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/mexnc/') % netcdf
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/snctools'); %netcdf
addpath('/uufs/chpc.utah.edu/common/home/horel_data/nlareau/ERA_northern_hemi'); %data directory
datadir='/uufs/chpc.utah.edu/common/home/horel_data/nlareau/ERA_northern_hemi/'; %data directory as a sting

%% start by loading lat/long 1.5 deg grid.
tfile='terrain.nc';
terrain=nc_varget(tfile,'z')./9.8; %divide by g 'cause it is geopotential
tlats=double(nc_varget(tfile,'latitude'));
tlons=double(nc_varget(tfile,'longitude'));
lons=tlons;
px=length(lons); %number of points in the zonal direction

lats=tlats;
py=length(lats); %number of points meridionally

[~,latmat]=meshgrid(lons,lats);

%define the earth radius and grid spacings
a=6.37*10^6;%meters
dy1=(1.5).*(pi./180).*a;% constant north-south grid spacing (y positive toward the pole)
dx=(1.5).*a.*cosd(latmat(2:end-1,:)).*(pi./180); %grid of zonal spacing (between lines of longitude);
dy=repmat(dy1,size(dx));
twodx=dx; %used to multiply by two... now this is just redudant naming
twody=dy;
dxsq=twodx.^2; %sqaured spacing for second derivatives
dysq=twody.^2;

% now we need to decide how many levels we'll actually be using:
lvls=9; % define # levs here
%now decide on the pressure level spacing (probably an even (100 hpa)
dp=100*10^2; % Pa
twodp=dp.*2;
twodp=ones(lvls-2,1)*twodp;
dpsqr=twodp.^2;
dxsq3d=repmat(dxsq,[1 1 lvls-2]); % tile the squared spacing into the height dimension
dysq3d=repmat(dysq,[1 1 lvls-2]);
dpsq3d=repmat(dpsqr,[1 py-2 px]);
dpsq3d=permute(dpsq3d,[2 3 1]);
dxflat=dxsq3d(:);
dyflat=dysq3d(2);
dpflat=dpsq3d(:);
clear dpsq3d;

% coriolis parameter
f0=10^-4;
beta=f0./a;
bplane=f0+beta.*(a*lats);
fsqr=(f0).^2;
rot=7.292*10^-5; %angular velocity of earth
latter=repmat(latmat(:,:),[1 1 lvls-2]);
clear latmat;
coriolis=(2*rot).*sind(latter);
cflat=coriolis(2:end-1,:,:).^2;
cflat=cflat(:);
clear coriolis latter;

% for constant sigma
sg=2.5*10^-6; % mean tropospheric value from Holton;

%% now construct the basic 3-d Laplacian differential operator which we will use lots later in the code
nx=py;
ny=px;
nz=lvls;

[~,~,L3]=laplacian([nx-2 ny nz-2],{'DD','P','DD'}); %shaped as (lats,lons,heights) %periodic in the zonal direction

%now pull out the diagnols and modify them to reflect irregular grids
[B2,d2]=spdiags(L3);
B2(:,4)=B2(:,4).*(1./dyflat); %use the dy spacing here
B2(:,6)=B2(:,6).*(1./dyflat); %use the dy spacing here
B2(:,1)=B2(:,1).*(0); % no vertical dervatives
B2(:,9)=B2(:,9).*(0); % no vertical dervatives
B2(:,3)=B2(:,3).*(1./dxflat); % dx spacing
B2(:,7)=B2(:,7).*(1./dxflat); % dx spacing
B2(:,8)=B2(:,8).*(1./dxflat);
B2(:,2)=B2(:,2).*(1./dxflat);
B2(:,5)=((2)./dxflat)+((2)./dyflat);
B2=sparse(B2); %make sure it is still sparse
L2d=spdiags(B2,d2,L3); %repopulate the laplacian

%% loop over years (or specify a year of interest)
for yr=syear:eyear
    mastertime=[]; %initialize these variables
    masterpsi=[];
    mastervort=[];
    if yr<2011
        halfnum=2; %use two 6 month netcdf data files
    else
        halfnum=1; % use just 1 6 month file 'cause the data ain't out yet for the second half of the year
    end
    for half=1:halfnum
        if half ==1
            monther='jan_jun_';
        else
            monther='jul_dec_';
        end
        
        
        
        
        file=strcat('output_',monther,num2str(yr),'.nc');
        
        %% Pull in the Data from the NETCDF files...
        
        time=double(nc_varget(file,'time'))./24 + datenum(1900,01,01,00,00,00); %convert to matlab serial time
        mastertime=[mastertime;time]; %concatenate the time into a master time file to be saved with the output data
        
        %mnlst=unique(mn); %list the unique month numbers
        %for mni=1:length(mnlst)
        %    mnidx=find(mn==mnlst(mni));
        
        rvort=double(nc_varget(file,'vo'));%,[mnidx(1)-1 0 0 0],[length(mnidx) Inf Inf Inf]));
        leveler1=size(rvort,2);
        hgt=double(nc_varget(file,'z'));%,[mnidx(1)-1 0 0 0],[length(mnidx) Inf Inf Inf]))./9.8;
        
        %pressure levels
        levels=double(nc_varget(file,'levelist'))*10^2;
        [tsteps,leveler,ltdim,lndim]=size(rvort);
        
        levmat=repmat(levels,[1 tsteps ltdim lndim]);
        levmat=permute(levmat,[2 1 3 4]);
        
        %% Now we need to get the stream function and non-divergent winds
        tic
        clear psi psii psi2
        
        rvort=permute(rvort,[1 3 4 2]); % move the levels dim to the end to be the same as the laplacian above
        psi2=zeros(size(rvort(:,2:end-1,:,2:end-1)));
        ruse=rvort(:,2:end-1,:,2:end-1);
        matlabpool open 8
        parfor ii=1:tsteps
            rnow=squeeze(ruse(ii,:,:,:)); %grab one time step and drop the time dim.
            rnow=rnow(:);%flatten it
            U=L2d\rnow; %this is where the magic happens... MATLAB LEFT DIVIDE INVERTS THE MATRIX
            psii=reshape(U,nx-2,ny,nz-2);
            psi2(ii,:,:,:)=psii;
        end
        matlabpool close
        
        toc %should take about 60 seconds
        
        psi=zeros(size(rvort))*nan;
        psi(:,2:end-1,:,2:end-1)=psi2;
        clear psi2 psii ruse;
        psi=permute(psi,[1 4 2 3]);
        masterpsi=cat(1,masterpsi,psi);
        mastervort=cat(1,mastervort,rvort);
        
    end
end

%%

[dpsidx,dpsidy]=met_gradientm_peri(masterpsi,lats,lons);
vpsi=-dpsidx;
upsi=dpsidy;


