% Sara Crepinsek 
% August 12th, 2014

%% Calling up program via script
%**************************************************************************
InputString = '20100921'; %InputString will be defined in Shupe script
%**************************************************************************
year = str2double(InputString(1:4)); 
month = str2double(InputString(5:6)); 
day = str2double(InputString(7:8)); 

[yday] = md2yd(year, month, day); %calculates Julian day given year, month, and day
doy = yday;
doy_string = num2str(doy);

%% Create file and calculate variables

if month == 1    Process_Month = 'Jan';
elseif month == 2     Process_Month = 'Feb';
elseif month == 3     Process_Month = 'Mar';
elseif month == 4     Process_Month = 'Apr';
elseif month == 5     Process_Month = 'May';
elseif month == 6     Process_Month = 'Jun';
elseif month == 7     Process_Month = 'Jul';
elseif month == 8     Process_Month = 'Aug';
elseif month == 9     Process_Month = 'Sep';
elseif month == 10    Process_Month = 'Oct';
elseif month == 11    Process_Month = 'Nov';
elseif month == 12    Process_Month = 'Dec';
end

Month = num2str(month,'%02d'); %adds leading zero
Process_Year = num2str(year,'%04d'); %adds leading zero);
Day = num2str(day,'%02d'); %adds leading zero);
Year_digit = year-2000; Process_Year_digit = num2str(Year_digit);
ly=@(yr)~rem(yr,400)|rem(yr,100)&~rem(yr,4);
ly=ly(year);

