clear all;
close all;
set_machine_et;

oceandir = [datadir 'DYNAMO/'];
load([oceandir 'DYNAMO_data_from_EThompson_July2020.mat']);

version_num = 'v1';
ncfilename =  [oceandir 'DYNAMO_2011_flux_met_sea_stable_10min_' version_num '.nc'];
matfilename =  [oceandir 'DYNAMO_2011_flux_met_sea_stable_10min_' version_num '.mat'];

pathreadme = oceandir;

%% vars I don't know... and am not saving
% not sure what dsea is? not in the readme.
% Smax doesn't make sense compar
% leave out SSS... it's not surface.

%% rename vars from met/flux file
fre = {'time';'Lat';'Lon';'yday';
'Tsea';'SOG';'COG';'stress';'shf';'lhf';'B_Qnet';...
'Qsea';'Cp';'sigDir';'sigH';'interped';...
'Heading';'TseaTSG';'SalTSG';'wdirR';...
'Pair10';'Ur10';'T02';'RH02';'Q02';'Q10';'RH10';'rhf';...
'IRdn';'IRup';'Solardn';'E';'P'; 'thetav_02';...
'Precip';'Evap';'Sclr';'Smax';'Solarup';'dSolar';...
'wspdR';'Qnet';'dTair';'dqair';'tw';'ustar';'rho0';'deltaT';...
'supp_ID';'dist_ID';'acti_ID';'wwb_ID';...
'hS';'hS_hat';'d_cham'};

fredone = {'t0';'lat';'lon';'jd';
'Tsnk';'sog';'cog';'tau';'hs';'hl';'B_Heat';...
'qsnk';'cp';'cdir';'Hsig';'U_flag';...
'hed';'Ttsg_s';'Stsg_s';'wdir_sfc';...
'P10';'U10_sfc';'T2';'rh2';'q2';'q10';'rh10';'hr';...
'lw_dn';'lw_up';'sw_dn';'erate';'prate';'thetav2';....
'paccum';'eaccum';'sw_dn_clr';'sw_dn_no_atm';'sw_up';'sw_dn_diffuse';...
'wspd_sfc';'hnet';'dT';'dq';'Tw';'usr';'SSRHO';'dT_rain';...
'in_supp';'in_dist';'in_acti';'in_wwb';...
'hStable';'hStable_hat';'depth'};

ab4 = a;

for i = 1:length(fre)
    a.(fredone{i}) = a.(fre{i});
end
a = rmfield(a,fre);

%%% fix time... rounding???
[year, month, day, hour, minute, second] = datevec(a.t0);
a.t = datenum(year, month, day, hour, roundn(minute,1), 0);
[a.year, a.month, a.day, a.hour, a.minute, second] = datevec(a.t);

% time elapsed in seconds since Jan 1 of that year
% note... don't add this to a structure. just call it time. used below in a
% special netcdf way. 

a.time = (a.t - datenum(2011,1,1,0,0,0))*24*3600;

% dimensions of files
ntime = length(a.t);
ndepth = length(a.depth);
nrl = a.nrl;
ndwl = a.ndwl;
nstable = a.nstable;
nboth = a.nboth;

%%% redo signs of fluxes: positive into ocean; same for air-sea gradients
a.sw_dn_clr = -a.sw_dn_clr;
a.sw_dn_no_atm = -a.sw_dn_no_atm;
a.sw_dn_diffuse = -a.sw_dn_diffuse;
a.sw_up = -a.sw_up;
a.lw_up = -a.lw_up;
a.sw_dn = -a.sw_dn;
a.lw_dn = -a.lw_dn;
a.hnet = -a.hnet;
a.hl = -a.hl;
a.hs = -a.hs;
a.hr = -a.hr;
a.dT = -a.dT;
a.dq = -a.dq;

