clear all; 
close all;


write_u=1;%if you want to write the velocities (true and background) to YAO DAT file
write_coefs=1;%if you want to write the filters coefs
%------------GRID-------------------------------------------------------%

%for plotting

Nlat=58; %Phi
Nlong=87; %Lambda
%The grid in degrees 
load '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('../obs_aviso/dt_ref_med_merged_madt_uv_20080423_20080423_20100503.nc','Grid_0001');
uaviso1=U1(1:58,end-86:end)./100; %choose the eatern med and make it in m/s

V1=ncread('../obs_aviso/dt_ref_med_merged_madt_uv_20080423_20080423_20100503.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('../obs_aviso/dt_ref_med_merged_madt_uv_20080430_20080430_20100503.nc','Grid_0001');
uaviso2=U2(1:58,end-86:end)./100;

V2=ncread('../obs_aviso/dt_ref_med_merged_madt_uv_20080430_20080430_20100503.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=7*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);
u_yao1=u_yao1
v_yao1=v_yao1
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/utrue.dat',Yao_u,'\t');
dlmwrite('../obs_float/vtrue.dat',Yao_v,'\t');
%save UV_0.mat uu vv;
dlmwrite('../obs_float/unext3d1.dat',Yao_up,'\t');
dlmwrite('../obs_float/vnext3d1.dat',Yao_vp,'\t');
%save UVnext3d.mat uu_pert vv_pert cc1 cc2;
end;

