source: trunk/src/readERA_LE.m @ 325

function[mat_3D,lat ,lon,jours,niveau]=readERA(file, z_deb, z_fin,nom_var,chemin)
%
%+
%
% ============
% readERA_LE.m
% ============
%
%
% .. function:: readERA(file, z_deb, z_fin,nom_var,chemin)
%
% DESCRIPTION
% ===========
%
% routine de lecture d un fichier de donnees ERA40 au format netcdf sur
% la zone definie par les points z_deb zfin de forme [lat lon time level]
% 
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20110705T125628Z cratos.locean-ipsl.upmc.fr (Linux)
%
%   * minimal header
%
%-
if(nargin<2)
z_deb=[ -30 -60 datenum(1957,1,1) 500]  
end
if(nargin<3)
z_fin=[ 45 50 datenum(2020,1,1) 500]   
end
if(nargin<4)
  [groupe,reste]=strtok(file,'.');
  [nom_var,reste]=strtok(reste,'.'); 
end
if(nargin<5)
     chemin='/homedata/eymard/OPERA_Trop/';
end
if(chemin(end)~='/')
chemin=[chemin,'/'];
end

file=[chemin,file];

nc=netcdf(file,'nowrite');
longitude=nc{'lon'}(:);
latitude=nc{'lat'}(:);
jours=nc{'time'}(:);
niveau=nc{'level'}(:);

jo=find(jours~=0);
jours=jours(jo)/24+datenum(1957,01,01); % En jours juliens
jo=find(jours>z_deb(3) & jours<z_fin(3));
jours=jours(jo);


lom=find(longitude-360>z_deb(2) & longitude-360<z_fin(2));
lop=find(longitude>z_deb(2) & longitude<z_fin(2));
lon=cat(1,longitude(lom)-360,longitude(lop));


%[lat,la]=sort(latitude,'ascend');
la=find(latitude>z_deb(1) & latitude<z_fin(1));
lat=sort(latitude(la),'ascend');

if(nom_var=='r')
ni=find(niveau>=z_deb(4) & niveau<=z_fin(4));

coef=[1 0.8 0.4 0.4];  % RH  la meilleure corr 2006 pour B5
coef=[0.6 0.8 1 0.8];  % RH  la meilleure corr 2006 pour B3 (mais moins sensible que B5)
%coef=[0 1 1 0];  % pour verif
coef
else
ni=[min(ni) max(ni)];  % pour geopt
end
niveau=niveau(ni);

 
variable_p=[];
variable_m=[];
 
if(isempty(niveau)) %variable 3D 
%extraction de la variable pour les longitudes positives et negatives
if(~isempty(lop))
  variable_p=nc{nom_var}(jo,la,lop);
end
if(~isempty(lom))
  variable_m=nc{nom_var}(jo,la,lom);
end
variable=cat(3,variable_m,variable_p); 
%Remise des latitudes dans le sens ascendant
variable=variable(:,end:-1:1,:);
%% mat_3D est la matrice (lat,lon, jours) de sortie 
mat_3D=permute(variable,[2 3 1]);  


 else %Variable 4D dont on moyenne la dimension "niveau" (nom_var=r) ou la difference (geopt)
%extraction de la variable pour les longitudes positives et negatives
if(~isempty(lop)) 
  variable_p=nc{nom_var}(jo,ni,la,lop);
end
if(~isempty(lom))
 variable_m=nc{nom_var}(jo,ni,la,lom);
end
variable=cat(4,variable_m,variable_p);
%Remise des latitudes dans le sens ascendant
variable=variable(:,:,end:-1:1,:);

%% mat_3D est la matrice (lat,lon,jours) de sortie 
if(nom_var=='r')
nom_var
  mat_3D_lev=permute(variable,[3 4 1 2]);
%mat_3D=squeeze(mean(mat_3D_lev,4));
mat_3D=squeeze(mat_3D_lev(:,:,:,1)*coef(1));
for i=2:size(niveau)
  mat_3D=mat_3D+mat_3D_lev(:,:,:,i)*coef(i);
end
 else
nom_var
   mat_3D=permute(variable,[2,3,4,1]);
mat_3D=squeeze(diff(mat_3D));

end
end