% reset paccum and eaccum to 0 at beginning of each cruise
octset = find(a.t < datenum(2011,11,1,0,0,0));
novset = find(a.t > datenum(2011,11,1,0,0,0));
a.paccum(octset) = a.paccum(octset) - a.paccum(1);
a.paccum(novset) = a.paccum(novset) - a.paccum(novset(1));
a.eaccum(octset) = a.eaccum(octset) - a.eaccum(1);
a.eaccum(novset) = a.eaccum(novset) - a.eaccum(novset(1));


% redo IDs to be 0/1 not nan/1
a.in_stable(isnan(a.in_stable) == 1) = 0;
a.in_clear(isnan(a.in_clear) == 1) = 0;
a.in_dwl(isnan(a.in_dwl) == 1) = 0;
a.in_rl(isnan(a.in_rl) == 1) = 0;
a.in_both(isnan(a.in_both) == 1) = 0;
a.in_orl(isnan(a.in_orl) == 1) = 0;
a.in_odwl(isnan(a.in_odwl) == 1) = 0;
a.in_supp(isnan(a.in_supp) == 1) = 0;
a.in_dist(isnan(a.in_dist) == 1) = 0;
a.in_acti(isnan(a.in_acti) == 1) = 0;
a.in_wwb(isnan(a.in_wwb) == 1) = 0;

a.stable_dwl(isnan(a.stable_dwl) == 1) = 0;
a.stable_rl(isnan(a.stable_rl) == 1) = 0;
a.stable_both(isnan(a.stable_both) == 1) = 0;
a.stable_start(isnan(a.stable_start) == 1) = 0;
a.stable_end(isnan(a.stable_end) == 1) = 0;

%%% flip dimensions
a.dwl_start     = a.dwl_start';
a.dwl_end       = a.dwl_end';
a.dwl_dur       = a.dwl_dur';
a.rl_start      = a.rl_start';
a.rl_end        = a.rl_end';
a.rl_dur        = a.rl_dur';
a.both_start    = a.both_start';
a.both_end      = a.both_end';
a.both_dur      = a.both_dur';

% make IDs out of wh_

a.in_mjo = zeros(1,ntime); a.in_mjo([a.octtime a.novtime]) = 1;

a.in_supp1= zeros(1,ntime); a.in_supp1(a.wh_supp1) = 1;
a.in_supp2= zeros(1,ntime); a.in_supp2(a.wh_supp2) = 1;
a.in_dist1= zeros(1,ntime); a.in_dist1(a.wh_dist1) = 1;
a.in_dist2= zeros(1,ntime); a.in_dist2(a.wh_dist2) = 1;
a.in_acti1= zeros(1,ntime); a.in_acti1(a.wh_acti1) = 1;
a.in_acti2= zeros(1,ntime); a.in_acti2(a.wh_acti2) = 1;
a.in_wwb1 = zeros(1,ntime); a.in_wwb1(a.wh_wwb1) = 1;
a.in_wwb2 = zeros(1,ntime); a.in_wwb2(a.wh_wwb2) = 1;
a.in_wwb3 = zeros(1,ntime); a.in_wwb3(a.wh_wwb3) = 1;
a.in_octtime = zeros(1,ntime); a.in_octtime(a.octtime) = 1;
a.in_novtime = zeros(1,ntime); a.in_novtime(a.novtime) = 1;

%       wh_supp_begin_RL: [1 19]
%       wh_dist_begin_RL: [2 3 4 5 6 7 8 9 10 20 21 22 23 24 25 26]
%       wh_acti_begin_RL: [11 12 13 14 15 16 17 27 28 29 30 31 32 33]
%        wh_wwb_begin_RL: [18 34]
%      wh_supp_begin_DWL: [1 2 3 4 5 6 7 8 19 20 21 22 23]
%      wh_dist_begin_DWL: [9 10 11 12 13 24 25 26 27 28]
%      wh_acti_begin_DWL: [14 15 16 17 18 29 30 31]
%       wh_wwb_begin_DWL: []
%     wh_supp_begin_BOTH: [1 11]
%     wh_dist_begin_BOTH: [2 3 4 5 6 12 13 14 15]
%     wh_acti_begin_BOTH: [7 8 9 10 16 17]

