% single-wavelength backscatter lidars (ceilometers)

% The volume backscatter coefficient represents the portion of light that is 
% reflected back towards the ceilometer from a distance z (e.g., by water droplets).

clear all;
close all;
    
%%% get the current UTC time for file creation
%check is MST or MDT. If equal 1, mean we are in Daylight time (MDT=UTC-6, MST=UTC-7)
if isdst(datetime(datestr(now),'TimeZone','America/Denver'))==1  %means MDT
    Tstring=datestr(now+datenum(0000,1,0,6,0,0),'yyyy-mm-dd HH:MM Z');
else %means MST
    Tstring=datestr(now+datenum(0000,1,0,7,0,0),'yyyy-mm-dd HH:MM Z');
end    

%% get flux file for  lat and lon

ffile = '/Users/eliz/DATA/PISTON_2019/Sally_Ride/flux/Processed/final/PISTON-nav-met-sea-flux-10min_RV-Sally-Ride_20190906_R1_thru_20190925.nc';
ncload2(ffile);
% ncdisp(ffile);
finfo = ncinfo(ffile);
fnames = {finfo.Variables.Name};
% other options in finfo: Attributes, Dimensions, Variables, Format, Filename

% ncload2 loads all variables in a netcdf file by variable name.
% Since multiple files have similar variables like time...
% Load all the variables from the flux file into a structure f, i.e.  
% f.$variable_name, and then clear the variable name from the workspace
nf = length(fnames);
for i = 1:nf
    eval(['f.' fnames{i} ' = ' fnames{i} ';']);
    eval(['clear ' fnames{i} ';']);
end
f = orderfields(f);
f.t = datenum(2019,1,1,0,0,f.time);
f.jd = f.t - datenum(2019,0,0,0,0,0);

%% get ceilo data

path_out='/Users/eliz/DATA/PISTON_2019/Sally_Ride/Ceilometer/Processed/all/';
path_in=path_out;
vnum =  'v2';
shipname = 'RV-Sally-Ride';
cruisename = 'PISTON';
cruise = 'PISTON_2019';
str_yr = cruise(end-3:end);               % year string


%%%%%%% 10-min or 1-hour file:

%                                 jd_10min: [4896×1 double]
%                                  t_10min: [4896×1 double]
%                                 te_10min: [4896×1 double]
%                               frac_cloud: [4896×1 double]
%                               frac_clear: [4896×1 double]
%                        count_cloud_total: [4896×1 double]
%                count_cloud_layers_1_to_3: [4896×1 double]
%                     count_cloud_obscured: [4896×1 double]
%     count_cloud_part_obsured_transparent: [4896×1 double]
%                              count_clear: [4896×1 double]
%                               cloud_base: [4896×1 double]
%                              prcnt_cld_5: [4896×1 double]
%                             prcnt_cld_10: [4896×1 double]
%                             prcnt_cld_16: [4896×1 double]
%                             prcnt_cld_50: [4896×1 double]
%                             prcnt_cld_84: [4896×1 double]

%%%%%%% 15-sec file

%                    jd: [195840×1 double]
%                     t: [195840×1 double]
%                    te: [195840×1 double]
%            cloud_code: [195840×1 double]
%      cloud_base_first: [195840×1 double]
%     cloud_base_second: [195840×1 double]
%      cloud_base_third: [195840×1 double]
%                ceil_b: [770×195840 double]


for k = 1:3
clear c;

% rename some variables for consistency and make units km for all vars
if k == 1  %%% 15-sec
    load([path_in cruise '_ceilometer_15sec.mat']);
    c = cc;
    timetype = '15sec';
    c.height = c.height/1000; % convert to km
    c.cloud_base_first = c.cloud_base_first/1000; % convert to km
    c.cloud_base_second = c.cloud_base_second/1000; % convert to km
    c.cloud_base_third = c.cloud_base_third/1000; % convert to km
    maxheight = max(c.height); % save this
elseif k == 2 %%% 10-min
    load([path_in cruise '_ceilometer_10min.mat']);
    c = c10;
    c.t = c.t_10min;
    c.jd = c.jd_10min;
    c.te = c.te_10min;
    timetype = '10min';
