Ajout de ecmwf2gotm.m, un script qui formate les donnees de re-analyse ECMWF...

Ajout de ecmwf2gotm.m, un script qui formate les donnees de re-analyse ECMWF en fichier meteo pour GOTM (airsea.nml).

scripts/matlab/ecmwf2gotm.m

+function ecmwf2gotm(region,file,lon,lat,year)
+% Format meteorological variables from ECMWF for GOTM.
+%
+% Variables:
+%
+%     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)
+%
+% Function arguments
+%
+%     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
+% 
+% cf. data-portal.ecmwf.int/data/d/interim_daily/
+% cf. GOTM
+% Uses the netcdf toolbox
+
+% Dany DUMONT
+% 14/05/2011: Creation
+% 24/01/2013: Standardized and adapted to the data file organization
+% ______________________________________________________________________
+
+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/dumoda01/shared/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 not 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);