a.dwl_start_supp= zeros(1,ndwl); a.dwl_start_supp(a.wh_supp_begin_DWL) = 1; %		index 1 = when warm layer started in suppressed period
a.dwl_start_dist= zeros(1,ndwl); a.dwl_start_dist(a.wh_dist_begin_DWL) = 1; %		index 1 = when warm layer started in disturbed period
a.dwl_start_acti= zeros(1,ndwl); a.dwl_start_acti(a.wh_acti_begin_DWL) = 1; %		index 1 = when warm layer started in active period
a.dwl_start_wwb= zeros(1,ndwl); a.dwl_start_wwb(a.wh_wwb_begin_DWL) = 1; %		index 1 = when warm layer started in WWB period

a.rl_start_supp= zeros(1,nrl); a.rl_start_supp(a.wh_supp_begin_RL) = 1;	%       index 1 = when rain layer started in suppressed period
a.rl_start_dist= zeros(1,nrl); a.rl_start_dist(a.wh_dist_begin_RL) = 1;	%       index 1 = when rain layer started in disturbed period
a.rl_start_acti= zeros(1,nrl); a.rl_start_acti(a.wh_acti_begin_RL) = 1;	%       index 1 = when rain layer started in active period
a.rl_start_wwb= zeros(1,nrl); a.rl_start_wwb(a.wh_wwb_begin_RL) = 1;	%           index 1 = when rain layer started in WWB period

a.both_start_supp= zeros(1,nboth); a.both_start_supp(a.wh_supp_begin_BOTH) = 1; %		index 1 = when combo rain+warm layer started in suppressed period
a.both_start_dist= zeros(1,nboth); a.both_start_dist(a.wh_dist_begin_BOTH) = 1; % 	index 1 = when combo rain+warm layer started in disturbed period
a.both_start_acti= zeros(1,nboth); a.both_start_acti(a.wh_acti_begin_BOTH) = 1; %		index 1 = when combo rain+warm layer started in active period
a.both_start_wwb= zeros(1,nboth); a.both_start_wwb(a.wh_wwb_begin_BOTH) = 1; %		index 1 = when combo rain+warm layer started in WWB period


%% add vars
a.slp = a.P10 + 0.125*10;



%% pad original Cham data with a zero depth value
a.S = vertcat(nan(1,ntime), a.S);
a.SIG = vertcat(nan(1,ntime), a.SIG);
a.T = vertcat(nan(1,ntime), a.T);
a.dSdz = vertcat(nan(1,ntime), a.dSdz);
a.dSIGdz = vertcat(nan(1,ntime), a.dSIGdz);
a.dTdz = vertcat(nan(1,ntime), a.dTdz);


