% function [xx,yy, field_out,spec] = make_psi(lon, lat,field, trunc);
clear all;close all;

Y=1980; MN=4; g=9.81;

cd /home/disk/viento/aaron/research/ncep_heat_trans_trentberth/data
fil='T63t_MRES_1979-2005_MM.nc'; 
nc=netcdf(fil);

lon=nc{'lon'}(:); %
lat=nc{'lat'}(:);
time=nc{'time'}(:);
year=rounddown(time/100); %kosher conversion
month=time-100*year; %ditto

dot=find(year == Y & month == MN)

MRES=squeeze(nc{'MRES'}(dot,:,:))/g;
sphcolor2(MRES, lat,lon,19); title('Trenberth Mass residual')


trunc=42;
field=MRES/(10^5/g); %mass averaged divergence

%routine to make a velocity potential given a divergence 


[NN,N]=size(field);


%if nargin<4; error('Error - Not enough input parameters specified'), end;
%if nargin>4; error('Error - to0 many input parameters specified'), end;
if length(lat)>NN;error('Error - vector of latitude values has too large of dimensions'),end;
if length(lat)<NN;error('Error - vector of latitude values has too small of dimensions'),end;
if length(lon)>N;error('Error - vector of longitude values has too large of dimensions'),end;
if length(lon)<NN;error('Error - vector of longitude values has too small of dimensions'),end;
if max(lat)>pi; lat=lat*pi/180; display('converting input latitudes to radians'); end;
if max(lon)>2*pi; lon=lon*pi/180; display('converting longitude input to radians'); end; 


[LON LAT]=meshgrid(lon,lat);
GPy=round((3*trunc+1)/2+.999);   %number of Gaussian latitudes
[yy, wy]=gaussquad(GPy);          %sin(lat) and weights of Gaussian quadrature points 

%now do the same in x
GPx=2*GPy;   
xx=0:360/GPx:360-360/GPx;
xx=xx*pi/180;
xx=xx';

[XX YY]=meshgrid(xx,asin(yy));
[XX WW]=meshgrid(xx,wy);

%now interpolate the field onto the gaussian grid
ifield=interp2(LON, LAT, field, XX, YY);


N=length(xx); %new dimensions of the guassian quadrature field
NN=length(yy);

a=6.37*10^6; %radius of the Earth

num=trunc+1:-1:1; 
NNN=sum(num); %total number of basis functions
nor=zeros(NNN,1); %pre-allocate for the normalization
yy=asin(yy); %redefine the y coordinate
le=zeros(NN, trunc+1, trunc+1); %pre-allocate for legendre polynomials 
dle=zeros(NN, trunc+1, trunc+1); % pre - allocate for the anayltically defined y derivative of legendre polynomials

%put together the legendre polynomials and derivatives  as a function of n
%and m
for n=0:1:(trunc)
    P=legendre(n,sin(yy));
    for m=0:1:n;
        le(:,n+1-m, m+1)=(P(m+1,:))';
        nnn(n+1-m, m+1)=n;
        dle(:,n+1-m,m+1)=divleg(m,n+1-m,yy);
    end
end


%now put them in a column list- just pre sort some stuff here
M=0;
xx=xx';
LE=zeros(NN, NNN);
dLE=zeros(NN, NNN);
for m=0:trunc
    for n=1:1:(1+trunc-abs(m))
        nor(M+n,1)=sum(sum(((le(:,n,abs(m)+1)*cos(m*xx)).^2).*WW));
        nn(M+n,1)=nnn(n,abs(m)+1);
        mm(M+n,1)=m;
        LE(:,M+n)=le(:,n, abs(m)+1)/mean(abs(le(:,n, abs(m)+1)));%%%%%%%%%%%%%%%%%%%%% added 7/28/08 AD- to prevent huge numbers in basis
        dLE(:,M+n)=dle(:,n,abs(m)+1)/mean(abs(le(:,n, abs(m)+1))); %also have to add this into the derivative of the legendre
    end
    M=M+n;
end

NNN=sum(num)+sum(num(2:length(num))); %would be 2*sum(num) but some are unneccesary


norm=zeros(NNN,1);
nnn=zeros(NNN,1);
offset=sum(num(2:length(num)));

