narr_daily2gotm.m 5.48 KB
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')

% 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);