source: altifloat/obs_external/read_davane.m

clear all;

%------------GRID-------------------------------------------------------%
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 


%--------------------TIME-------------------------------------------------%
dt=2*3600; %in seconds small delta t corresponding to Nume Scheme for advection

%---------------------U---------------------------------------------------%
load 'velbck.dat';
u=reshape(velbck(:,1),Nlong,Nlat,[]);
v=reshape(velbck(:,2),Nlong,Nlat,[]);
u1=u(:,:,1);
v1=v(:,:,1);


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);




u_pers=u1./delta_x; %u pers second
v_pers=v1./delta_y;



%note in fortran he does the nonuniform grid:

%rad = 2. * asin( 1. ) / 180. this is pi/180=2*pi/360
%     rae = 6371229.
%      DO jj = 1, jpj-1
%      DO ji = 1, jpi-1
%     e1t(ji,jj) = rae*rad*cos(rad*phi(ji,jj))*(lam(ji+1,jj)-lam(ji,jj))
%     e2t(ji,jj) = rae*rad                    *(phi(ji,jj+1)-phi(ji,jj))
%      END DO
%      END DO




uu=u_pers.';
vv=v_pers.';

%-----------------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);
%create the vectors rehsaped
c1=[]; for k=1:Nlong; c1=[c1; k*ones(Nlat,1)]; end
c2=repmat([1:1:Nlat].',Nlong,1);

%Yao_u=[ones(Nlong*Nlat,1) ones(Nlong*Nlat,1)  c1 c2 u_yao1];
%Yao_v=[ones(Nlong*Nlat,1) ones(Nlong*Nlat,1)  c1 c2 v_yao1];
Yao_u=[ ones(Nlong*Nlat,1)  c1 c2 u_yao1];
Yao_v=[ ones(Nlong*Nlat,1)  c1 c2 v_yao1];

%dlmwrite('../obs_float/uzero.dat',Yao_u,'\t');
%dlmwrite('../obs_float/vzero.dat',Yao_v,'\t');



% xx=R_earth*(2*pi/360)*GR_long.*cos(0.5*(lat_min+lat_max)); % this is an approximation 
% %that gives a "uniform" rectangular grid but at these scales, 
% %this approximation is not so good! see fig 1 
% 
% xx1=R_earth*(2*pi/360)*GR_long.*cos(2*pi*GR_lat/360); %this gives 
% %the non uniform grid that is not a rectangle
% 
% yy=R_earth*(2*pi/360)*GR_lat;