%% fields to save
%%% leave out cog -- way too noisy because speed was slow
%%% leave out cdir -- all nan (?)
keep = {...
'dwl_start';'dwl_end';'dwl_dur';'rl_start';'rl_end';'rl_dur';...
'both_start';'both_end';'both_dur';'stable_start';'stable_end';'stable_dur';...
'stable_dwl';'stable_rl';'stable_both';...
'dwl_start_supp';'dwl_start_dist';...
'dwl_start_acti';'dwl_start_wwb';...
'rl_start_supp';'rl_start_dist';...
'rl_start_acti';'rl_start_wwb';...
'both_start_supp';'both_start_dist';...
'both_start_acti';'both_start_wwb';...
'in_stable';'in_clear';'in_dwl';'in_rl';'in_both';'in_orl';'in_odwl';...
'in_supp';'in_dist';'in_acti';'in_wwb';'in_mjo';...
'in_supp1';'in_supp2';'in_dist1';'in_dist2';'in_acti1';'in_acti2';...
'in_wwb1';'in_wwb2';'in_wwb3';'in_octtime';'in_novtime';...
'depth';'time';'year';'month';'day';'hour';'minute';...
'S';'T';'SIG';'dSdz';'dTdz';'dSIGdz';'alpha';'beta';...
'N2T_smooth';'N2S_smooth';'N2_smooth';'N2T';'N2S';'N2';...
'PR';'EV';'hStable';'hStable_hat';'hTC';'hB';'hZ';'US';'B';'B_Salt';'B_Heat';'B_Train';...
'Tw';'dT_rain';'thetav2';'c_sw';'Fatm';'slp';...
'SST';'SSRHO';'jd';'sw_dn_diffuse';'Tsnk';'sog';'wdir';'tau';'hs';'hl';...
'qsnk';'Hsig';'cp';'cdir';'cspd';'U_flag';'hed';'Ttsg_s';'Stsg_s';'wdir_sfc';...
'U10';'P10';'U10_sfc';'T2';'rh2';'q2';'T10';'q10';'rh10';'hr';...
'SSQ';'lw_dn';'lw_up';'sw_dn';'erate';'prate';'paccum';'eaccum';...
'sw_dn_clr';'sw_dn_no_atm';'zenith';'sw_up';...
'wspd';'wspd_sfc';'hnet';'dT';'dq';'Pcombo_1min';'rmm1';'rmm2';'mjo_amp';'mjo_phase'};

clear q;
for i = 1:length(keep)
   q.(keep{i}) = a.(keep{i}); 
end

q = orderfields(q);
%%% redo field names in reordered position
f_q = fields(q);
nq = length(f_q);

[lf, uf, lnf] = dynamo_assign_long_info(f_q, nq);
for j = 1:nq
    if isempty(lf{j}) == 1
        disp(['missing info for field ' sprintf('%i',j) ' = ' string(f_q(j))]);
    end
    if isempty(uf{j}) == 1
        disp(['missing unit for field ' sprintf('%i',j) ' = ' string(f_q(j))]);
    end
    if isempty(lnf{j}) == 1
        disp(['missing long name for field ' sprintf('%i',j) ' = ' string(f_q(j))]);
    end
end

save(matfilename,'q');

%% readme
% readme_f_narrative  = [pathreadme '_narrative_' cruise '_' version_num '.txt'];
% readme_full_f10     = [pathreadme '_readme_' cruise '_' version_num '.txt'];
% readme_data_f10     = [pathreadme '_fields_' cruise '_' version_num '.txt'];
% header_data_f10     = [pathreadme '_header_' cruise '_' version_num '.txt'];
% fileID_d_f10 = fopen(readme_data_f10,'w');
% fileID_h_f10 = fopen(header_data_f10,'w');
% 
% for j = 1:nq
%     formatSpec = 'i(%i) = %15s %-60s\n';
%     fprintf(fileID_d_f10, formatSpec, j, string(f_q(j)), string(lf(j)) );
%     fprintf(fileID_h_f10,'%20s',string(f_q(j)) );
% end
% fclose(fileID_d_f10);
% fclose(fileID_h_f10);
% 
% command_f10_a= ['cat ' readme_f_narrative ' ' readme_data_f10 ' > ' readme_full_f10];
% system(command_f10_a);


%% netcdf file

thedata = q;
ff = fields(thedata);
nf = length(ff);

%==========================================================================
% prep the data (from ODM save style)
%==========================================================================
% ff = fieldnames(DSD);
% ff = vars_to_save;

%%% if present, get rid of single-value or readme like metadata fields.
%%% These should be saved as attributes not variables.


vartypes = cell(nf,1);
drow = nan(nf,1);
dcol = nan(nf,1);
for i = 1:nf
    eval(['thisvar = thedata.' ff{i} ';']);
    varname = char(ff(i));
    vartype = class(thisvar);
    vartypes(i) = {vartype};
    [drow(i), dcol(i)] = size(thisvar);
    disp([' saving ' varname ' = ' vartype]);
