Commit ade7be80 authored by dumoda01's avatar dumoda01

Ajout de scripts matlab de preparation des forcages

parent eabac40a
function narr_daily2gotm(year)
% cf. ftp://ftp.cdc.noaa.gov/Datasets/NARR/monolevel/
% cf. GOTM
% Dany DUMONT
% 17/06/2008: Creation
% 23/02/2011: Put year as an argument
% ______________________________________________________________________
yrstr = num2str(year);
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:1/24:init_date+nod-1/24]';
timevec = datevec(time);
tlength = length(time);
% ----------------------------------------------------------------------
% Momentum flux calculated using
% wind speed at 10m from NARR 3-hourly reanalysis
% ----------------------------------------------------------------------
ncfile = ['uwnd_',yrstr,'_narr_amdgulf.nc'];
uwnd = nc_varget(ncfile,'UWND'); % m s^-1
ncfile = ['vwnd_',yrstr,'_narr_amdgulf.nc'];
vwnd = nc_varget(ncfile,'VWND'); % m s^-1
% moyenne journali�re des donn�es
window = 8;
uwnd = filter(ones(1,window)/window,1,uwnd);
vwnd = filter(ones(1,window)/window,1,vwnd);
% interpolation des donn�es aux heures
uwnd = interp(uwnd,3);
vwnd = interp(vwnd,3);
% Wind stress calculation
% cf. Stewart, R.H. (2002) Introduction to Physical Oceanography
rho = 1.3; % density of air (kg m^-3)
U = sqrt(uwnd.^2 + vwnd.^2); % velocity amplitude (m s^-1)
C10 = nan .* ones(length(U),1);
for i = 1:length(U)
if (U(i) >= 3 && U(i) < 6)
C10(i) = (0.29 + 3.1./U(i) + 7.7./U(i).^2).*1e-3;
elseif U(i) >= 6
C10(i) = (0.60 + 0.070.*U(i)).*1e-3; % Yelland and Taylor (1996)
else
C10(i) = 2.18e-3;
end
end
% C10 = (0.75 + 0.065.*U).*1e-3; % Garratt (1977)
dir = atan2(vwnd,uwnd); % wind direction (rad)
tau = rho.*C10.*U.*U; % wind stress amplitude (N m^-2)
taux = tau.*cos(dir); % zonal wind stress (N m^-2)
tauy = tau.*sin(dir); % meridional wind stress (N m^-2)
% figure; plot(U,C10,'ob')
% figure; plot(U,tau,'ok')
clear U dir tau C10 rho
% Creation of the file 'narr_momentumflux.dat'
% file with tx and ty given in N/m^2
momentumflux = nan .* ones(tlength,8);
momentumflux(:,1:6) = timevec;
momentumflux(:,7) = taux;
momentumflux(:,8) = tauy;
fid = fopen(['narr_daily_momentumflux_',yrstr,'.dat'], 'wt');
fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
momentumflux');
fclose(fid);
% ----------------------------------------------------------------------
% Incoming shortwave radiative flux and surface heat flux
% from NARR 3-hourly reanalysis
% ----------------------------------------------------------------------
ncfile = ['dswrf_',yrstr,'_narr_amdgulf.nc'];
dswrf = nc_varget(ncfile,'DSWRF'); % W/m^2
% moyenne journali�re des donn�es
dswrf = filter(ones(1,window)/window,1,dswrf);
% interpolation des donn�es � chaque heure
dswrf = interp(dswrf,3);
a = find(dswrf < 0.0);
dswrf(a) = 0.0;
ncfile = ['uswrf_',yrstr,'_narr_amdgulf.nc'];
uswrf = nc_varget(ncfile,'USWRF'); % W/m^2
% moyenne journali�re des donn�es
uswrf = filter(ones(1,window)/window,1,uswrf);
% interpolation des donn�es aux heures
uswrf = interp(uswrf,3);
ncfile = ['dlwrf_',yrstr,'_narr_amdgulf.nc'];
dlwrf = nc_varget(ncfile,'DLWRF'); % W/m^2
% moyenne journali�re des donn�es
dlwrf = filter(ones(1,window)/window,1,dlwrf);
% interpolation des donn�es � chaque heure
dlwrf = interp(dlwrf,3);
ncfile = ['ulwrf_',yrstr,'_narr_amdgulf.nc'];
ulwrf = nc_varget(ncfile,'ULWRF'); % W/m^2
% moyenne journali�re des donn�es
ulwrf = filter(ones(1,window)/window,1,ulwrf);
% interpolation des donn�es � chaque heure
ulwrf = interp(ulwrf,3);
ncfile = ['shtfl_',yrstr,'_narr_amdgulf.nc'];
shtfl = nc_varget(ncfile,'SHTFL'); % W/m^2
% moyenne journali�re des donn�es
shtfl = filter(ones(1,window)/window,1,shtfl);
% interpolation des donn�es � chaque heure
shtfl = interp(shtfl,3);
ncfile = ['lhtfl_',yrstr,'_narr_amdgulf.nc'];
lhtfl = nc_varget(ncfile,'LHTFL'); % W/m^2
% moyenne journali�re des donn�es
lhtfl = filter(ones(1,window)/window,1,lhtfl);
% interpolation des donn�es � chaque heure
lhtfl = interp(lhtfl,3);
% Transmitted shortwave radiative flux
% that include the extinction due to the ice cover
% (Data from the Canadian Ice Service)
eval(['load cis_icec_',yrstr])
swr = (dswrf - uswrf).*(1-icec);
% No ice effect
% swr = (dswrf - uswrf);
figure
t=time-init_date;
plot(t,dswrf-uswrf,'LineWidth',2,'Color',[0.5 0.5 0.5]); hold on
plot(t,swr,'-k','LineWidth',2); hold off
xlim([1 nod]);
xlabel(['day of year ',yrstr])
ylabel('W/m^2')
h = legend('Over ice','Under ice')
legend(h,'boxoff')
saveas(gcf,'swr_icec.tiff', 'tiffn')
% Total heat flux (excluding the shortwave)
% positive = towards the ocean surface
% qt = dlwrf - ulwrf + shtfl + lhtfl; (????)
% that include the insulation due to the ice cover
% (Data from the Canadian Ice Service)
qt = (dlwrf - ulwrf + shtfl + lhtfl) .* (1-icec);
% No ice effect
% qt = dlwrf - ulwrf + shtfl + lhtfl;
% Creation of the file 'heatflux.dat'
% file with swr and total heat flux in W/m^2
init_date = datenum([year 01 01 00 00 00]);
time = [init_date:1/24:init_date+nod-1/24]';
timevec = datevec(time);
tlength = length(time);heatflux = nan .* ones(tlength,8);
heatflux(:,1:6) = timevec;
heatflux(:,7) = swr;
heatflux(:,8) = qt;
fid = fopen(['narr_daily_heatflux_ice_',yrstr,'.dat'], 'wt');
fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
heatflux');
fclose(fid);
function narr_hourly2gotm(year)
% cf. ftp://ftp.cdc.noaa.gov/Datasets/NARR/monolevel/
% cf. GOTM
% Uses the netcdf toolbox
% Dany DUMONT
% 14/05/2008: Creation
% 23/02/2011: Put year as an argument
% ______________________________________________________________________
yrstr = num2str(year)
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:1/24:init_date+nod-1/24]';
timevec = datevec(time);
tlength = length(time);
% ----------------------------------------------------------------------
% Momentum flux calculated using
% wind speed at 10m from NARR 3-hourly reanalysis
% ----------------------------------------------------------------------
ncfile = ['uwnd_',yrstr,'_narr_amdgulf.nc'];
uwnd = nc_varget(ncfile,'UWND');
ncfile = ['vwnd_',yrstr,'_narr_amdgulf.nc'];
vwnd = nc_varget(ncfile,'VWND');
% interpolation des donn�es � chaque heure
uwnd = interp(uwnd,3);
vwnd = interp(vwnd,3);
% Wind stress calculation
% cf. Stewart, R.H. (2002) Introduction to Physical Oceanography
rho = 1.3; % density of air (kg m^-3)
U = sqrt(uwnd.^2 + vwnd.^2); % velocity amplitude (m s^-1)
C10 = nan .* ones(length(U),1);
for i = 1:length(U)
if (U(i) >= 3 && U(i) < 6)
C10(i) = (0.29 + 3.1./U(i) + 7.7./U(i).^2).*1e-3;
elseif U(i) >= 6
C10(i) = (0.60 + 0.070.*U(i)).*1e-3; % Yelland and Taylor (1996)
else
C10(i) = 2.18e-3;
end
end
% C10 = (0.75 + 0.065.*U).*1e-3; % Garratt (1977)
dir = atan2(vwnd,uwnd); % wind direction (rad)
tau = rho.*C10.*U.*U; % wind stress amplitude (N m^-2)
taux = tau.*cos(dir); % zonal wind stress (N m^-2)
tauy = tau.*sin(dir); % meridional wind stress (N m^-2)
% figure; plot(U,C10,'ob')
% figure; plot(U,tau,'ok')
clear U dir tau C10 rho
% Creation of the file 'narr_momentumflux.dat'
% file with tx and ty given in N/m^2
momentumflux = nan .* ones(tlength,8);
momentumflux(:,1:6) = timevec;
momentumflux(:,7) = taux;
momentumflux(:,8) = tauy;
fid = fopen(['narr_hourly_momentumflux_',yrstr,'.dat'], 'wt');
fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
momentumflux');
fclose(fid);
% ----------------------------------------------------------------------
% Incoming shortwave radiative flux and surface heat flux
% from NARR 3-hourly reanalysis
% ----------------------------------------------------------------------
ncfile = ['dswrf_',yrstr,'_narr_amdgulf.nc'];
dswrf = nc_varget(ncfile,'DSWRF'); % W/m^2
% interpolation des donn�es � chaque heure
dswrf = interp(dswrf,3);
a = find(dswrf < 0.0);
dswrf(a) = 0.0;
ncfile = ['uswrf_',yrstr,'_narr_amdgulf.nc'];
uswrf = nc_varget(ncfile,'USWRF'); % W/m^2
% interpolation des donn�es � chaque heure
uswrf = interp(uswrf,3);
ncfile = ['dlwrf_',yrstr,'_narr_amdgulf.nc'];
dlwrf = nc_varget(ncfile,'DLWRF'); % W/m^2
% interpolation des donn�es � chaque heure
dlwrf = interp(dlwrf,3);
ncfile = ['ulwrf_',yrstr,'_narr_amdgulf.nc'];
ulwrf = nc_varget(ncfile,'ULWRF'); % W/m^2
% interpolation des donn�es � chaque heure
ulwrf = interp(ulwrf,3);
ncfile = ['shtfl_',yrstr,'_narr_amdgulf.nc'];
shtfl = nc_varget(ncfile,'SHTFL'); % W/m^2
% interpolation des donn�es � chaque heure
shtfl = interp(shtfl,3);
ncfile = ['lhtfl_',yrstr,'_narr_amdgulf.nc'];
lhtfl = nc_varget(ncfile,'LHTFL'); % W/m^2
% interpolation des donn�es � chaque heure
lhtfl = interp(lhtfl,3);
% Transmitted shortwave radiative flux
% that include the extinction due to the ice cover
% (Data from Canadian Ice Service)
eval(['load cis_icec_',yrstr])
swr = (dswrf - uswrf).*(1.0-icec);
for ii = 1:length(swr)
if swr(ii) < 0.0;
swr(ii) = 0.0;
end
end
% No ice effect
% swr = (dswrf - uswrf);
% Total heat flux (minus the shortwave)
% positive = towards the ocean surface
% qt = dlwrf - ulwrf + shtfl + lhtfl; (????)
% that include the insulation due to the ice cover
% (Data from D. Barber group)
qt = (dlwrf - ulwrf + shtfl + lhtfl) .* (1.0-icec);
% No ice effect
% qt = dlwrf - ulwrf + shtfl + lhtfl;
figure; plot(time,qt)
figure; plot(time,swr)
% Creation of the file 'heatflux.dat'
% file with swr and total heat flux in W/m^2
init_date = datenum([year 01 01 00 00 00]);
time = [init_date:1/24:init_date+nod-1/24]';
timevec = datevec(time);
tlength = length(time);heatflux = nan .* ones(tlength,8);
heatflux(:,1:6) = timevec;
heatflux(:,7) = swr;
heatflux(:,8) = qt;
fid = fopen(['narr_hourly_heatflux_ice_',yrstr,'.dat'], 'wt');
fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
heatflux');
fclose(fid);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment