source: trunk/src/cal_hcl.m @ 325

function [tps,alt,lat,lon,tpt,ep_cl]=cal_hcl(tps_min,tps_max,lat_min,lat_max,lon_min,lon_max)

%CAL_HCL Calcul de la hauteur de la couche limite ECMWF opera

%+
%
% .. _cal_hcl.m:
%
% =========
% cal_hcl.m
% =========
%
% DESCRIPTION
% ===========
%
% Calcul de la hauteur de la couche limite
% a partir des analyses operationnelles ECMWF
%
% Entrees :
%
% tps_min 
%    borne min pour la date au format matlab
% tps_max 
%    borne max pour la date au format matlab
% lat_min 
%    borne min pour la latitude
% lat_max 
%    borne max pour la latitude
% lon_min 
%    borne min pour la longitude
% lon_max 
%    borne max pour la longitude
%
% Sorties :
%
% tps
%    vecteur temps format matlab
% alt
%    vecteur altitude
% lat 
%    vecteur latitude
% lon
%    vecteur longitude
% tpt 
%    temperature potentielle
% ep_cl
%    hauteur de couche limite
%
% EXAMPLES
% ========
%
% ::
%
%   >> tps_min=datenum([2006 07 01])
%   >> tps_max=datenum([2006 07 10]);
%   >> lat_min=-5; lat_max=5;
%   >> lon_min=-10; lon_max=10;
%   >> [tps,alt,lat,lon,tpt,ep_cl]=cal_hcl(tps_min,tps_max,lat_min,lat_max,lon_min,lon_max);
%   >> moytps_epcl=squeeze(nanmean(ep_cl,1));
%   >> % Graphique
%   >> nb_color=20;
%   >> xmin=min(min(moytps_epcl));
%   >> xmax=max(max(moytps_epcl));
%   >> dx=(xmax-xmin)/nb_color;
%   >> figure;
%   >> contourf(lon,lat,moytps_epcl,[xmin:dx:xmax]);
%   >> colormap(jet(nb_color));
%   >> colorbar;
%
% SEE ALSO
% ========
% Appelle :
%
% :ref:`extract_nc.m` -> lecture fichier .nc
%
% :ref:`hspe2rmix.m`  -> calcul le rapport de mélange avec l'humidité spécifique
%
% :ref:`tpot.m`       -> calcul de la température potentielle
%
% :ref:`h_clim.m`     -> calcul de la hauteur de couche limite
%
% TODO
% ====
%
% pour l'instant fonction sur climserv donc matlab 7.10
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20100713T110004Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * header in ReStructured Text
%
% - mllod 20100713
%
%   * 1er commit
%
%-
[y_min m_min d_min hh_min mm_min ss_min]=datevec(tps_min);
[y_max m_max d_max hh_max mm_max ss_max]=datevec(tps_max);

%% Charge les analyses operationnelles ECMWF

tps_min_mod=(tps_min-datenum([1957 01 01 00 00 00]))*24;
tps_max_mod=(tps_max-datenum([1957 01 01 00 00 00]))*24;

if m_min == m_max

  [tps,level,lat,lon,ta]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/ta.2006%02d.amh.WESTAFR_025.nc',m_min),{'ta'},'time',tps_min_mod,tps_max_mod,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
  [tps,level,lat,lon,q]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/q.2006%02d.amh.WESTAFR_025.nc',m_min),{'q'},'time',tps_min_mod,tps_max_mod,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
  [tps,lat,lon,msl]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_SF/2006/msl.2006%02d.ash.WESTAFR_025.nc',m_min),{'msl'},'time',tps_min_mod,tps_max_mod,'lat',lat_min,lat_max,'lon',lon_min,lon_max);

else

  [tps,level,lat,lon,ta]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/ta.2006%02d.amh.WESTAFR_025.nc',m_min),{'ta'},'time',tps_min_mod,NaN,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
  [tps,level,lat,lon,q]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/q.2006%02d.amh.WESTAFR_025.nc',m_min),{'q'},'time',tps_min_mod,NaN,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
  [tps,lat,lon,msl]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_SF/2006/msl.2006%02d.ash.WESTAFR_025.nc',m_min),{'msl'},'time',tps_min_mod,NaN,'lat',lat_min,lat_max,'lon',lon_min,lon_max);

  i=1;
  while m_min+i ~= m_max

    [tps_tmp,level,lat,lon,ta_tmp]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/ta.2006%02d.amh.WESTAFR_025.nc',m_min+i),{'ta'},'time',NaN,NaN,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
    [tps_tmp,level,lat,lon,q_tmp]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/q.2006%02d.amh.WESTAFR_025.nc',m_min+i),{'q'},'time',NaN,NaN,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
    [tps_tmp,lat,lon,msl_tmp]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_SF/2006/msl.2006%02d.ash.WESTAFR_025.nc',m_min+i),{'msl'},'time',NaN,NaN,'lat',lat_min,lat_max,'lon',lon_min,lon_max);

    tps=[tps;tps_tmp];
    ta=[ta;ta_tmp];
    q=[q;q_tmp];
    msl=[msl;msl_tmp];

    i=i+1;
  end

  [tps_tmp,level,lat,lon,ta_tmp]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/ta.2006%02d.amh.WESTAFR_025.nc',m_min+i),{'ta'},'time',NaN,tps_max_mod,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
  [tps_tmp,level,lat,lon,q_tmp]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_ML/2006/q.2006%02d.amh.WESTAFR_025.nc',m_min+i),{'q'},'time',NaN,tps_max_mod,'level',61.,91.,'lat',lat_min,lat_max,'lon',lon_min,lon_max);
  [tps_tmp,lat,lon,msl_tmp]=extract_nc(sprintf('/bdd/OPERA/NETCDF/WESTAFR_025/4xdaily/AN_SF/2006/msl.2006%02d.ash.WESTAFR_025.nc',m_min+i),{'msl'},'time',NaN,tps_max_mod,'lat',lat_min,lat_max,'lon',lon_min,lon_max);


  tps=[tps;tps_tmp];
  ta=[ta;ta_tmp];
  q=[q;q_tmp];
  msl=[msl;msl_tmp];
end
tps=tps/24+datenum([1957 01 01 00 00 00]);

% Calcul des altitudes des niveaux de modele en passnat par la pression
p=flipud(load('/homedata/mleduc/OPERA/ml91.dat')); % col 1 : numero du niveau ; col 5 : pression a mi-couche
for i=1:length(level)+1
  phalf(:,i,:,:)=p(i,2)+p(i,3)*msl;
end
for i=1:length(level)
  p_mod(:,i,:,:)=squeeze(0.5*(phalf(:,i,:,:)+phalf(:,i+1,:,:)));
  alt(:,i,:,:)=(-8.32*((squeeze(ta(:,1,:,:))+squeeze(ta(:,i,:,:)))/2).*(log(squeeze(p_mod(:,i,:,:)))-log(msl))/(29*9.81))*1000;
end

% Calcul rapport de melange
r=hspe2rmix(q,p_mod);

% Calcul teperature potentielle
[tpt tv]=tpot(ta,p_mod,r);

% Calcul hauteur de couche limite
for i= 1:size(tpt,3)
  for j=1:size(tpt,4)
    [ep_cl(:,i,j),grad_tpt(:,i,j)]=h_clim(tpt(:,:,i,j),alt(:,:,i,j));
  end
end