%     nccreate(dsdfile,varname,'Datatype',vartype,'Dimensions',{'r',drow,'c',dcol});
%    ncwrite(dsdfile,varname,thisvar)
end

d2 = find(drow == ndepth & dcol == ntime);
dt = find(drow == 1 & dcol == ntime);
dd = find(drow == 1 & dcol == ndepth);
ds = find(drow == 1 & dcol == nstable);
dr = find(drow == 1 & dcol == nrl);
dw = find(drow == 1 & dcol == ndwl);
db = find(drow == 1 & dcol == nboth);
nd = length(d2) + length(dt) + length(dd) + length(ds) + length(dr) + length(dw) + length(db);
disp([' out of ' strnd(nf,0) ' total vars, ' strnd(nd,0) ' are accounted for']);

%%% redo all of this so that its in order by dimension
%%% these are the variables we are saving to nc
fff = [ff(dd); ff(d2); ff(dt); ff(ds); ff(dr); ff(dw); ff(db)]; 
% nf = length(fff);
vartypesf = [vartypes(dd); vartypes(d2); vartypes(dt); vartypes(ds); vartypes(dr); vartypes(dw); vartypes(db)]; 

ddrow = nan(nf,1);
ddcol = nan(nf,1);
for i = 1:nf
    eval(['thisvar = thedata.' fff{i} ';']);
    [ddrow(i), ddcol(i)] = size(thisvar);
end
dd2 = find(ddrow == ndepth & ddcol == ntime);
ddt = find(ddrow == 1 & ddcol == ntime);
ddd = find(ddrow == 1 & ddcol == ndepth);
dds = find(ddrow == 1 & ddcol == nstable);
ddr = find(ddrow == 1 & ddcol == nrl);
ddw = find(ddrow == 1 & ddcol == ndwl);
ddb = find(ddrow == 1 & ddcol == nboth);

nf2 = length(dd2);
nft = length(ddt);
nfd = length(ddd);
nfs = length(dds);
nfr = length(ddr);
nfw = length(ddw);
nfb = length(ddb);

%==========================================================================
% save this data to a netcdf file
%==========================================================================

% [nc dataID timeID] = nc_init_dsd(dsdfile,DSD.t,DSD.bin_sizes,units(d2),ff(d2),longnames(d2));

% initialize file:
nc = netcdf.create(ncfilename, 'clobber');


% create ID for (define) the dimensions
time_dimID=netcdf.defDim(nc,'time',ntime);
depth_dimID=netcdf.defDim(nc,'depth',ndepth);
stable_dimID=netcdf.defDim(nc,'stable_layers',nstable);
rl_dimID=netcdf.defDim(nc,'rain_layers',nrl);
dwl_dimID=netcdf.defDim(nc,'diurnal_warm_layers',ndwl);
both_dimID=netcdf.defDim(nc,'combo_rain_warm_layers',nboth);

theID = nan(nf,1);

%% create ID (define) for the variables using these dimensions

%%% time variables
for i = 1:nft
    j = ddt(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),time_dimID);
end

%%% depth variables
for i = 1:nfd
    j = ddd(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),depth_dimID);
end

%%% depth x time variables
for i = 1:nf2
    j = dd2(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),[depth_dimID time_dimID]);
end

%%% stable variables
for i = 1:nfs
    j = dds(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),stable_dimID);
end
%%% rl variables
for i = 1:nfr
    j = ddr(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),rl_dimID);
end
%%% dwl variables
for i = 1:nfw
    j = ddw(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),dwl_dimID);
end
%%% both variables
for i = 1:nfb
    j = ddb(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),both_dimID);
end

