%{
Reader of 1hr Flux database
%}
clear all;

%% loads data and identify variables
% x=importdata('/Volumes/LACIE_SHARE/FluxDataBase/FluxBase/FluxBase_1hr_v09162016.txt');
y=importdata('c:\data\cwf\dropbox\doc\climate\oceobs\flux_base\FluxBase_1hr_v12032018.txt');
ii=find(y(:,55)>10); % select dynamo
x=y(ii,:);
jdyx=x(:,1);%julian day at beginning of time average
ushp=x(:,2);%doppler log, SCS (m/s)
U=x(:,3);%true wind,PSD sonic (m/s)
dir=x(:,4);%true wind direction, PSD sonic (deg)
urel=x(:,5);%relative wind speed, PSD (m/s)
reldir=x(:,6);%relative wind direction (from),clockwise rel ship's bow, PSD sonic (deg)
head=x(:,7);%ship heading, deg clockwise rel north, SCS laser ring gyro
tsnk=x(:,8);%sea snake temperature, PSD, 0.05 m depth (C)
ta=x(:,9);%air temperature, PSD (C)%
qse=x(:,10);%sea surface specific humidity, from snake (g/kg)
%qse=qsee([tsnk press]);
qa=x(:,11);%air specific humidity, PSD (g/kg)
hsc=x(:,12);%sensible heat flux, covariance, PSD sonic anemometer(W/m^2)
hsib=x(:,13);%sensible heat flux, ID, PSD sonic anemometer(W/m^2)
hsb=x(:,14);%bulk sensible heat flux, (W/m^2)
x(:,15)=x(:,15)*1.00;%acount for licor-sonic separation
hlc=x(:,15);%latent heat flux, covariance, (W/m^2)
hlib=x(:,16);%latent heat flux, ID, (W/m^2)
hlb=x(:,17);%bulk latent heat flux, W/m^2 (includes Webb et al. correction)
taucx=x(:,18);%covariance streamwise stress, PSD sonic anemometer (N/m^2)
taucy=x(:,19);%covariance cross-stream stress, PSD sonic anemometer (N/m^2)
tauib=x(:,20);%ID streamwise stress, PSD sonic anemometer (N/m^2)
taub=x(:,21);%bulk wind stress along mean wind, (N/m^2)
rs=x(:,22);%downward solar flux, PSD units (W/m^2)
rl=x(:,23);%downward IR flux, PSD units (W/m^2)
org=x(:,24);%rainrate, PSD STI optical rain gauge, uncorrected (mm/hr)
J=x(:,25);%ship plume contamination index
tiltx=x(:,26);%flow tilt at PSD sonic anemometer, earth frame
Jm=x(:,27);%ship maneuver index
ct=x(:,28);%ct^2 (K^2/m^.667)
cq=x(:,29);%cq^2 ((g/kg)^2/m^.667)
cu=x(:,30);%cu^2 ((m/s)^2/m^.667)
cw=x(:,31);%cw^2 ((m/s)^2/m^.667)
hrain=x(:,32);%rain heat flux,Gosnell et al 1995, JGR, 18437-18442, (W/m^2)
hlwebb=x(:,33);%correction to measured latent heat flux, Webb et al.
lat=x(:,34);%latitude, deg (SCS pcode)
lon=x(:,35);%longitude, deg (SCS pcode)
zu=x(:,36);%height of mean wind sensor
zt=x(:,37);%height of mean air temperature sensor
zq=x(:,38);%height of mean air humidity sensor
%***** ships imet and scs data
sog=x(:,39);%speed over ground, SCS gps, (m/s)
U_scs=x(:,40); %true wind speed, imet propvane anemometer (m/s)
dir_scs=x(:,41);%true wind direction (from),clockwise rel north, imet,(deg)
cog=x(:,42);%%course over ground, SCS gps, (m/s)
tsg=x(:,43);%tsg water temperature, 5 m depth, (C)
ta_im=x(:,44);%imet air temperature (C)
qs_tsg=x(:,45);%imet bulk water specific humidity (g/kg)
qa_im=x(:,46);%imet air specific humidity, (g/kg)
press=x(:,47);
RH=x(:,48);
qlic=x(:,49);     % Specific Humidity from Licor (g/kg)
rlcler=x(:,50);   % Longwave downwelling clear sky model (W/m2)
Rscl=x(:,51);     % Shortwave solar radiance clear sky model (W/mccc2)
rnl=x(:,52);      % net longwave heat flux
hnet=x(:,53);     % net heat flux, bulk, W/m2 hourly
year=x(:,54);
pcode=x(:,55);    %project code;  
% Metz =  77; Epic = 69; WHOTS = 87; STRATUS = 83; NEAQS = 78; HiWInGS = 72; GasEx = 71
% DYNAMO = 68; Calwater = 67; Capricorn = 73;

% uuj=x(:,56);  %variance of u (m/s)
% vvj=x(:,57);  %variance of v (m/s)
% wwj=x(:,58);  %variance of w (m/s)
% tt=x(:,59);   %variance of t (degC)
% qq=x(:,60);   %std of q  (g/kg)

jdy=jdyx+datenum(year,12,31)-datenum(min(year),12,31);
jdyloc=jdy+lon/15/24;

%% fratini correction
ffc=(1.58+1.95*log(qa))/10;
   ii=find(~isnan(qa)& ~isnan(qlic) & pcode~=68);
   hlc(ii)=hlc(ii)+(qa(ii)-qlic(ii))./max(1,qa(ii)).*ffc(ii).*hlc(ii);
   hlib(ii)=hlib(ii)+(qa(ii)-qlic(ii))./max(1,qa(ii)).*ffc(ii).*hlib(ii);


%% do a quick map of all tracks
Lonmin = -180;       % longitudes are 0-360, going E
Lonmax = 180;
Latmin = -80;        % deg N is positive
Latmax = 80;
LatLonLim = [Lonmin Lonmax Latmin Latmax];
Marker=3;

figure;
m_proj('mercator','lon',[LatLonLim(1) LatLonLim(2)],'lat',[LatLonLim(3) LatLonLim(4)]);
%m_grid('box','fancy','tickdir','in');
m_grid('box','on','tickdir','in');
m_coast('linewidth',2,'color',0.5+[0 0 0]);
m_coast('patch',[.7 .7 .7],'edgecolor','none');
hw09=m_line(lon,lat,'color','r','linestyle','none','marker','.','markersize',Marker);


%% filtering
rdir_lo     = 300;	% usually -60 to 60 but BB adds sometimes +360 to negative values
rdir_hi     = 60;
rain_lim    = 5;    % rain rate, mm/hr
ushp_lim    = 3.5;  % ship speed, m/s


%% Selects DYNAMO data for plot
toto = find (pcode==68);
hw09=m_line(lon(toto),lat(toto),'color','g','linestyle','none','marker','.','markersize',Marker);


%% Plot the transfer coefficients
Plot_TransferCoeff2