elseif k == 3 %%% 1-hour
    load([path_in cruise '_ceilometer_1hour.mat']);
    timetype = '1hour';
    %%% take _hr off of all the field names
    chrf = fields(chr);
    for j = 1:length(chrf)
        varname = char(chrf{j});
        newname = varname(1:end-3);
        c.(newname) = chr.(chrf{j});
    end
end


%%% take out eez time periods so that it matches the flux file time base
eez_in     = datenum(2019,9, 6, 5,0,0);
eez_out    = datenum(2019,9,25,12,0,0);

wh_eez = find(c.t < eez_in | c.t > eez_out);

%%% remove eez
fc = fields(c);

ntime = length(c.t);
for j = 1:length(fc)
    x = c.(fc{j});
    [nrow,ncol] = size(x);
    if ncol == 1 && nrow == ntime % for time series, but not height
        x(wh_eez) = [];
    elseif nrow < ncol % for vertical profiles in 15-sec data
        x(:,wh_eez) = [];
    end
    c.(fc{j}) = x; 
end

jdy=c.jd;

lat=interp1(f.jd, f.lat, c.jd);
lon=interp1(f.jd, f.lon, c.jd);

startday = datestr(c.t(1),'YYYYmmDD');
endday = datestr(c.t(end),'YYYYmmDD');

file_out=[cruisename '-ceilometer-' timetype '_' shipname '_' startday '_thru_' endday '_' vnum '.nc']
cmode = netcdf.getConstant('NETCDF4');
cmode = bitor(cmode,netcdf.getConstant('CLASSIC_MODEL'));
ncid = netcdf.create(strcat(path_out,file_out),cmode);


%% Define the dimensions (note in NetCDF, first dimension is 0)

if k == 1
    [z_dimlen,time_dimlen]=size(c.ceil_b);
    z_dimID= netcdf.defDim(ncid,'z',z_dimlen); 
    time_dimID = netcdf.defDim(ncid,'time',time_dimlen);
else
    [time_dimlen,z_dimlen]=size(c.t);
    time_dimID = netcdf.defDim(ncid,'time',time_dimlen);
end
    

%% Define the coordinate and dimension variables (time,latitude,longitude,trajectory,...): IDs and attributes
% profile_ID=netcdf.defVar(ncid,'profile','int',prof_dimID);
%     netcdf.putAtt(ncid, profile_ID, 'long_name', 'Backscatter profiles from ceilometer');
%     netcdf.putAtt(ncid, profile_ID, 'cf_role', 'profile_id');  
time_ID=netcdf.defVar(ncid,'time','double',time_dimID);
%     netcdf.putAtt(ncid, time_ID, 'long_name', 'Time');
    netcdf.putAtt(ncid, time_ID, 'standard_name', 'time');    
%     netcdf.putAtt(ncid, time_ID, 'units', 'decimal day of year');
%     netcdf.putAtt(ncid, time_ID, 'units', sprintf('seconds since %s UTC', datestr(datenum(datestr('01-Jan-2020 00:00:00')),0)));
    netcdf.putAtt(ncid, time_ID, 'units', ['seconds since ' str_yr '-01-01 00:00:00 UTC'])
    netcdf.putAtt(ncid, time_ID, 'calendar', 'standard');
    netcdf.putAtt(ncid, time_ID, 'comment', 'Time marks the start point of each duration interval');
    netcdf.putAtt(ncid, time_ID , 'coverage_content_type', 'coordinate');
%     netcdf.putAtt(ncid, time_ID, 'cf_role', 'profile_id');
    netcdf.defVarFill(ncid, time_ID, false, -9999.);
latitude_ID=netcdf.defVar(ncid,'lat','double',time_dimID);
%     netcdf.putAtt(ncid, latitude_ID, 'long_name', 'latitude');
    netcdf.putAtt(ncid, latitude_ID, 'standard_name', 'latitude');    
    netcdf.putAtt(ncid, latitude_ID, 'units', 'degree_north');
    netcdf.putAtt(ncid, latitude_ID, 'actual_range', [min(lat),max(lat)]);
    netcdf.putAtt(ncid, latitude_ID, 'valid_range', [-90., 90.]);
    netcdf.putAtt(ncid, latitude_ID , 'coverage_content_type', 'coordinate');
    netcdf.putAtt(ncid, latitude_ID , 'instrument', 'Ship GPS');
    netcdf.defVarFill(ncid, latitude_ID, false, -9999.);