% variable names and units
[lf, uf, lnf] = dynamo_assign_long_info(fff, nf);
for j = 1:nf
    if isempty(lf{j}) == 1
        disp(['missing info for field ' sprintf('%i',j) ' = ' string(fff{j})]);
    end
    if isempty(uf{j}) == 1
        disp(['missing unit for field ' sprintf('%i',j) ' = ' string(fff{j})]);
    end
    if isempty(lnf{j}) == 1
        disp(['missing long name for field ' sprintf('%i',j) ' = ' string(fff{j})]);
    end
end

% define the units and long names
for i = 1:nf
    netcdf.putAtt(nc,theID(i),'units',char(uf(i)));
    netcdf.putAtt(nc,theID(i),'longname',char(lnf(i)));
    netcdf.putAtt(nc,theID(i),'info',char(lf(i)));
end

% define some global attributes
varid = netcdf.getConstant('GLOBAL');
netcdf.putAtt(nc,varid,'creation_date',datestr(now));
netcdf.putAtt(nc,varid,'creator','Elizabeth J. Thompson');
netcdf.putAtt(nc,varid,'contact_email','elizabeth.thompson@noaa.gov');
netcdf.putAtt(nc,varid,'please_acknowledge_1:','The PIs responsible for collecting ship met, flux, SST, navigation data: Chris Fairall PSL, Simon de Szoeke OSU, James Edson WHOI/UConn... with much help from Ludovic, Sergio, Byron');
netcdf.putAtt(nc,varid,'please_acknowledge_2:','The PI responsible for collecting ocean data: James Moum OSU... with much help from OSU Ocean Mixing Group including Aurelie Moulin');
netcdf.putAtt(nc,varid,'please_acknowledge_3:','The PI responsible for radar rain datasets: Steven Rutledge CSU');
netcdf.putAtt(nc,varid,'please_acknowledge_4:','The PI responsible for ocean stable layer analysis and interpolation of datasets: Elizabeth Thompson PSL, as part of her PhD at CSU (advisors: Steven Rutledge CSU, James Moum OSU). Aurelie Moulin OSU provided much help with data processing.');
netcdf.putAtt(nc,varid,'please_acknowledge_NOAA_PSL:','Chris Fairall, Sergio Pezoa, Elizabeth J. Thompson at NOAA Physical Sciences Laboratory');
netcdf.putAtt(nc,varid,'please_acknowledge_CIRES_CU_NOAA_PSL:','Ludovic Bariteau, Byron Blomquist at CIRES cooperative institute between University of Colorado and NOA Physical Sciences Laboratory');
netcdf.putAtt(nc,varid,'please_acknowledge_OSU:','Simon de Szoeke, James Moum, Aurelie Moulin at Oregon State University (Aurelie now at APL-UW)');
netcdf.putAtt(nc,varid,'please_acknowledge_WHOI:','James Edson currently at Woods Hole Oceanographic Institution (at University of Connecticut during experiment)');
netcdf.putAtt(nc,varid,'please_acknowledge_CSU:','Steven Rutledge at Colorado State University (E Thompson was also with CSU during experiment)');
netcdf.putAtt(nc,varid,'reference','This reference, and those therein, explain this dataset: Thompson, E. J., Moum, J. N., Fairall, C. W., & Rutledge, S. A. (2019). Wind limits on rain layers and diurnal warm layers. J. Geophysical Research: Oceans, 124, 897? 924. https://doi.org/10.1029/2018JC014130');
netcdf.putAtt(nc,varid,'sensor_info_1', 'zu = 17;       % NOAA bow mast sonic anemometer sensor height (m) ... for wind');
netcdf.putAtt(nc,varid,'sensor_info_2', 'zt = 18;       % NOAA bow mast aspirated air T sensor height (m)');
netcdf.putAtt(nc,varid,'sensor_info_3', 'zq = 18;       % NOAA bow mast aspirated air rh sensor height (m)');
netcdf.putAtt(nc,varid,'sensor_info_4', 'zsnk = 0.05;   % NOAA sea snake depth below water (m)');
netcdf.putAtt(nc,varid,'sensor_info_5','ztsg_ship = 5; % ship intake depth below water for thermosalinograph (TSG) T, S (m)');
netcdf.putAtt(nc,varid,'sensor_info_n','please read the readme and variable info attributes for more information');
netcdf.putAtt(nc,varid,'sign_conventions_1','wind directions are in meteorological convention: direction from');
netcdf.putAtt(nc,varid,'sign_conventions_2','current directions are in oceanographic convention: direction toward');
netcdf.putAtt(nc,varid,'sign_conventions_3','fluxes are positive if heating the ocean: lw_dn and sw_dn are positive by definition; hs is sometimes positive');
netcdf.putAtt(nc,varid,'sign_conventions_4','net heat flux positive heating the ocean: hnet = sw_dn + sw_up + lw_dn + lw_up + hs + hl + hr');
netcdf.putAtt(nc,varid,'sign_conventions_5','B surace buoyancy flux into ocean is positive when making ocean more buoyant');
netcdf.putAtt(nc,varid,'sign_conventions_6','Fatm surface buoyancy flux into atmosphere is positive when making atmosphere more buoyant');
netcdf.putAtt(nc,varid,'mjo_time_periods_1',['supp1 start ... end = ' datestr(a.supp1(1),0) ' ... ' datestr(a.supp1(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_2',['supp2 start ... end = ' datestr(a.supp2(1),0) ' ... ' datestr(a.supp2(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_3',['dist1 start ... end = ' datestr(a.dist1(1),0) ' ... ' datestr(a.dist1(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_4',['dist2 start ... end = ' datestr(a.dist2(1),0) ' ... ' datestr(a.dist2(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_5',['acti1 (contains WWB 1)    start ... end = ' datestr(a.acti1(1),0) ' ...' datestr(a.acti1(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_6',['acti2 (contains WWBs 2+3) start ... end = ' datestr(a.acti2(1),0) ' ...' datestr(a.acti2(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_7',['wwb1 start ... end = ' datestr(a.wwb1(1),0) ' ... ' datestr(a.wwb1(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_8',['wwb2 start ... end = ' datestr(a.wwb2(1),0) ' ... ' datestr(a.wwb2(2),0)]);
netcdf.putAtt(nc,varid,'mjo_time_periods_9',['wwb3 start ... end = ' datestr(a.wwb3(1),0) ' ... ' datestr(a.wwb3(2),0)]);
netcdf.putAtt(nc,varid,'note_1','Course over ground, cog, not offered because it was too noisy due to slow speeds. Use heading, hed, instead.');
netcdf.putAtt(nc,varid,'note_2','Current direction, cdir, had nans. Not yet available, but will be included in next version. From experience, current was from west. ');
netcdf.putAtt(nc,varid,'pr@_FillValue','NaN');

% % tell the netcdf file that you are done defining how it will work
netcdf.endDef(nc);
% 
% % write the data to the variables you just prepared.
for i = 1:nf
%    disp(['dataID = ' fff{i} 'ID;']);
%    disp(['data = DSD.' fff{i} ';']);
%    eval(['dataID = ' fff{i} 'ID;']);
   eval(['data = thedata.' fff{i} ';']);
   netcdf.putVar(nc,theID(i),data);
end

% 
% % IT IS VERY IMPORTANT TO CLOSE THE FILE!
netcdf.close(nc)

%% check the nc file
disp('----------------------------------');

ncload2(ncfilename);
% ncdisp(ncfilename);



% %% check with plots
% plot_it = 1;
% if plot_it == 1
%     
% figure;
% plot(q.t, q.prate,'o', q.t, q.Pcombo_1min,'.');
% grid on;
% legend('prate','Pcombo_1min');
% datetick('x','mm/dd');
% 
% 
% end