source: altifloat/matlab_toolbox/read_aviso_file.m @ 199

function [ output_args ] = read_aviso_file( aviso_netcdf_file_bck, aviso_netcdf_file_next, Nlat, Nlong, nbofdays )
%%%%%%%%%% This function is to read the AVISO netcdf files that are
%%%%%%%%%% downloaded from the Aviso website (aviso_netcdf_file).
%%%%%%%%%% The user should enter the name of the file of AVISO at the date 
%%%%%%%%%% needed (aviso_netcdf_file_bck) and the name of the file of 
%%%%%%%%%% Aviso a week later of the date (aviso_netcdf_file_next).
%%%%%%%%%% Note that you should not add the extension while entering the 
%%%%%%%%%% name of the file;
%%%%%%%%%% The user should also specify the dimention of the grid that  
%%%%%%%%%%  he wants (Nlat is for the latitude and Nlong is for the 
%%%%%%%%%% longitude). It then extracts the
%%%%%%%%%% horizontal u and vertical v terms of the velocity fields in the
%%%%%%%%%% dimentions needed and stores them in an output file that could
%%%%%%%%%% be read by the YAO software. This code gives you the velocity
%%%%%%%%%% fields at the 2 dates: the date needed and a week later.
%%%%%%%%%% The nbof dqys is the number of days between the 2 AVISO datas.

%%%%%%%%%% Usually use:
%%%%%%%%%% Nlat=58; %Phi
%%%%%%%%%% Nlong=87; %Lambda

write_u=1;

 
%for plotting


%The grid in degrees 
load 'obs_float/meshgrid2.dat';

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:



%UTRUE=the first week

U1=ncread(strcat('obs_aviso/',aviso_netcdf_file_bck,'.nc'),'Grid_0001');
uaviso1=U1(1:58,end-86:end)./100; %choose the eatern med and make it in m/s

V1=ncread(strcat('obs_aviso/',aviso_netcdf_file_bck,'.nc'),'Grid_0002');
vaviso1=V1(1:58,end-86:end)./100; %choose the eatern med and make it in m/s



%U the second week

U2=ncread(strcat('obs_aviso/',aviso_netcdf_file_next,'.nc'),'Grid_0001');
uaviso2=U2(1:58,end-86:end)./100;

V2=ncread(strcat('obs_aviso/',aviso_netcdf_file_next,'.nc'),'Grid_0002');
vaviso2=V2(1:58,end-86:end)./100;


%do the mask manually:
II1=find(isfinite(uaviso1)==0);
uaviso1(II1)=0;
II2=find(isfinite(uaviso2)==0);
uaviso2(II2)=0;

JJ1=find(isfinite(vaviso1)==0);
vaviso1(JJ1)=0;
JJ2=find(isfinite(vaviso2)==0);
vaviso2(JJ2)=0;


%interpolate to get the field every 2 h during the first week, in case need velocity in between:
Nt=nbofdays*24/2; %one week in units of delta t
for jt=1:Nt+1
    umod(:,:,jt)=(uaviso1*(Nt+1-jt)./Nt)+ (uaviso2*(jt-1)/Nt);
    vmod(:,:,jt)=(vaviso1*(Nt+1-jt)./Nt)+ (vaviso2*(jt-1)/Nt);
end
    





%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=umod(:,:,1)./delta_x.'; %u has units of pers second
vv=vmod(:,:,1)./delta_y.';

%Choose Uback=Uafter 3days

uu_pert=umod(:,:,35)./delta_x.'; 
vv_pert=vmod(:,:,35)./delta_y.';

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


%%%%%%%%%%%%%%%%%%%PUT U in YAO format%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

up_yao1=reshape(uu_pert, Nlong*Nlat,1);
vp_yao1=reshape(vv_pert, Nlong*Nlat,1);
%create the vectors rehsaped
cd1=[]; for k=1:Nlong; cd1=[cd1; k*ones(Nlat,1)]; end
cd2=repmat([1:1:Nlat].',Nlong,1);

Yao_u=[ ones(Nlong*Nlat,1)  cd1 cd2 u_yao1];
Yao_v=[ ones(Nlong*Nlat,1)  cd1 cd2 v_yao1];

Yao_up=[ ones(Nlong*Nlat,1)  cd1 cd2 up_yao1];
Yao_vp=[ ones(Nlong*Nlat,1)  cd1 cd2 vp_yao1];

if (write_u==1)

dlmwrite('obs_float/uzero.dat',Yao_u,'\t');
dlmwrite('obs_float/vzero.dat',Yao_v,'\t');
%save UV_0.mat uu vv;
dlmwrite('obs_float/unext3d.dat',Yao_up,'\t');
dlmwrite('obs_float/vnext3d.dat',Yao_vp,'\t');
%save UVnext3d.mat uu_pert vv_pert cc1 cc2;
end;


end