source: altifloat/obs_external/read_davane.m @ 88

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 Nlat by Nlong matrix 
%containing the long of the points 

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


lat_min=2*pi*min(min(GR_lat))/360;
lat_max=2*pi*max(max(GR_lat))/360;

%transforming from degrees to metres with an approximation of avergaging the
%cosine since latitude doesnt seem to change much in region of interest

R_earth=6371229; %in meters

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;

figure(1);
plot(xx,yy,'ro')
hold on
plot(xx1,yy,'k+')

%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


dx=xx(1,2)-xx(1,1)
dy=yy(2,1)-yy(1,1)


%--------------------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).'.*dt/dx;
v1=v(:,:,1).'.*dt/dx;

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

%create a Nlat*Nlong by 5 matrix first 
u_yao1=reshape(u1, Nlong*Nlat,1);
v_yao1=reshape(v1, 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');