longitude_ID=netcdf.defVar(ncid,'lon','double',time_dimID);
%     netcdf.putAtt(ncid, longitude_ID, 'long_name', 'longitude');
    netcdf.putAtt(ncid, longitude_ID, 'standard_name', 'longitude');    
    netcdf.putAtt(ncid, longitude_ID, 'units', 'degree_east');
    netcdf.putAtt(ncid, longitude_ID, 'actual_range', [min(lon),max(lon)]);%)[min(lon), max(lon)]);
    netcdf.putAtt(ncid, longitude_ID, 'valid_range', [-180., 180.]);
    netcdf.putAtt(ncid, longitude_ID , 'coverage_content_type', 'coordinate');
    netcdf.putAtt(ncid, longitude_ID , 'instrument', 'Ship GPS');
    netcdf.defVarFill(ncid, longitude_ID, false, -9999.);
    
%% Define global_attributes of the file
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'title','PISTON (Propagation of Intra-Seasonal Tropical Oscillations) 2019 Field campaign, : ceilometer 10min statistics from the RV Sally Ride by the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Lab (PSL)');
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'title','PISTON (Propagation of Intra-Seasonal Tropical Oscillations) 2019 Field campaign, : ceilometer 10min statistics from the RV Sally Ride by the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Lab (PSL) in collaboration with Joe (Harindra) Fernando at University of Notre Dame');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'summary','The temporal statistics in this file were computed from the 15-sec data that are output automatically by the Vaisala software. This covers different cloud conditions categorized by a cloud code : 0 = clear sky, 1 = one cloud layer detected, 2 = two cloud layers detected, 3 = three cloud layers detected,4 = totally obscured, 5 = partly obscured but determined transparent. Different percentile values were tested to determine how to most accurately characterize cloud base height from the lowest layer for this experiment. For this cruise the 5th percentile seemed to work best for cloud base height. All percentiles calculated are saved and explained in the comment global attribute.') 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'keywords','Earth Science > Atmosphere > Clouds > Cloud Properties > Cloud Base Height, Earth Science > Atmosphere > Clouds > Cloud Properties > Cloud Fraction, Earth Remote Sensing Instruments > Passive Remote Sensing > Photon/Optical Detectors >	Ceilometers, In Situ Ocean-based Platforms > SHIPS ');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'keywords_vocabulary', 'GCMD Earth Science Keywords. Version 8.6') ;
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'program','Funding agencies: NOAA Climate Program Office Climate Variability and Predictability, ONR Marine Meteorology ');

% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'campaign','EUREC4A');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'platform','Research Vessel Sally Ride');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'project','PISTON 2019');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'sea_name','Philippine Sea');

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'institution',sprintf('NOAA Physical Sciences Laboratory, 325 Broadway, Boulder, CO 80305'));

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lat_min', min(lat)); %Describes a simple lower latitude limit. 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lat_max', max(lat)); %Describes a simple upper latitude limit. 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lon_min', min(lon)); %Describes a simple lower longitude limit. 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lon_max', max(lon)); %Describes a simple upper longitude limit. 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'time_coverage_start', datestr(c.t(1),0)); %Describes the time of the first data point in the data set.
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'time_coverage_end', datestr(c.t(end),0)); % Describes the time of the last data point in the data set.

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'publisher_type','person'); %Specifies type of creator with one of the following: 'person', 'group', 'institution', or 'position
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'publisher_institution','NOAA'); %The institution that presented the data file or equivalent product to users; should uniquely identify the institution
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'publisher_name','Elizabeth J. Thompson'); %The name of the person (or other entity specified by the publisher_type attribute) responsible for publishing the data file or product to users, with its current metadata and format
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'publisher_email','elizabeth.thompson@noaa.gov'); %The email address of the person (or other entity specified by the publisher_type attribute) responsible for publishing the data file or product to users, with its current metadata and format;
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'contributor_name','Chris Fairall, Elizabeth J. Thompson, Sergio Pezoa, Byron Blomquist, Simon de Szoeke, Ludovic Bariteau'); %The name of any individuals, projects, or institutions that contributed to the creation of this data;

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'instrument','Ceilometer: Vaisala CL31 single-wavelength backscatter lidar');

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'product_version','v1');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'history','v1: processed data');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'featureType','profile');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'cdm_data_type','profile'); 

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'date_created' ,Tstring); %Creation date of this version of the data(netCDF). Could be datestr(now,'yyyy-mm-ddTHH:MM:SS-0700')); for MDT
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'date_modified' ,Tstring); %The date on which the data was last modified. Note that this applies just to the data, not the metadata. Could be datestr(now,'yyyy-mm-ddTHH:MM:SS-0700')); for MDT

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'standard_name_vocabulary','CF Standard Name Table v41');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'ncei_template_version','NCEI_NetCDF_Trajectory_Template_v2.0');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'Conventions','CF-1.6, ACDD-1.3');

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'comment', sprintf('The ceilometer produced vertical profiles of backscatter intensity produced by refractive index gradients, such as differences between dry air above the trade wind inversion vs. the atmospheric boundary layer air that is laden with aerosol, sea spray, and water vapor. The principal use of the ceilometer is to track the strong backscatter returns associated with cloud base height. The base level was computed as the 5th percentile value of the heights detected from the first visible cloud level (cloud_code detections of first_cloud_base at 15-sec resolution) in a given hour or 10-min time period. The percentiles of first_cloud_base heights that were computed, tested, and saved in this file are: prcnt_cld_5, _10, _16, _50, _84. Comparing these heights can potentially yield some information on how sheared or how vertically-developed the non-raining lowest cloud layer was. However, these values are not all indicative of cloud base. The higher percentiles of first_cloud_base heights are mostly indicative of sheared lateral edges of clouds and detached tops of clouds whose bases have evaporated and dissipated. The signal is completely attenuated by rainfall, which was rare but did occur in this experiment. In the absence of rain, strong backscatter returns are produced by the cloud base, then the signal is rapidly attenuated above that level. The ceilometer data should never be used to determine cloud top height. The ceilometer will also detect a strong scattering return from the horizontally-sheared lateral edges of clouds overhead of the ship. For example, during situations when the cloud top becomes sheared laterally away from the location of the cloud base, the ceilometer return at that time will be higher than the actual cloud base. These times should be flagged and interpreted differently. Actual cloud base changes very slowly in these data. Detrained leftover upper-level cloud debris from mature dissipating clouds can also produce a detectable signal from the ceilometer within a km or so above the cloud base, such that the original cloud base has already evaporated. In these cases where the cloud base appears to change relatively rapidly, the lowest altitude signals from the ceilometer should be interpreted as cloud base. The ceilometer also detected middle-layer clouds in the atmosphere near 4-6 km that are separate from the low-level cumulus cloud field. A more thorough check of ceilometer statistics of cloud base and cloud fraction should be conducted between ceilometer, lidar, and LCL height calculated from mixed layer air parcels (not surface parcels).'))


