% program to save level 2 ODM disdrometer output form SPURS-2, PISTON, and
% TGT datasets

% history / authors
%%%     --- code updated Oct 2020 and run by E. Thompson
%              Elizabeth Thompson
%              NOAA Physical Sciences Lab
%              elizabeth.thompson@noaa.gov

%%%     --- code originally developed by K. Drushka
%              Kyla Drushka
%              Applied Physics Lab, University of Washington
%              kdrushka@apl.uw.edu
%              see nc_init.m and nc_init_dsd.m


clear all;
close all;
set_machine_et;



% for f = 1:3
for f = 1:2
    if f == 1
        cruise = 'PISTON'; yrstr = '2018'; 
        odm_location_on_ship = 'on top of pilot house, starboard side';
    elseif f == 2
        cruise = 'SPURS2'; yrstr = '2016-2017';
        odm_location_on_ship = '2016: on top of forward O2 deck van; 2017: bow foremast near 17 m';
    elseif f == 3
        cruise = 'TGT'; yrstr = '2017-2018';
        odm_location_on_ship = 'on top of pilot house, starboard side';
    end
    
cruise_str = [cruise '_' yrstr];

thepath = ['/Users/eliz/Dropbox/DSD/' cruise '/final/'];

version_num = 'v1';


% load latest DSD data file (might be updating if on the cruise currently)
these_files = dir([thepath cruise '_dsd_level_2_' version_num '*.mat']);
ncfilename = [thepath cruise '_' yrstr '_dsd_level_2_' version_num '.nc'];
matfilename = [thepath cruise '_' yrstr '_dsd_level_2_' version_num '.mat'];
files = char(these_files(:).name);
[nfiles, nchars] = size(files);
load([thepath files(end,:)]);
disp(['latest DSDN file: ' files(end,:)]);

% load fortran data... but threshold it with same wh_good
if strcmp(cruise,'SPURS2') == 1
    kfile_2 = '/Users/eliz/Dropbox/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_mod.mat';
    if isfile(kfile_2) == 0
        k_s16 = load('/Users/eliz/Dropbox/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_2016.mat');
        k_s17 = load('/Users/eliz/Dropbox/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_2017.mat');
        k = cat_2mats(k_s16.dsd, k_s17.dsd);
        k.README = [k_s16.dsd.README ; k_s17.dsd.README];
    else
        load(kfile_2)    
    end    
    DSDN.Pfortran = k.p(DSDN.ID_good_gamma == 1);
else
    DSDN.Pfortran = DSDN.t*nan;
end


vars_to_save = {'bad_data';'good_data';'good_gamma';'t';'year';'month';'day';'hour';'minute';...
    'Pfortran';'P';'z';'Z';'W';'Pd';'zd';'Zd';'Wd';'voltage';'urel_anem';'nsumdrops';'ndrops';...
    'Dmax';'Dmax_i';'bin_nums';'D';'dD';'edges';'vfall';'mass';'Nd';'nd';...
    'odm_length';'odm_width';'t_step';'nbins';'ntimes';'Nw';'D0';'mu'};

% change rain to drops, since that's what we are using and saving 
DSDN.nsumdrops = DSDN.nsumrain;
DSDN.ndrops = DSDN.nrain;

% good / bad data flag
DSDN.bad_data = zeros(DSDN.ntimes,1);
DSDN.good_data = zeros(DSDN.ntimes,1);

if strcmp(cruise,'TGT') == 1
    DSDN.bad_data(DSDN.t > datenum(2018,5,15,0,0,0)) = 1;
end

wh_bad1 = find(DSDN.voltage < 0.12);
wh_bad2 = find(DSDN.Dmax > 15);
wh_bad = [wh_bad1; wh_bad2];
if ~isempty(wh_bad)
    DSDN.bad_data(wh_bad) = 1;
end

wh_good = find(DSDN.bad_data == 0);
DSDN.good_data(wh_good) = 1;


% add more human readable time arrays
% DSDN.tsince1950 = DSDN.t - datenum(1950, 1, 0, 0, 0, 0);
[DSDN.year, DSDN.month, DSDN.day, DSDN.hour, DSDN.minute] = datevec(DSDN.t);
DSDN.ndrops = DSDN.nrain;
DSDN.nsumdrops = DSDN.nsumrain;

% save this final version of the .mat file as well
for i = 1:length(vars_to_save)
    dsd.(vars_to_save{i}) = DSDN.(vars_to_save{i});
end
save(matfilename,'dsd');

%==========================================================================
% prep the data
%==========================================================================
% ff = fieldnames(DSDN);
ff = vars_to_save;

%%% if present, get rid of single-value or readme like metadata fields.
%%% These should be saved as attributes not variables.
ff(strcmp(ff, 'README') == 1) = [];
ff(strcmp(ff, 'nbins') == 1) = [];
ff(strcmp(ff, 'ntimes') == 1) = [];
ff(strcmp(ff, 'odm_width') == 1) = [];
ff(strcmp(ff, 'odm_length') == 1) = [];
ff(strcmp(ff, 't_step') == 1) = [];

nf = length(ff);

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

d2 = find(dr == DSDN.ntimes & dc == DSDN.nbins);
dt = find(dr == DSDN.ntimes & dc == 1);
db = find(dr == 1 & dc == DSDN.nbins);
dd = find(dr == 1 & dc == DSDN.nbins+1);
d1 = find(dr == 1 & dc == 1);
dread = find(dr > 5 & dc == 1 & strcmp(vartypes, 'cell') == 1);
nd = length(d2) + length(d1) + length(dt) + length(db) + length(dread) + length(dd);
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(dt); ff(db); ff(dd); ff(d2)]; 
% nf = length(fff);
vartypesf = [vartypes(dt); vartypes(db); vartypes(dd); vartypes(d2)];