for n=1:sum(num)
    basis(:,:,n)=LE(:,n)*cos(mm(n)*xx);
    norm(n,1)=sum(sum((basis(:,:,n).^2).*WW)); 
    nnn(n,1)=-nn(n)*(nn(n)+1)/(a^2);
    tn(n,1)=nn(n); %totalwavenumber
    vvpbasis(:,:,n)=dLE(:,n)*cos(mm(n)*xx); %defines the meridional velocity associated with a velocity potential spectral coefficient 
    for m=1:1:NN
        uvpbasis(m,:,n)=(-mm(n)/(a*cos(yy(m))))*LE(m,n)*sin(mm(n)*xx); %defines the meridional velocity associated with a velocity potential spectral coefficient -code copied from prefilter (simple loop over y ? AD 4/2009)? 
    end
    
    basis(:,:,n+offset)=LE(:,n)*sin(mm(n)*xx);
    norm(n+offset,1)=sum(sum((basis(:,:,n).^2).*WW));
    nnn(n+offset,1)=-nn(n)*(nn(n)+1)/(a^2);
    tn(n+offset,1)=nn(n); %total wavenumber
    vvpbasis(:,:,n+offset)=dLE(:,n)*sin(mm(n)*xx); %defines the meridional velocity associated with a velocity potential spectral coefficient 
    for m=1:1:NN
        uvpbasis(m,:,n+offset)=(mm(n)/(a*cos(yy(m))))*LE(m,n)*cos(mm(n)*xx); %defines the meridional velocity associated with a velocity potential spectral coefficient -code copied from prefilter (simple loop over y ? AD 4/2009)?
    end
    
    
end

for n=2:1:NNN
    %spec(n,1)=sum(sum(basis(:,:,n).*WW.*ifield))/(norm(n)); %just divide by the total wavenumber instead of multiplying  
    spec(n,1)=(1/nnn(n))*sum(sum(basis(:,:,n).*WW.*ifield))/(norm(n)); %just divide by the total wavenumber instead of multiplying    
end

field_out=zeros(NN,N);
V_out=zeros(NN,N);
U_out=zeros(NN,N);


for n=2:1:NNN
    field_out=field_out+spec(n,1)*basis(:,:,n);
    V_out=V_out+spec(n,1)*vvpbasis(:,:,n);
    U_out=U_out+spec(n,1)*uvpbasis(:,:,n);
 
end


figure(2); clf;
subplot(2,1,1)
load coast;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
plotm(lat,long, 'k');hold on;
hold on;
pcolorm(YY*180/pi, XX*180/pi, ifield);
color_trans; colormap(cmap);
MAX=max(max(abs(field)));
caxis([-MAX MAX]);
colorbar;

div_out=sphdiv(U_out,V_out,yy*180/pi,xx*180/pi);

figure(2); hold on;
subplot(2,1,2)
load coast;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
plotm(lat,long, 'k');hold on;
hold on;
color_trans; colormap(cmap);
MAX=max(max(field_out));
caxis([-MAX MAX]);
pcolorm(YY*180/pi,XX*180/pi, field_out);
colorbar;




% figure(2); hold on;
% subplot(2,1,2)
% load coast;
% axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
% load coast
% plotm(lat,long, 'k');hold on;
% hold on;
% color_trans; colormap(cmap);
% %MAX=max(max(div_out));
% caxis([-MAX MAX]);
% pcolorm(YY*180/pi,XX*180/pi, div_out);
% colorbar;

[LON2,LAT2]=meshgrid(xx*180/pi, yy*180/pi);

lon3=[LON2(1,:) 360];
lat3=[-89;LAT2(:,1);89];
[LON3, LAT3]= meshgrid(lon3,lat3);


V_ext=[ [V_out(1,:) V_out(1,1)]; [V_out V_out(:,1)]; [V_out(end,:) V_out(end,1)]]; 
U_ext=[ [U_out(1,:) U_out(1,1)]; [U_out U_out(:,1)]; [U_out(end,:) U_out(end,1)]];

%interpolate back to the original grid
V_OUT=interp2(LON3,LAT3,V_ext, LON*180/pi, LAT*180/pi, 'spline');
U_OUT=interp2(LON3,LAT3,U_ext, LON*180/pi, LAT*180/pi, 'spline');

DIV_OUT=sphdiv(U_OUT, V_OUT, lat,lon);

sphcolor2(DIV_OUT*10^5/g, lat,lon,20);title('DIVERGENCE OF CORRECTED WINDS (kg/(m^2 s)') % the divergence of the corrected winds

%and the residual from the whole process
sphcolor2(DIV_OUT*10^5/g-MRES, lat,lon,21);title('THE RESIDUAL') % the divergence of the corrected winds


break
figure(3); clf
load coast;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
plotm(lat,long, 'k');hold on;
hold on; quiverm(LAT2,LON2,U_out,V_out)




