function [x_grid,y_grid,z_grid,u_grid,v_grid,w_grid,ws_grid,XLAT,XLONG,PSFC,varargout]=wrf_windfield(input_data, varargin )

if length(input_data.u_var)==0
    input_data.u_var='U_1';
    input_data.v_var='V_1';
    input_data.w_var='W_1';
    input_data.p_var='PHB';
    input_data.lat_var='XLAT';
    input_data.long_var='XLONG';
    input_data.psfc='PSFC';
end

path1=input_data.path1;
file_number=input_data.file_number;

path0=(['C:\Documents and Settings\sylvie.lorsolo\My Documents\research\HRD research\',path1]); 
filename=dir([path0,'\*.nc']);
ncload([path0,'\',filename(file_number).name],input_data.u_var,input_data.v_var,input_data.w_var,input_data.p_var, input_data.lat_var, input_data.long_var,input_data.psfc,varargin{:});
U=eval(input_data.u_var);
V=eval(input_data.v_var);
W=eval(input_data.w_var);
PHB=eval(input_data.p_var);
XLAT=eval(input_data.lat_var);
XLONG=eval(input_data.long_var);
PSFC=eval(input_data.psfc);

filename(file_number).name

clear('input_data');




[li,col]=size(XLAT);
temp_p=reshape(PSFC,li*col,1);
[center_pressure,center_grid]=(min(temp_p));
lat0=XLAT(center_grid)*ones(li*col,1);
lon0=XLONG(center_grid)*ones(li*col,1);
xlati=reshape(XLAT,li*col,1);
xlongi=reshape(XLONG,li*col,1);
range=distdim(distance(lat0, lon0, xlati, xlongi),'deg','m');
azi=azimuth(lat0, lon0, xlati, xlongi);
azi2=azi;
uu=find(azi>=270 & azi <360);
azi2(uu)=360-azi(uu)+90;
uu=find(azi>=0 & azi <270);
azi2(uu)=90-azi2(uu);
azi=azi2;
[x,y]=pol2cart(deg2rad(azi),range);
xx=reshape(x,li,col);
yy=reshape(y,li,col);
g= 9.780327*(1 + .0053024*(sin(deg2rad(XLAT))).^2 - .0000058*(sin(2*deg2rad(XLAT))).^2 );


U=permute(U,[2 3 1]);
V=permute(V,[2 3 1]);
W=permute(W,[2 3 1]);
PHB=permute(PHB,[2 3 1]);

[li,col,lev]=size(U);
u_grid=NaN*zeros(li,li,lev);
v_grid=NaN*zeros(li,li,lev);
w_grid=NaN*zeros(li,li,lev);
geop=NaN*zeros(li,li,lev);

for j=1:li
    u_grid(:,j,:)=0.5*(U(:,j,:)+U(:,j+1,:));
end

for i=1:li
    v_grid(i,:,:)=0.5*(V(i,:,:)+V(i+1,:,:));
end

for k=1:lev
    w_grid(:,:,k)=0.5*(W(:,:,k)+W(:,:,k+1));
    geop(:,:,k)=0.5*(PHB(:,:,k)+PHB(:,:,k+1));
end

ws_grid=sqrt(u_grid.^2 + v_grid.^2);

nout=max(nargout,1)-10;    

for nn=1:nout
    temp=eval(varargin{nn});
    
    dim_temp=ndims(temp);
    if dim_temp >2
        temp_final=NaN*zeros(li,li,lev);
        temp=permute(temp,[2 3 1]);
        [li_temp,col_temp,lev_temp]=size(temp);
        if li_temp > li 
            for i=1:li
                
                temp_final(i,:,:)=0.5*(temp(i,:,:)+temp(i+1,:,:));
            end
            
        elseif col_temp > li
            
            for j=1:li
                temp_final(:,j,:)=0.5*(temp(:,j,:)+temp(:,j+1,:));
            end
        elseif lev_temp > lev
            
            for k=1:lev
                temp_final(:,:,k)=0.5*(temp(:,:,k)+temp(:,:,k+1));
            end
        else
            temp_final=temp;
        end
    else
        [li_temp,col_temp]=size(temp);
        if li_temp > li
            temp_final=NaN*zeros(li,li);
            for i=1:li
                
                temp_final(i,:)=0.5*(temp(i,:)+temp(i+1,:));
            end
        elseif col_temp > li 
            temp_final=NaN*zeros(li,li);
            for j=1:li
                temp_final(:,j)=0.5*(temp(:,j)+temp(:,j+1));
            end
        elseif col_temp == lev+1
            temp_final=NaN*zeros(li_temp,lev);
            for k=1:lev
                temp_final(:,k)=0.5*(temp(:,k)+temp(:,k+1));
            end
        else
            temp_final=temp;
        end
    end
    varargout(nn)={temp_final};
end

[li,col,z_level]=size(geop);
g1=repmat(g,[1 1 z_level]);
z_grid=geop./g1;
x_grid= repmat(xx,[1 1 z_level]); %linspace(-221250,221250,L);;%0:20:4180;  %0:10:2090;%-5200:50:5250;%400:50:10450;
y_grid= repmat(yy,[1 1 z_level]); %linspace(-221250,221250,L);%0:20:4180;   %0:10:2090;%-5200:50:5250;%-10400:100:10500;
%z_grid=1:z_level;
%[x_3dgrid,y_3dgrid,z2_3dgrid]=meshgrid(x_grid,y_grid,z_grid);