ddr = nan(nf,1);
ddc = nan(nf,1);
for i = 1:nf
    eval(['thisvar = DSDN.' fff{i} ';']);
    [ddr(i), ddc(i)] = size(thisvar);
end
dd2 = find(ddr == DSDN.ntimes & ddc == DSDN.nbins);
ddt = find(ddr == DSDN.ntimes & ddc == 1);
ddb = find(ddr == 1 & ddc == DSDN.nbins);
ddd = find(ddr == 1 & ddc == DSDN.nbins+1);
nf2 = length(dd2);
nft = length(ddt);
nfb = length(ddb);
nfd = length(ddd);

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

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

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

% define the dimensions
Ltime=DSDN.ntimes;
Lbin=DSDN.nbins;
Lbind=DSDN.nbins+1;

% create ID for (define) the dimensions
timedimID=netcdf.defDim(nc,'time',Ltime);
bindimID=netcdf.defDim(nc,'D_bin_centers',Lbin);
binddimID=netcdf.defDim(nc,'D_bin_edges',Lbind);

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)),timedimID);
end

%%% bin variables
for i = 1:nfb
    j = ddb(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),bindimID);
end

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

%%% time x bin 2D variables
for i = 1:nf2
    j = dd2(i);
    theID(j) = netcdf.defVar(nc,char(fff(j)),char(vartypes(j)),[timedimID bindimID]);
end

% variable names and units
[uf, lnf] = ODM_assign_long_info_DSD(fff, nf);
for j = 1:nf
    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)));
end

% define some global attributes
varid = netcdf.getConstant('GLOBAL');
netcdf.putAtt(nc,varid,'creation_date', datestr(now));
netcdf.putAtt(nc,varid,'cruise', cruise_str);
netcdf.putAtt(nc,varid,'version_number', version_num);
netcdf.putAtt(nc,varid,'creator','Elizabeth J. Thompson');
netcdf.putAtt(nc,varid,'contact_email_1','elizabeth.thompson@noaa.gov');
netcdf.putAtt(nc,varid,'contact_email_2','drushka@apl.uw.edu');
netcdf.putAtt(nc,varid,'please_acknowledge:','Elizabeth J. Thompson NOAA PSL and Kyla Drushka APL-UW and Klepp (2015, 2018) papers listed below');
netcdf.putAtt(nc,varid,'acknowledge_institution_1:','Applied Physics Laboratory at the University of Washington (APL-UW)');
netcdf.putAtt(nc,varid,'acknowledge_institution_2:','NOAA Physics Sciences Laboratory (PSL)');
netcdf.putAtt(nc,varid,'odm_width','0.022 m = 22 mm = width of odm sampling volume between optical sensors. Used to compute cross sectional area for n(d) and N(d).');
netcdf.putAtt(nc,varid,'odm_length','0.120 m = 120 mm = length of odm sampling volume between optical sensors. Used to compute cross sectional area for n(d) and N(d).');
netcdf.putAtt(nc,varid,'odm_time_step','60 seconds = time interval of instrument. Drops are counted and logged for this time interval. Used to compute n(d) and N(d).');
netcdf.putAtt(nc,varid,'odm_location_on_ship',odm_location_on_ship);
netcdf.putAtt(nc,varid,'bin_arrangemnts','The diameters provided in D are the bin centerpoints (dimension = nbins). The edges are provided as well (dimension = nbins+1).');
netcdf.putAtt(nc,varid,'reference1','Klepp, C., 2015: The oceanic shipboard precipitation measurement network for surface validation- OceanRAIN. Atmos. Res., 163, 74?90. http://dx.doi.org/10.1016/j.atmosres.2014.12.014');
netcdf.putAtt(nc,varid,'reference2','Klepp, C., 2018: Data Descriptor: OceanRAIN, a new in-situ shipboard global ocean surface-reference dataset of all water cycle components. Nature Scientific Data. DOI: 10.1038/sdata.2018.122');
netcdf.putAtt(nc,varid,'instrument','ODM-470 is an optical disdrometer developed by Christian Klepp (Klepp 2015, 2018)');
netcdf.putAtt(nc,varid,'about_instrument','ODM-470 was specifically designed for ships in that it rotates freely with the wind using a vane. A cup anemometer attached to the disdrometer measures the relative wind speed used to calculate the drop density, rain rate, etc... (see Klepp 2015).)');
% netcdf.putAtt(nc,varid,'data_quality_note','due to noise or imperfections, drops are often detected by unit even when it is definitely not raining. These do not produce a significant rain rate. Suggest thresholding data by a minimum number of total drops per time step, i.e. with the N variable');
netcdf.putAtt(nc,varid,'data_quality_note_1','the first 14 size bins, below 0.44 mm, do not contain data because these drops are too small to produce voltage drops greater than the noise floor');
netcdf.putAtt(nc,varid,'data_quality_note_2','128 size bins total, which are logarithmically spaced');
netcdf.putAtt(nc,varid,'data_quality_note_3','additional precipitation quantities were calculated with equations from Klepp 2015');
netcdf.putAtt(nc,varid,'data_quality_note_4','data are flagged as bad when voltage < 0.12 when sensor was cleaned, and when Dmax > 15 mm because that is unphysical for rain.');
netcdf.putAtt(nc,varid,'data_quality_note_5','Compared to level 1 data file, this level 2 file has fewer data points (times) because they are only included/available when Nw, D0, mu retrievals were possible. The index of the original time series when a retrieval was possible is in the variable good_gamma.');


% % 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 = DSDN.' fff{i} ';']);
%    eval(['dataID = ' fff{i} 'ID;']);
   eval(['data = DSDN.' 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);

end