function sfa = ENRR_read_surfa_ecmwf(surfa_nc_path,stjd,endjd)

%{
    read surfa fields for ENRR 2016
    
    input start and end day-of-year

    returns structure 'sfa' of model fields
%}

%% initialize

% set DOY range (subtract 1 since ECMWF forecast initializes 12 hr earlier)
stjd = stjd-1; endjd = endjd-1;  

tdk = 273.15;

%% generate cell array of file names and compute timestamp
jd_vec = [];
for ii = stjd:endjd  % assuming 8 files per day @ 3 hr intervals
    jd_vec = [jd_vec;ii;ii;ii;ii;ii;ii;ii;ii];
end

Nfiles = length(jd_vec); % number of files to load
dfl = cell(Nfiles,1);
jd = NaN(Nfiles,1);
for ii = 1:8:Nfiles
    [m,d] = yd2md(2015,jd_vec(ii)); mstr = sprintf('%02i',m); dstr = sprintf('%02i',d);
    dfl{ii} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_12.nc']);
    dfl{ii+1} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_15.nc']);
    dfl{ii+2} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_18.nc']);
    dfl{ii+3} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_21.nc']);
    dfl{ii+4} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_24.nc']);
    dfl{ii+5} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_27.nc']);
    dfl{ii+6} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_30.nc']);
    dfl{ii+7} = fullfile(surfa_nc_path,['ECMWF_2016',mstr,dstr,'12_33.nc']);
    jd(ii:ii+7) = jd_vec(ii:ii+7) + 12/24 + [12;15;18;21;24;27;30;33]/24;
end

%% load surfa data
% read grid units, lon & lat
lon = ncread(dfl{1},'lon'); % deg E
lat = ncread(dfl{1},'lat'); % deg N
A = length(lon); B = length(lat);

% preallocate empty arrays
rsi0 = NaN(A,B,Nfiles);
rli0 = NaN(A,B,Nfiles);
rsin = NaN(A,B,Nfiles);
rlin = NaN(A,B,Nfiles);
rsi = NaN(A,B,Nfiles);
rli = NaN(A,B,Nfiles);
hsi = NaN(A,B,Nfiles);
hli = NaN(A,B,Nfiles);
evap = NaN(A,B,Nfiles);
lowcld = NaN(A,B,Nfiles);
midcld = NaN(A,B,Nfiles);
tauni = NaN(A,B,Nfiles);
tauwi = NaN(A,B,Nfiles);
wcloud = NaN(A,B,Nfiles);
P = NaN(A,B,Nfiles);
ts = NaN(A,B,Nfiles);
tdew = NaN(A,B,Nfiles);
ta = NaN(A,B,Nfiles);
ux = NaN(A,B,Nfiles);
uy = NaN(A,B,Nfiles);

% loop through nc files to load individual fields for all dates/times
for ii = 1:Nfiles
    rsi0(:,:,ii) = ncread(dfl{ii},'Clear_sky_downward_solar_flux');             % W/m2 at surface
    rli0(:,:,ii) = ncread(dfl{ii},'Clear_sky_downward_long_wave_flux');         % W/m2 at surface
    rsin(:,:,ii) = ncread(dfl{ii},'Net_short_wave_radiation_surface');          % W/m2 at surface
    rlin(:,:,ii) = ncread(dfl{ii},'Net_long_wave_radiation_surface');           % W/m2 at surface
    rsi(:,:,ii) = ncread(dfl{ii},'Short_wave_radiation');                       % W/m2 at surface
    rli(:,:,ii) = ncread(dfl{ii},'Long_wave_radiation');                        % W/m2 at surface
    hsi(:,:,ii) = ncread(dfl{ii},'Sensible_heat_flux');                         % W/m2 at surface
    hli(:,:,ii) = ncread(dfl{ii},'Latent_heat_flux');                           % W/m2 at surface
    evap(:,:,ii) = ncread(dfl{ii},'Evaporation')*1e3;                           % g/m2
    lowcld(:,:,ii) = ncread(dfl{ii},'Low_level_cloud_cover');                   % percent
    midcld(:,:,ii) = ncread(dfl{ii},'Mid_level_cloud_cover');                   % percent
    tauni(:,:,ii) = ncread(dfl{ii},'Meridional_component_of_momentum_flux');    % N/m2
    tauwi(:,:,ii) = ncread(dfl{ii},'Zonal_component_of_momentum_flux');         % N/m2
    wcloud(:,:,ii) = ncread(dfl{ii},'Cloud_water')*1e3;                         % g/m2
    P(:,:,ii) = ncread(dfl{ii},'Pressure')/100;                                 % mb @ 0m
    ts(:,:,ii) = ncread(dfl{ii},'Soil_temperature')-tdk;                        % C @ 0m
    tdew(:,:,ii) = ncread(dfl{ii},'Dew_point_temperature')-tdk;                 % C @ 2m
    ta(:,:,ii) = ncread(dfl{ii},'Temperature')-tdk;                             % C @ 2m
    ux(:,:,ii) = ncread(dfl{ii},'u_wind');                                      % m/s @ 10m
    uy(:,:,ii) = ncread(dfl{ii},'v_wind');                                      % m/s @ 10m