if k > 1
    
    %%%%%% FOR 10-min and 1-hr files only

    %%% Define the main variables
    vari=c.cloud_base;
        cloud_base_ID=add_mainVar_attributes_v2020(ncid,'cloud_base_height','double',time_dimID,'cloud base height','cloud_base_altitude','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, cloud_base_ID , 'comment', 'Computed as the 5th percentile of first_cloud_base, cb_5_prct, the base height of first cloud layer detected, from 15 sec raw data over this time period.');
    vari=c.frac_cloud;
        frac_cloud_ID=add_mainVar_attributes_v2020(ncid,'frac_cloud','double',time_dimID,'cloud fraction','cloud_area_fraction_in_atmosphere_layer','1','physicalMeasurement',[0., 1.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, frac_cloud_ID , 'comment', 'cloud fraction of time period when cloud was detected in 15 sec data cloud_code, including codes 1-5 ranging from lower to upper level clouds and obscured clouds. The layer here ranges from surface to 7km');
    vari=c.frac_clear;
        frac_clear_ID=add_mainVar_attributes_v2020(ncid,'frac_clear','double',time_dimID,'clear sky fraction','','1','physicalMeasurement',[0., 1.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, frac_clear_ID , 'comment', 'fraction of time period when clear sky was detected in 15 sec data cloud_code, including code 0');
    vari=c.count_cloud_total;
        count_total_ID=add_mainVar_attributes_v2020(ncid,'count_total','int',time_dimID,'count of total cloud detections','','1','physicalMeasurement',[0., 100.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, count_total_ID , 'comment', 'count of any cloud classifications in this time period from 15 sec cloud_code data, either codes 1-5 ranging from lower to upper level clouds and obscured clouds');
    vari=c.count_cloud_layers_1_to_3;
        count_layers_1_to_3_ID=add_mainVar_attributes_v2020(ncid,'count_layers_1_to_3','int',time_dimID,'count of cloud detections in layers 1 to 3','','1','physicalMeasurement',[0., 100.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, count_layers_1_to_3_ID , 'comment', 'count of any cloud classifications from codes 1, 2, 3 for first, second, and third cloud layers detected in this time period from 15 sec cloud_code data');
    vari=c.count_cloud_obscured;
        count_cloud_obscured_ID=add_mainVar_attributes_v2020(ncid,'count_cloud_obscured','int',time_dimID,'count of obscured cloud detections','','1','physicalMeasurement',[0., 100.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, count_cloud_obscured_ID , 'comment', 'count of obscured cloud classifications in this time period from 15 sec cloud_code data, from code 4')
    vari=c.count_cloud_part_obsured_transparent;
        count_part_obsured_transparent_ID=add_mainVar_attributes_v2020(ncid,'count_part_obscured_transparent','int',time_dimID,'count of partly obscured cloud detections','','1','physicalMeasurement',[0., 100.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, count_part_obsured_transparent_ID , 'comment', 'count of partly obscured and transparent cloud classifications in this time period from 15 sec cloud_code data, from code 5')
    vari=c.count_clear;
        count_clear_ID=add_mainVar_attributes_v2020(ncid,'count_clear','int',time_dimID,'count of clear detections','','1','physicalMeasurement',[0., 100.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, count_clear_ID , 'comment', 'count of clear sky cloud classifications in this time period from 15 sec cloud_code data, from code 0')

    vari=c.prcnt_cld_5;
        prcnt_cld_5_ID=add_mainVar_attributes_v2020(ncid,'cb_5_prct','double',time_dimID,'cloud base height from 5th percentile','cloud_base_altitude','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, prcnt_cld_5_ID , 'comment', 'This is the same as cloud_base_height for this dataset since appeared to be the most consistent threshold. This is the 5th percentile of first_cloud_base, the base height of first cloud layer detected, from 15 sec raw data over this time period.');
    vari=c.prcnt_cld_10;
        prcnt_cld_10_ID=add_mainVar_attributes_v2020(ncid,'cb_10_prct','double',time_dimID,'cloud base height from 10th percentile','cloud_base_altitude','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, prcnt_cld_10_ID , 'comment', 'This is the 10th percentile of first_cloud_base, base height of first cloud layer detected, from 15 sec raw data over this time period.');
    vari=c.prcnt_cld_16;
        prcnt_cld_16_ID=add_mainVar_attributes_v2020(ncid,'cb_16_prct','double',time_dimID,'cloud base height from 16th percentile','cloud_base_altitude','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, prcnt_cld_16_ID , 'comment', 'This is the 16th percentile of first_cloud_base, base height of first cloud layer detected, from 15 sec raw data over this time period. Therefore this is one standard deviation below the median.');
    vari=c.prcnt_cld_50;
        prcnt_cld_50_ID=add_mainVar_attributes_v2020(ncid,'cb_50_prct','double',time_dimID,'cloud base height from 50th percentile','cloud_base_altitude','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, prcnt_cld_50_ID , 'comment', 'This is the 50th percentile or median of first_cloud_base, base height of first cloud layer detected, from 15 sec raw data over this time period.');
    vari=c.prcnt_cld_84;
        prcnt_cld_84_ID=add_mainVar_attributes_v2020(ncid,'cb_84_prct','double',time_dimID,'cloud base height from 84th percentile','cloud_base_altitude','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, prcnt_cld_84_ID , 'comment', 'This is the 84th percentile of first_cloud_base, base height of first cloud layer detected, from 15 sec raw data over this time period. Therefore this is one standard deviation above the median. ');



    %%% We are done defining the NetCdf
    netcdf.endDef(ncid);    

    %%% To be consistent with nomencalture, replace NaNs with -9999. Note: could use NaN but not conventional across all platforms
    c.frac_cloud(isnan(c.frac_cloud) == 1)=-9999;
    c.frac_clear(isnan(c.frac_clear) == 1)=-9999;
    c.count_cloud_total(isnan(c.count_cloud_total) == 1)=-9999;
    c.count_cloud_layers_1_to_3(isnan(c.count_cloud_layers_1_to_3) == 1)=-9999;
    c.count_cloud_obscured(isnan(c.count_cloud_obscured) == 1)=-9999;
    c.count_cloud_part_obsured_transparent(isnan(c.count_cloud_part_obsured_transparent) == 1)=-9999;
    c.count_clear(isnan(c.count_clear) == 1)=-9999;
    c.cloud_base(isnan(c.cloud_base) == 1)=-9999;
    c.prcnt_cld_5(isnan(c.prcnt_cld_5) == 1)=-9999;
    c.prcnt_cld_10(isnan(c.prcnt_cld_10) == 1)=-9999;
    c.prcnt_cld_16(isnan(c.prcnt_cld_16) == 1)=-9999;
    c.prcnt_cld_50(isnan(c.prcnt_cld_50) == 1)=-9999;
    c.prcnt_cld_84(isnan(c.prcnt_cld_84) == 1)=-9999;


else
    
    %%%%%% For 15-sec file
    
    %%% Define the main variables
    vari=c.ceil_b;  
        ceil_b_ID=add_mainVar_attributes_v2020(ncid,'ceil_b','double',[z_dimID, time_dimID],'backscatter coefficient from ceilometer','volume_attenuated_backwards_scattering_function_in_air','km-1 sr-1','physicalMeasurement',[-10., 10.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999.);
    vari=c.height;  
        height_ID=add_mainVar_attributes_v2020(ncid,'height','double',z_dimID,'height','','km','physicalMeasurement',[0 max(c.height)],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999.);
    vari=c.cloud_code; 
        cloud_code_ID=add_mainVar_attributes_v2020(ncid,'cloud_code','int',time_dimID,'cloud detection code','','','thematicClassification',[0, 5],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
        netcdf.putAtt(ncid, cloud_code_ID , 'comment', '0 = clear sky, 1 = one cloud layer detected, 2 = two cloud layers detected, 3 = three cloud layers detected, 4 = totally obscured, 5 = partly obscured but determined transparent');
    vari=c.cloud_base_first; 
        cloud1_ID=add_mainVar_attributes_v2020(ncid,'first_cloud_base','double',time_dimID,'base height of first cloud layer','','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
    vari=c.cloud_base_second;  
        cloud2_ID=add_mainVar_attributes_v2020(ncid,'second_cloud_base','double',time_dimID,'base height of second cloud layer','','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);
    vari=c.cloud_base_third;  
        cloud3_ID=add_mainVar_attributes_v2020(ncid,'third_cloud_base','double',time_dimID,'base height of third cloud layer','','km','physicalMeasurement',[0., 7.690],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999);

    %%% We are done defining the NetCdf
    netcdf.endDef(ncid);    

    %%% To be consistent with nomencalture, replace NaNs with -9999. Note: could use NaN but not conventional across all platforms
    c.ceil_b(isnan(c.ceil_b) == 1)=-9999; 
    c.height(isnan(c.height) == 1)=-9999; 
    c.cloud_base_first(isnan(c.cloud_base_first) == 1)=-9999; 
    c.cloud_base_second(isnan(c.cloud_base_second) == 1)=-9999; 
    c.cloud_base_third(isnan(c.cloud_base_third) == 1)=-9999; 
    c.cloud_code(isnan(c.cloud_code) == 1)=-9999; 

end


%%% get rid of nans

lat((isnan(lat) == 1))=-9999;
lon((isnan(lon) == 1))=-9999;    


%%% store the dimension variables in
% vari=(c.jd_ceil-1)*86400;
vari=c.te;
    netcdf.putVar(ncid,time_ID,vari);
vari=lat;
    netcdf.putVar(ncid,latitude_ID,vari);
vari=lon;
    netcdf.putVar(ncid,longitude_ID,vari);
    

if k == 1
        
    %%% store main variables
    vari=c.ceil_b;
        netcdf.putVar(ncid,ceil_b_ID,vari);     
    vari=c.height;
        netcdf.putVar(ncid,height_ID,vari);      
    vari=c.cloud_code;  
        netcdf.putVar(ncid,cloud_code_ID,vari);      
    vari=c.cloud_base_first;  
        netcdf.putVar(ncid,cloud1_ID,vari);      
    vari=c.cloud_base_second;  
        netcdf.putVar(ncid,cloud2_ID,vari);      
    vari=c.cloud_base_third;  
        netcdf.putVar(ncid,cloud3_ID,vari);  
    
else
        
    %%% store main variables
    vari=c.frac_cloud;
        netcdf.putVar(ncid,frac_cloud_ID,vari);      
    vari=c.frac_clear;
        netcdf.putVar(ncid,frac_clear_ID,vari);      
    vari=c.count_cloud_total;
        netcdf.putVar(ncid,count_total_ID,vari);      
    vari=c.count_cloud_layers_1_to_3;
        netcdf.putVar(ncid,count_layers_1_to_3_ID,vari);      
    vari=c.count_cloud_obscured;
        netcdf.putVar(ncid,count_cloud_obscured_ID,vari);      
    vari=c.count_cloud_part_obsured_transparent;
        netcdf.putVar(ncid,count_part_obsured_transparent_ID,vari);      
    vari=c.count_clear;
        netcdf.putVar(ncid,count_clear_ID,vari);      
    vari=c.cloud_base;
        netcdf.putVar(ncid,cloud_base_ID,vari);      
    vari=c.prcnt_cld_5;
        netcdf.putVar(ncid,prcnt_cld_5_ID,vari);      
    vari=c.prcnt_cld_10;
        netcdf.putVar(ncid,prcnt_cld_10_ID,vari); 
    vari=c.prcnt_cld_16;
        netcdf.putVar(ncid,prcnt_cld_16_ID,vari);     
    vari=c.prcnt_cld_50;
        netcdf.putVar(ncid,prcnt_cld_50_ID,vari);     
    vari=c.prcnt_cld_84;
        netcdf.putVar(ncid,prcnt_cld_84_ID,vari);     
end
    
%% We're done, close the netcdf
netcdf.close(ncid) 


ncdisp([path_out file_out]);

end % for k = 1, 2, 3 for different time scales

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

function Var_ID=add_mainVar_attributes_v2020(ncid,ID,datatype,dimIDs,long_name,standard_name,units,cdm_data_type,valid_range,actual_range,fillNvalue)

Var_ID = netcdf.defVar (ncid, ID, datatype, dimIDs);
    if ~isempty(long_name); netcdf.putAtt(ncid, Var_ID , 'long_name', long_name); end
    if ~isempty(standard_name);netcdf.putAtt(ncid, Var_ID , 'standard_name', standard_name);end
    if ~isempty(units);netcdf.putAtt(ncid, Var_ID , 'units', units);end
    if ~isempty(cdm_data_type);netcdf.putAtt(ncid, Var_ID , 'cdm_data_type', cdm_data_type);end    
    if ~isempty(valid_range);netcdf.putAtt(ncid, Var_ID , 'valid_range', valid_range);end
    if ~isempty(actual_range);netcdf.putAtt(ncid, Var_ID , 'actual_range', actual_range);end
    if ~isempty(fillNvalue);netcdf.defVarFill(ncid, Var_ID , false, fillNvalue); end
end