clear all;
close all;
set_machine_et;

%%%     --- code 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

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

cruise_str = [cruise '_' yrstr];

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

version_num = 'v3';

fillval = -9999;

% load latest DSD data file (might be updating if on the cruise currently)
these_files = dir([thepath cruise '_dsd_level_1_all_data_' version_num '*.mat']);
ncfilename = [thepath cruise '_' yrstr '_dsd_level_1_' version_num '.nc'];
matfilename = [thepath cruise '_' yrstr '_dsd_level_1_' version_num '.mat'];

files = char(these_files(:).name);
[nfiles, nchars] = size(files);
load([thepath files(end,:)]);
disp(['latest DSD file: ' files(end,:)]);

% archive formatted name
revision_num = 'R0';
day1_str = datestr(DSD.t(1),'yyyymmdd');
dayn_str = datestr(DSD.t(end),'yyyymmdd');
ncfilename_new =[thepath cruise '-disdrometer_RV-' ship '_' day1_str '_' revision_num '_thru_' dayn_str '.nc'];


if strcmp(cruise,'SPURS2') == 1
    kfile_2 = '/Users/eliz/DATA/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_mod.mat';
    if isfile(kfile_2) == 0
        k_s16 = load('/Users/eliz/DATA/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_2016.mat');
        k_s17 = load('/Users/eliz/DATA/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];
        k.lat = [k_s16.lat ; k_s17.lat];
        k.lon = [k_s16.lon ; k_s17.lon];
    else
        load(kfile_2)    
    end    
    DSD.Pfortran = k.p;
    DSD.lat = k.lat;
    DSD.lon = k.lon;
else
    DSD.Pfortran = ones(length(DSD.t),1)*fillval;
    DSD.lat = ones(length(DSD.t),1)*fillval;
    DSD.lon = ones(length(DSD.t),1)*fillval;
end


DSD.time = (DSD.t - datenum(theyear,1,1,0,0,0))*24*3600; 

%%% removed 'qc_flag'

vars_to_save = {'time';'lat';'lon';...
    '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'};

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

flag_used = 0; % only for TGT
if flag_used == 1
    % good / bad data flag
    DSD.bad_data = zeros(DSD.ntimes,1);
    DSD.qc_flag = zeros(DSD.ntimes,1);  % e.g. good data

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

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

    wh_good = find(DSD.bad_data == 0);
    DSD.qc_flag(wh_good) = 1; % e.g. good data

end

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

% save this final version of the .mat file as well, reset any nans to
% fillval
for i = 1:length(vars_to_save)
    thevar = DSD.(vars_to_save{i});
    thevar(find(isnan(thevar) == 1)) = fillval; 
    dsd.(vars_to_save{i}) = thevar;
end
save(matfilename,'dsd');


%==========================================================================
% prep the data
%==========================================================================
% 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.
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);
vminf = nan(nf,1);
vmaxf = nan(nf,1);
for i = 1:nf
    eval(['thisvar = DSD.' 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 == DSD.ntimes & dc == DSD.nbins);
dt = find(dr == DSD.ntimes & dc == 1);
db = find(dr == 1 & dc == DSD.nbins);
dd = find(dr == 1 & dc == DSD.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 = DSD.' fff{i} ';']);
    [ddr(i), ddc(i)] = size(thisvar);
    vminf(i) = min(min(thisvar));
    vmaxf(i) = max(max(thisvar));
end
dd2 = find(ddr == DSD.ntimes & ddc == DSD.nbins);
ddt = find(ddr == DSD.ntimes & ddc == 1);
ddb = find(ddr == 1 & ddc == DSD.nbins);
ddd = find(ddr == 1 & ddc == DSD.nbins+1);
nf2 = length(dd2);
nft = length(ddt);
nfb = length(ddb);
nfd = length(ddd);

% redo weird reflectivity value so it's not negative infinity
vminf(vminf < -1e20) = -40;

%==========================================================================
% 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, 'CLASSIC_MODEL');

% define the dimensions
Ltime=DSD.ntimes;
Lbin=DSD.nbins;
Lbind=DSD.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


%%% special for time
netcdf.putAtt(nc,0,'long_name','time');
netcdf.putAtt(nc,0,'axis','T');
netcdf.putAtt(nc,0,'standard_name','time');
netcdf.putAtt(nc,0,'units','seconds since 2016-01-01T00:00:00Z');
netcdf.putAtt(nc,0,'coverage_content_type','coordinate');
netcdf.putAtt(nc,0,'comments','to convert to matlab datetime: t = datenum(year, 1, 1, 0, 0, time)');


% %%% special for the qc_flag
% netcdf.putAtt(nc,2,'long_name','quality control flag');
% netcdf.putAtt(nc,2,'flag_values','[0 1]');
% netcdf.putAtt(nc,2,'flag_meanings','0 = bad data, 1 = good data. 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,2,'coordinates','time latitude longitude');
% netcdf.putAtt(nc,2,'_FillValue', NaN	);
% netcdf.putAtt(nc,2,'axis','T');
% netcdf.putAtt(nc,2,'coverage_content_type','auxiliaryInformation');


%%% special for the lat
netcdf.putAtt(nc,1,'long_name','latitude of measurement');
netcdf.putAtt(nc,1,'axis','Y');
netcdf.putAtt(nc,1,'standard_name','latitude');
netcdf.putAtt(nc,1,'units','degrees_north');
netcdf.putAtt(nc,1,'coverage_content_type','coordinate');
netcdf.putAtt(nc,1,'valid_min','-90');
netcdf.putAtt(nc,1,'valid_max','90');




%%% special for the lon
netcdf.putAtt(nc,2,'long_name','longitude of measurement');
netcdf.putAtt(nc,2,'axis','X');
netcdf.putAtt(nc,2,'standard_name','longitude');
netcdf.putAtt(nc,2,'units','degrees_east');
netcdf.putAtt(nc,2,'coverage_content_type','coordinate');
netcdf.putAtt(nc,2,'valid_min','-180');
netcdf.putAtt(nc,2,'valid_max','180');



% variable names and units
[snf, lnf, uf, cf, ccf, crdf] = ODM_assign_long_info_DSD(fff, nf, yrstr);
% for j = 1:nf
%     if isempty(snf{j}) == 1
%         disp(['missing standard name 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
%     if isempty(uf{j}) == 1
%         disp(['missing unit for field ' sprintf('%i',j) ' = ' string(fff{j})]);
%     end
% end

% define the units and long names
for i = 1:nf
    if i >= 4
        disp(fff(i));
        if isempty(snf{i})==0 
            netcdf.putAtt(nc,theID(i),'standard_name',char(snf{i}));
        end
        netcdf.putAtt(nc,theID(i),'long_name',char(lnf(i)));
        netcdf.putAtt(nc,theID(i),'units',char(uf(i)));
        netcdf.putAtt(nc,theID(i),'coverage_content_type',char(ccf(i)));
        if isempty(crdf{i})==0 
            netcdf.putAtt(nc,theID(i),'coordinates',char(crdf(i)));
        end
        netcdf.putAtt(nc,theID(i),'valid_min',vminf(i));
        netcdf.putAtt(nc,theID(i),'valid_max',vmaxf(i));
        if isempty(cf{i})==0 
            netcdf.putAtt(nc,theID(i),'comment',char(cf(i)));
        end
        netcdf.defVarFill(nc,theID(i),false,fillval);
    end
end



% define some global attributes
varid = netcdf.getConstant('GLOBAL');
netcdf.putAtt(nc,varid,'DOI','10.5067/SPUR2-DISDR');
netcdf.putAtt(nc,varid,'title','SPURS2 Disdrometer Data');
netcdf.putAtt(nc,varid,'file_version',1.0);
netcdf.putAtt(nc,varid,'Conventions','CF-1.8, ACDD-1.3');
netcdf.putAtt(nc,varid,'keywords_vocabulary','NASA Global Change Master Directory (GCMD) Science Keywords');
netcdf.putAtt(nc,varid,'keywords','Earth Science > Atmosphere > Precipitation > Droplet Size; Earth Science > Atmosphere > Clouds > Cloud Microphysics > Particle Size Distribution; Earth Science > Atmosphere > Clouds > Cloud Microphysics > Cloud Droplet Concentration/Size; Earth Science > Atmosphere > Precipitation > Precipitation Rate'); 

netcdf.putAtt(nc,varid,'summary','ODM-470 optical disdrometer data collected during SPURS-2 field campaign in the Eastern tropical Pacific between 2016-2017 in two different research cruises separated by almost a year');
netcdf.putAtt(nc,varid,'uuid','c35e5960-e3fa-4124-b157-d2a7fc7d6495');
netcdf.putAtt(nc,varid,'id','SPURS2-DISDR');
netcdf.putAtt(nc,varid,'naming_authority','nasa.gov');
netcdf.putAtt(nc,varid,'feature_type','trajectory');
netcdf.putAtt(nc,varid,'cdm_data_type','Trajectory');
netcdf.putAtt(nc,varid,'source','shipboard distrometer data');
netcdf.putAtt(nc,varid,'platform','Water-based Platforms>Vessels>Surface>R/V Roger Revelle');
netcdf.putAtt(nc,varid,'platform_vocabulary','GCMD platform keywords');
netcdf.putAtt(nc,varid,'instrument','In Situ/Laboratory Instruments > Recorders/Loggers > DISDROMETERS');
netcdf.putAtt(nc,varid,'instrument_vocabulary','NASA Global Change Master Directory (GCMD) Instrument Keywords');
netcdf.putAtt(nc,varid,'processing_level','L2');
netcdf.putAtt(nc,varid,'standard_name_vocabulary','CF Standard Name Table v72');
netcdf.putAtt(nc,varid,'references','Klepp et al. 2015 and 2018: http://dx.doi.org/10.1016/j.atmosres.2014.12.014 and https://doi.org/10.1038/sdata.2018.122');
netcdf.putAtt(nc,varid,'acknowledgement','Elizabeth J. Thompson NOAA PSL and Kyla Drushka APL-UW and Klepp (2015, 2018) papers listed in references global attribute');
netcdf.putAtt(nc,varid,'comment',['The location on the ship each year was: ' odm_location_on_ship '. ODM-470 is an optical disdrometer developed by Christian Klepp (Klepp 2015, 2018) but is no longer supported. It was specifically designed for ships but still has known sampling issues on this platform. The device 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. 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. additional precipitation quantities were calculated with equations from Klepp 2015. The ODM cannot detect drops < 0.44 diameter and this significantly impairs the ability of this data to accurately estimate DSD parameters such as Nw, D0, Dm with precision. Therefore, these parameters are not provided here and ODM data are not recommended to be used for this tropical warm rain regime. Tests were conducted in PISTON and SPURS-2 compared to 2DVD results on Gan and Manus Islands. ODM data appear most problematic in distinguishing shallow weak convection from stratiform rain becauese there difference is normally distinguishable by a high concentration of small drops in the former, but this cannot be measured well be ODM. The mean D0 for tropical rain in Manus Is was 1.1 mm (Thompson et al. 2015).']);
netcdf.putAtt(nc,varid,'license','creative commons ccby4.0');
netcdf.putAtt(nc,varid,'product_version','1.0');
netcdf.putAtt(nc,varid,'history','final processed data. Data have been removed (fill value applied) where voltage < 0.12 V because these signals cannot be distinguished from noise.');
netcdf.putAtt(nc,varid,'metadata_link','https:doi.org/10.5067/SPUR2-DISDR');
netcdf.putAtt(nc,varid,'creator_name','Elizabeth J. Thompson, Haonan Chen, Kyla Drushka');
netcdf.putAtt(nc,varid,'creator_email','elizabeth.thompson@noaa.gov');
netcdf.putAtt(nc,varid,'creator_url','https://noaa.gov');
netcdf.putAtt(nc,varid,'creator_type','person');
netcdf.putAtt(nc,varid,'creator_institution','GOVERNMENT AGENCIES-U.S. FEDERAL AGENCIES>DOC>NOAA>DOC/NOAA/OAR/ESRL/PSL, ACADEMIC > UWA/APL, ACADEMIC > COLOSTATE/CIRA');
netcdf.putAtt(nc,varid,'institution','GOVERNMENT AGENCIES-U.S. FEDERAL AGENCIES>DOC>NOAA>DOC/NOAA/OAR/ESRL/PSL, ACADEMIC > UWA/APL, ACADEMIC > COLOSTATE/CIRA');
netcdf.putAtt(nc,varid,'institution_vocabulary','NASA Global Change Master Directory (GCMD) Institution Keywords');
netcdf.putAtt(nc,varid,'contributor_name','Frederick Bingham');
netcdf.putAtt(nc,varid,'contributor_role','Project Data Manager');
netcdf.putAtt(nc,varid,'contributor_type','person');
netcdf.putAtt(nc,varid,'contributor_institution','University of North Carolina Wilmington');
netcdf.putAtt(nc,varid,'project','Salinity Processes in the Upper Ocean Regional Study (SPURS)');
netcdf.putAtt(nc,varid,'program','NASA Ocean Salinity Science Team (OSST)');
netcdf.putAtt(nc,varid,'publisher_name','Physical Oceanography Distributed Active Archive Center, Jet Propulsion Laboratory, NASA');
netcdf.putAtt(nc,varid,'publisher_email','podaac@podaac.jpl.nasa.gov');
netcdf.putAtt(nc,varid,'publisher_url','https://podaac.jpl.nasa.gov'); 
netcdf.putAtt(nc,varid,'publisher_type','institution');
netcdf.putAtt(nc,varid,'publisher_institution','NASA/JPL/PODAAC');
netcdf.putAtt(nc,varid,'sea_name','OCEAN > PACIFIC OCEAN > EASTERN PACIFIC OCEAN');
netcdf.putAtt(nc,varid,'date_created', datestr(now,'yyyy-mm-ddTHH:MM:SSZ'));
netcdf.putAtt(nc,varid,'geospatial_vertical_positive','down');
netcdf.putAtt(nc,varid,'geospatial_lat_min',min(DSD.lat));
netcdf.putAtt(nc,varid,'geospatial_lat_max',max(DSD.lat));
netcdf.putAtt(nc,varid,'geospatial_lat_units','degrees_north');
netcdf.putAtt(nc,varid,'geospatial_lon_min',min(DSD.lon));
netcdf.putAtt(nc,varid,'geospatial_lon_max',max(DSD.lon));
netcdf.putAtt(nc,varid,'geospatial_lon_units','degrees_east');
netcdf.putAtt(nc,varid,'time_coverage_start', datestr(min(DSD.t), 'yyyy-mm-ddTHH:MM:SSZ'));
netcdf.putAtt(nc,varid,'time_coverage_end',datestr(max(DSD.t), 'yyyy-mm-ddTHH:MM:SSZ'));


%%%%%%%%%%%%%%%


% % 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 = DSD.' 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);

% copy the modified file to the new name that matches archive formatting
command = ['cp ' ncfilename ' ' ncfilename_new];
system(command);

end