end

u = sqrt(ux.^2 + uy.^2); % compute wind speed

qa = td2q(tdew,P);       % specific humidity @ 2m

%% take derviatives of accumulated flux variables

xx = find(mod(jd,1)==0); % find surfa integral reset indices - 0000Z each day

taunx = diff(tauni,1,3)/3600/3;         % time interval is 3 hrs
taun = cat(3, tauni(:,:,1)/3600/12, taunx);
taun(:,:,xx) = tauni(:,:,xx)/3600/12;   % first value each day is 12 hr integral
tauwx = diff(tauwi,1,3)/3600/3;
tauw = cat(3, tauwi(:,:,1)/3600/12, tauwx);
tauw(:,:,xx) = tauwi(:,:,xx)/3600/12;
tau = sqrt(taun.^2 + tauw.^2);
clear taunx tauwx

hsx = diff(hsi,1,3)/3600/3;
hs = cat(3,hsi(:,:,1)/3600/12, hsx);
hs(:,:,xx) = hsi(:,:,xx)/3600/12;
hs = -hs;  % sign convention
clear hsx

hlx = diff(hli,1,3)/3600/3;
hl = cat(3,hli(:,:,1)/3600/12, hlx);
hl(:,:,xx) = hli(:,:,xx)/3600/12;
hl = -hl;  % sign convention
clear hlx

rsx = diff(rsi,1,3)/3600/3;
rs = cat(3,rsi(:,:,1)/3600/12, rsx);
rs(:,:,xx) = rsi(:,:,xx)/3600/12;
clear rsx

rlx = diff(rli,1,3)/3600/3;
rl = cat(3,rli(:,:,1)/3600/12, rlx);
rl(:,:,xx) = rli(:,:,xx)/3600/12;
clear rlx

rs0x = diff(rsi0,1,3)/3600/3;
rs0 = cat(3,rsi0(:,:,1)/3600/12, rs0x);
rs0(:,:,xx) = rsi0(:,:,xx)/3600/12;
clear rs0x

rl0x = diff(rli0,1,3)/3600/3;
rl0 = cat(3,rli0(:,:,1)/3600/12, rl0x);
rl0(:,:,xx) = rli0(:,:,xx)/3600/12;
clear rl0x

rsnx = diff(rsin,1,3)/3600/3;
rsn = cat(3,rsin(:,:,1)/3600/12, rsnx);
rsn(:,:,xx) = rsin(:,:,xx)/3600/12;
clear rsnx

rlnx = diff(rlin,1,3)/3600/3;
rln = cat(3,rlin(:,:,1)/3600/12, rlnx);
rln(:,:,xx) = rlin(:,:,xx)/3600/12;
clear rlnx

evpx = diff(evap,1,3)/3600/3;
evp = cat(3,evap(:,:,1)/3600/12, evpx);
evp(:,:,xx) = evap(:,:,xx)/3600/12;
clear evpx


%% create output structure

sfa = struct; % initialize output struture

vars = {'dfl','jd','lat','lon','evp','lowcld','midcld','wcloud','P','ts','tdew','ta','ux','uy','u',...
        'tau','hs','hl','rs','rl','rs0','rl0','rsn','rln','qa'};

for ii = 1:length(vars)
    name = vars{ii};
    sfa = setfield(sfa, name, eval(name));
end

end