ecmwf2gotm.m 3.27 KB
function ecmwf2gotm(region,file,lon,lat,year)
% ECMWF2GOTM - Format meteorological variables from ECMWF for GOTM. It uses
% data archived on brandypot.uqar.ca. Variables are:
%
%     T2M       2 meter temperature (degK)
%     U10       10 meter zonal wind component (m/s)
%     V10       10 meter meriodional wind component (m/s)
%     D2M       2 meter dewpoint temperature (degK)
%     MSL       Mean sea level pressure (Pa)
%     TCC       Total cloud cover (0-1)
%
% Syntax:  ecmwf2gotm(region,file,lon,lat,year)
%
% Inputs:
%
%     region     'pc'  - Potter Cove
%                'ap'  - Antarctic Peninsula and Weddell Sea
%                'gsl' - Golfe Saint-Laurent
%                'now' - North Water
%                'amg' - Amundsen Gulf
%                'gsj' - Golfo San Jorge
%     file       String to be usd in the output file name
%     lon        longitude in degrees_east (0-360)
%     lat        latitude in degrees_north (-90-90)
%     year       year to be extracted
%
% Outputs:
%
%    output1 - Description
%    output2 - Description
%
% Example: 
%    ecmwf2gotm('gsl',
%
% Other m-files required: Netcdf Toolbox
% Subfunctions: none
% MAT-files required: none
%
% See also: NARR2GOTM, NCEP2GOTM, ICEC2GOTM
%
% cf. data-portal.ecmwf.int/data/d/interim_daily/
% cf. GOTM
% Uses the netcdf toolbox

% Author: Dany Dumont
% email: dany_dumont@uqar.ca
% Website: http://www.ismer.ca/dumont-dany
% May 2011; Last revision: 15-October-2014
% ______________________________________________________________________

yrstr = num2str(year);
disp(yrstr)

if mod(year,400) == 0
   nod = 366;
elseif mod(year,100) == 0
   nod = 365;
elseif mod(year,4) == 0
   nod = 366;
else
   nod = 365;
end

% Time scale given in Matlab serial date number
init_date = datenum([year 01 01 00 00 00]);
time      = (init_date:6/24:init_date+nod-6/24)';
timevec   = datevec(time);
%tlength   = length(time);

datadir = ['/sas/usagers/share_lasso/data/ecmwf/',region,'/6hourly/'];

ncfile = [datadir,'ecmwf_',region,'_meteo_',yrstr,'.nc'];
hlon = double(nc_varget(ncfile,'longitude'));
hlat = double(nc_varget(ncfile,'latitude'));
u10  = nc_varget(ncfile,'u10');
v10  = nc_varget(ncfile,'v10');
msl  = nc_varget(ncfile,'msl');
t2m  = nc_varget(ncfile,'t2m');
d2m  = nc_varget(ncfile,'d2m');
tcc  = nc_varget(ncfile,'tcc');

% Unit conversion
t2m  = t2m - 273.15;     % degK to degC
d2m  = d2m - 273.15;     % degK to degC
msl  = msl./100;         % Pa to hPa

% Extract index position
tol  = 0.8;
i    = find(hlon < lon + tol & hlon > lon - tol,1);
j    = find(hlat < lat + tol & hlat > lat - tol,1);
if isempty(i) || isempty(j)
    disp('** Error : Can''t find the specified lat/lon location.        ')
    disp('**         Verify that you have specified the longitude in the')
    disp('**         [0-360] degree interval (e.g. 296 instead of -64). ')
    return
end

% Extract one point and form the array
u10  = squeeze(u10(:,j,i));
v10  = squeeze(v10(:,j,i));
msl  = squeeze(msl(:,j,i));
t2m  = squeeze(t2m(:,j,i));
d2m  = squeeze(d2m(:,j,i));
tcc  = squeeze(tcc(:,j,i));

array = cat(2,timevec,u10,v10,msl,t2m,d2m,tcc);

% Write the file in the GOTM format
fid  = fopen(['ecmwf_',file,'_meteo_',yrstr,'.dat'], 'w');
fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%3.2f\t%3.2f\t%4.1f\t%3.2f\t%3.2f\t%1.2f\n', ...
    array');
fclose(fid);