source: altifloat/matlab_toolbox/aviso2yao.m @ 160

function [Yao_uv,mask1] = aviso2yao( aviso_netcdf, meshmask,outfile,outmask,convlon)
% [Yao_uv,mask1] = aviso2yao( aviso_netcdf, meshmask,outfile,outmask)
% Read AVISO file aviso_netcdf_file and save to be read by read_bck (time step=1)
% Interpolate on meshmask (meshmask : two columns lon,lat)
% Column 1 : time step (=1)
% Column 2 : lon (Yi)
% Column 3 : lat (Yj)
% Column 4 : u value
% Column 5  : v value
% Values u and v are unitless
  

% Usually use:
% Nlat=58; %Phi
% Nlong=87; %Lambda
% aviso_netcdf='../../data/AVISO/dt_med_allsat_madt_uv_20130828_20140704.nc'

% meshmask='../obs_float/meshgrid3.dat'
% outfile='../obs_float/uv_bck_20130828.dat';
% outmask='../obs_float/mask_20130828.dat';

write_u=1;
if ~exist('convlon')
    convlon=true;
end

 
%for plotting


%The grid in degrees 
meshgrid2=load(meshmask);
Nlong=length(unique(meshgrid2(:,1)));
Nlat=length(unique(meshgrid2(:,2)));



GR_long=reshape(meshgrid2(:,1),Nlong,Nlat); %will give a Nlong by Nlat matrix 
%containing the long of the points 

GR_lat=reshape(meshgrid2(:,2),Nlong,Nlat); %will give a Nlong by Nlat matrix 
%containing the lat of the points 






%-----------The velocities U and V---------------------------------------%



%Load directly from AVISO folder where the u s for the whole mediterranean:



U=ncread(aviso_netcdf,'u');
V=ncread(aviso_netcdf,'v');
lon=ncread(aviso_netcdf,'lon');
lat=ncread(aviso_netcdf,'lat');

[LAT,LON]=meshgrid(lat,lon); %C'est à l'envers à cause de matlab

if convlon
%Conversion des longitudes
LON=LON-360;
end

%Conversion des vitesses inutile (malgré le champ scale_factor, ce que je trouve bizarre,
%mais les ordres de grandeurs semblent bon

Uint=interp2(LAT,LON,U,GR_lat,GR_long,'nearest');
Vint=interp2(LAT,LON,V,GR_lat,GR_long,'nearest');


mask_u=isfinite(Uint);
mask_v=isfinite(Vint);
mask=mask_u.*mask_v;

mask1=mask(:);

if (exist('outmask') & ~isempty(outmask))
dlmwrite(outmask,mask1,'precision',1);
end

%Masked velocity are zero
Uint(mask==0)=0;
Vint(mask==0)=0;


%change grid from degrees to uniteless

R_earth=6371229; %in meters
delta_x=zeros(Nlong,Nlat);
delta_y=zeros(Nlong,Nlat);

for ilon=1:Nlong-1;
    for jlat=1:Nlat-1;
        
    delta_x(ilon,jlat)=R_earth*(2*pi/360)*(GR_long(ilon+1,jlat)-GR_long(ilon,jlat))...
        .*cos(2*pi*GR_lat(ilon,jlat)/360); 
    delta_y(ilon,jlat)=R_earth*(2*pi/360)*(GR_lat(ilon,jlat+1)-GR_lat(ilon,jlat));  
    
    
end 
end



delta_x(Nlong,1:Nlat)=delta_x(Nlong-1,1:Nlat);
delta_y(Nlong,1:Nlat)=delta_y(Nlong-1,1:Nlat);


delta_x(1:Nlong,Nlat)=delta_x(1:Nlong,Nlat-1);
delta_y(1:Nlong,Nlat)=delta_y(1:Nlong,Nlat-1);



%choose UTRUE=Ufirst week

uu=Uint./delta_x; %u has units of pers second
vv=Vint./delta_y;

%Choose Uback=Uafter 3days


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%PUT U in YAO format%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Transpose everything


%create a Nlat*Nlong by 5 matrix first 
u_yao1=reshape(uu', Nlong*Nlat,1);
v_yao1=reshape(vv', Nlong*Nlat,1);

%cd1=[]; for k=1:Nlong; cd1=[cd1; k*ones(Nlat,1)]; end
%cd2=repmat([1:1:Nlat].',Nlong,1);

cd_lon=repmat([1:Nlong],Nlat,1);
cd_lat=repmat([1:Nlat]',1,Nlong);


Yao_uv=[ ones(Nlong*Nlat,1)  cd_lon(:) cd_lat(:)  u_yao1 v_yao1];

if (write_u==1 & exist('outfile') & ~isempty(outfile))

dlmwrite(outfile,Yao_uv,'\t');

end;


end