%Call up file to append/write to
cd (['W:\data\Summit\sodar\nc\'])
%cd (['\psd3data\arctic\summit\sodar\raw\'])
InputFile = ([Process_Year_digit,doy_string,'_SodarPC.nc.gz']);
gunzip(InputFile)
InputFileName = ([Process_Year_digit,doy_string,'_SodarPC.nc']);
%ncidIn = netcdf.open(InputFileName,'NC_WRITE');
inputfile = rd_netcdf(InputFileName);
inputatts = attribute_list(InputFileName);

%Convert file name convention
cd (['W:\data\Summit\sodar\nc\'])
%cd (['\psd3data\arctic\summit\sodar\raw\'])
OutputFileName = ['smtsodarX1.a1.',Process_Year,Month,Day,'.000000.cdf'];
ncid = netcdf.create(OutputFileName,'NOCLOBBER');
netcdf.close(ncid);
%copyfile(InputFileName,OutputFileName) %copy old file information/data to new file
%vinfo = ncinfo(OutputFileName);  vinfo = ncinfo(InputFileName);

%Create time arrays
base_time = inputfile.Time;
%base_time = ncread(OutputFileName,'Time'); %this is the time since 1970-01-01 00:00:00 0:00
base_time1 = base_time(1);
time_offset=base_time-base_time(1); %time since 1970 - time of start of file = time since start of file; %time_offset = seconds since base_time

% time to convert
BaseTime = base_time1; %UTC
epoch = datenum(1970,1,1,0,0,0); % Get the serial date number representing the epoch
ConvertDays = BaseTime/(86400); %60*60*24=86400, Divide the timestamp you want to convert by the number of seconds in a day
StringTime = epoch +ConvertDays; % Add this value to the serial date number representing the epoch
StringTime = double(StringTime);
formatOut = 'yyyy-mm-dd HH:MM:SS';
stringOutput = datestr(StringTime,formatOut);
string = [stringOutput,' 0:00'];

%create height array
sodarheight1 = 0.00:1.0:511.0; %float vector range
sodarheight = sodarheight1*0.3125+0.3125; %convert float vector heights to meters


%% Define global attributes
ncid = netcdf.open([OutputFileName],'NC_WRITE');%Open the file

ncwriteatt(OutputFileName,'/','date_created',string);
ncwriteatt(OutputFileName,'/','product_name','smtsodarX1.a1');
ncwriteatt(OutputFileName,'/','proc_level','a1');
ncwriteatt(OutputFileName,'/','number_input_data_streams','1');
ncwriteatt(OutputFileName,'/','input_data_stream_01','daily raw sodar files (i.e., YYJJJ_SodarPC.nc)');
ncwriteatt(OutputFileName,'/','description','NOAA sodar raw measurements data set converted from raw netCDF files');
ncwriteatt(OutputFileName,'/','PI_names','Matthew Shupe');
ncwriteatt(OutputFileName,'/','contact','matthew.shupe@noaa.gov');
ncwriteatt(OutputFileName,'/','project_name','ICECAPS: Integrated Characterization of Energy, Clouds, Atmospheric state, and Precipitation at Summit');
ncwriteatt(OutputFileName,'/','references','Shupe,M.D., D.D.Turner, V.P.Walden, R.Bennartz, and Coauthors, 2013:High and Dry:New observations of tropospheric and cloud properties above the Greenland ice Sheet. Bulletin of the American Meteorological Society, 94, 169-186.');
ncwriteatt(OutputFileName,'/','site_id','smt');
ncwriteatt(OutputFileName,'/','facility_id','Summit Station, Greenland');
ncwriteatt(OutputFileName,'/','instrument','NOAA sodar');
ncwriteatt(OutputFileName,'/','support','The NOAA Earth System Research Laboratory has provided this instrument for the NSF-funded ICECAPS project.');


%% Add/define Variables

%open the new file to write to
%ncid = netcdf.open([OutputFileName],'NC_WRITE');%Open the file

%base_time
netcdf.reDef(ncid)%put netcdf file in DEFINE mode
dimidB = netcdf.defDim(ncid,'base_time',1);%Define the dimensions
BaseTime_ID = netcdf.defVar(ncid,'base_time','double',[dimidB]);%Define the main variable ()
varidAttB = netcdf.inqVarID(ncid,'base_time');
netcdf.putAtt(ncid,varidAttB,'string',string);%Define IDs for the dimension variables
netcdf.putAtt(ncid,varidAttB,'long_name','Base time in Epoch');%Define IDs for the dimension variables
netcdf.putAtt(ncid,varidAttB,'units','seconds since 1970-1-1 0:00:00 0:00');%Define IDs for the dimension variables
netcdf.endDef(ncid);%We are done defining the NetCdf
netcdf.putVar(ncid,BaseTime_ID,base_time1);%Then store my main variable
    
%time_offset
netcdf.reDef(ncid)%put netcdf file in DEFINE mode
DimensionT = length(time_offset);
dimidT = netcdf.defDim(ncid,'time_offset',DimensionT);
TimeOffset_ID = netcdf.defVar(ncid,'time_offset','double',[dimidT]);%Define the main variable ()
varidAttT = netcdf.inqVarID(ncid,'time_offset');
netcdf.putAtt(ncid,varidAttT,'long_name','Time offset from base_time');%Define IDs for the dimension variables
netcdf.putAtt(ncid,varidAttT,'units',['seconds since ',string]);%Define IDs for the dimension variables
netcdf.endDef(ncid);%We are done defining the NetCdf
netcdf.putVar(ncid,TimeOffset_ID,time_offset);%Then store my main variable

% Time
netcdf.reDef(ncid)%put netcdf file in DEFINE mode
cdfPut(inputfile.Time,'Time','NC_INT',inputatts.Time,length(inputfile.Time),ncid);
netcdf.endDef(ncid);%We are done defining the NetCdf

%Height
netcdf.reDef(ncid)%put netcdf file in DEFINE mode
DimensionH = length(sodarheight);
dimidH = netcdf.defDim(ncid,'Height',DimensionH);%Define the dimensions
Height_ID = netcdf.defVar(ncid,'Height','double',[dimidH]);%Define the main variable ()
varidAttH = netcdf.inqVarID(ncid,'Height');
netcdf.putAtt(ncid,varidAttH,'long_name','height');%Define IDs for the dimension variables
netcdf.putAtt(ncid,varidAttH,'units','m, AGL');%Define IDs for the dimension variables
netcdf.endDef(ncid);%We are done defining the NetCdf
netcdf.putVar(ncid,Height_ID,sodarheight);%Then store my main variable

% Latitude
Atts.long_name = 'latitude';
Atts.units = 'degrees North';
Atts.value = '72.6';
netcdf.reDef(ncid)%put netcdf file in DEFINE mode
cdfPut(inputfile.Latitude,'Latitude','NC_INT',Atts,length(inputfile.Latitude),ncid);

% Lonitude
Atts.long_name = 'longitude';
Atts.units = 'degrees East';
Atts.value = '-38.42';
cdfPut(inputfile.Longitude,'Longitude','NC_INT',Atts,length(inputfile.Longitude),ncid);

% Elevation
Atts.long_name = 'elevation';
Atts.units = 'Meters, MSL';
Atts.value = '3255';
cdfPut(inputfile.Elevation,'Elevation','NC_INT',Atts,length(inputfile.Elevation),ncid);

% scfPut = (datain,'varname','prec',ATTs,size)
% if two dims: [size(x,1) size(x,2)] for size

% Azimuth_Beam1
cdfPut(inputfile.Azimuth_Beam1,'Azimuth_Beam1','NC_INT',inputatts.Azimuth_Beam1,length(inputfile.Azimuth_Beam1),ncid);

% Elevation_Beam1
cdfPut(inputfile.Elevation_Beam1,'Elevation_Beam1','NC_INT',inputatts.Elevation_Beam1,length(inputfile.Elevation_Beam1),ncid);

% TxFrequency
cdfPut(inputfile.TxFrequency,'TxFrequency','NC_INT',inputatts.TxFrequency,length(inputfile.TxFrequency),ncid);

% PulseLength
cdfPut(inputfile.PulseLength,'PulseLength','NC_INT',inputatts.PulseLength,length(inputfile.PulseLength),ncid);

% SampleRate
cdfPut(inputfile.SampleRate,'SampleRate','NC_INT',inputatts.SampleRate,length(inputfile.SampleRate),ncid);

% RecordDelay
cdfPut(inputfile.RecordDelay,'RecordDelay','NC_INT',inputatts.RecordDelay,length(inputfile.RecordDelay),ncid);

% SampleLength
cdfPut(inputfile.SampleLength,'SampleLength','NC_INT',inputatts.SampleLength,length(inputfile.SampleLength),ncid);

% AirTemperatureC
cdfPut(inputfile.AirTemperatureC,'AirTemperatureC','NC_INT',inputatts.AirTemperatureC,length(inputfile.AirTemperatureC),ncid);

% HeaderIndex
cdfPut(inputfile.HeaderIndex,'HeaderIndex','NC_INT',inputatts.HeaderIndex,length(inputfile.HeaderIndex),ncid);

% QcFlag
%atts replaces []
%atts.long_name = '';
%att.units = '';
cdfPut(inputfile.QcFlag,'QcFlag','NC_INT',[],length(inputfile.QcFlag),ncid);

% ReflectivityData
cdfPut(inputfile.ReflectivityData,'ReflectivityData','float',inputatts.ReflectivityData,[size(inputfile.ReflectivityData,1) size(inputfile.ReflectivityData,2)],ncid);


%We're done, close the netcdf
netcdf.close(ncid)

%% Plotting

vinfo_new = ncinfo(OutputFileName);
Refdata = ncread(OutputFileName,'ReflectivityData');
Ref2 = transpose(Refdata);
X=time_offset./3600; %puts plot in 24hr per day
Y=sodarheight;
Z=transpose(Ref2);

colormap(jet(256)); %code-line from Chris Cox
imagesc(X,Y,log(Z)); %Code-line from Chris Cox
set(gca,'ydir','normal'); %flips y-axis because imagesc plots it backwards
caxis([1 7]); %rough estimate of the limits of log(ref)
t = colorbar;
set(get(t,'xlabel'),'string','log(dRef)'); %maybe this is 'ylabel'
ylabel('Height [m]')
hold on
txtX=(['Time [hours, UTC]']);
xlabel(txtX); hold on
% Create textbox
annotation(gcf,'textbox',...
    [0.00775675675675676 0.0053475935828877 0.140047297297297 0.0374331550802139],...
    'String',{'NOAA/ESRL/PSD/Weather&Climate Physics'},...
    'FontSize',6,...
    'FitBoxToText','on','LineStyle','none');
site = 'Summit Sodar Reflectivity - ';
txt2={strcat([site, Day, Process_Month, Process_Year])};
title(txt2)      
grid on

cd (['W:\data\Summit\sodar\nc\'])
%cd (['\psd3data\arctic\summit\sodar\quicklooks\'])
FileName = strcat(Process_Year,Month,Day);
saveas(gcf,['smtsodarX1.a1.',FileName,'.000000.jpg'])

close all
