function out = SURFA_avg(stjd,surfa_nc_path,lon,lat,conv_drops)
% computes means for SURFA fields within a polygon
% defined by vectors lon(conv_drops) and lat(conv_drops).
% out = [u,ta,rh,P,ts,rs,rl,qa,hs,hl];

% first, read SURFA data fields for the given day of year
sfa = ENRR_read_surfa_ecmwf(surfa_nc_path,stjd,stjd);

% define  rough lat/lon bounds for ENRR flight area and
% find SURFA indices corresponding to these limits
lat_lo = -2;        lat_lo_idx = find(sfa.lat==lat_lo);
lat_hi = 20;        lat_hi_idx = find(sfa.lat==lat_hi);
lon(lon<0) = lon(lon<0)+360;
lon_lo = floor(min(lon))-1;   lon_lo_idx = find(sfa.lon==lon_lo);
lon_hi = ceil(max(lon))+1;   lon_hi_idx = find(sfa.lon==lon_hi);

% define vertices of polygon from location of dropsondes
% that encircle the area of interest
Xlon = lon(conv_drops)'; Ylat = lat(conv_drops)';

% duplicate fisrt vertex at end of array to complete a polygon
Xlon = [Xlon; Xlon(1)]; Ylat = [Ylat; Ylat(1)];

% build array of SURFA lat/lon indices for the rough flight region
Nlon = length(sfa.lon(lon_lo_idx:lon_hi_idx)); Nlat = length(sfa.lat(lat_hi_idx:lat_lo_idx));
SFAlon = repmat(sfa.lon(lon_lo_idx:lon_hi_idx),Nlat,1);
SFAlat = NaN(Nlon*Nlat,1);
for ii = 1:Nlat
    SFAlat(Nlon*(ii-1)+1:Nlon*ii) = sfa.lat(lat_hi_idx+ii-1);
end

% find SURFA indices that fall within the study area polygon
idx = inpolygon(SFAlon,SFAlat,Xlon,Ylat);
lon_valid = SFAlon(idx); lat_valid = SFAlat(idx); N = length(lon_valid);
lon_idx = NaN(N,1); lat_idx = NaN(N,1);
for ii = 1:N  % SURFA indices for points falling within polygon
    lon_idx(ii) = find(sfa.lon==lon_valid(ii));
    lat_idx(ii) = find(sfa.lat==lat_valid(ii));
end

% now, compute means of grid cells within polygon

% sensible heat, W/m2
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.hs(lon_idx(ii),lat_idx(ii),1);
end
hs = mean(temp);

% latent heat, W/m2
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.hl(lon_idx(ii),lat_idx(ii),1);
end
hl = mean(temp);

% SST
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.ts(lon_idx(ii),lat_idx(ii),1);
end
ts = mean(temp);

% T-air
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.ta(lon_idx(ii),lat_idx(ii),1);
end
ta = mean(temp);

% q-air
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.qa(lon_idx(ii),lat_idx(ii),1);
end
qa = mean(temp);

% rh
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.rh(lon_idx(ii),lat_idx(ii),1);
end
rh = mean(temp);

% wind speed
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.u(lon_idx(ii),lat_idx(ii),1);
end
u = mean(temp);

% pressure
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.P(lon_idx(ii),lat_idx(ii),1);
end
P = mean(temp);

% surface downwelling solar radiation
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.rs(lon_idx(ii),lat_idx(ii),1);
end
rs = mean(temp);

% surface downwelling IR
temp = NaN(N,1);
for ii = 1:N
    temp(ii) = sfa.rl(lon_idx(ii),lat_idx(ii),1);
end
rl = mean(temp);

out = [u,ta,rh,P,ts,rs,rl,qa,hs,hl];

end