% 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-ddTHH:MM:SSZ');
else %means MST
    Tstring=datestr(now+datenum(0000,1,0,7,0,0),'yyyy-mm-ddTHH:MM:SSZ');
end    

%% Open the data file and .nc file to be created (netcdf4_classic format)
% data_in=load(['/Volumes/MyPassport/ATOMIC_2020/RHB/Ceilometer/Processed/all/ATOMIC_2020_ceilometer_15sec.mat']);
% temp=load('/Volumes/MyPassport/ATOMIC_2020/RHB/Ceilometer/Processed/backscatter/ceil_2019_0109.mat'); %to get range...
% % cloudbase=load(['/Volumes/MyPassport/ATOMIC_2020/RHB/Ceilometer/Processed/cloudbase/ATOMIC_2020_cloudbase_' sprintf('%03i',datenum(2020,m,d,0,0,0)-datenum(2020,0,0,0,0,0)) '.txt']);
% path_out='/Volumes/MyPassport/ATOMIC_2020/RHB/Ceilometer/Processed/all/';
% file_out=strcat(['EUREC4A_ATOMIC_RonBrown_Ceilometer_15s_20200109-20200212_v1.0.nc']);
data_in=load(['/Volumes/Data_1/PISTON_2019/Sally_Ride/Ceilometer/Processed_Mar2021/all/PISTON_2019_ceilometer_15sec.mat']);
 temp=load('/Volumes/Data_1/PISTON_2019/Sally_Ride/Ceilometer/Processed_Mar2021/all/PISTON_2019_ceilometer_0904_DOY_248.mat'); %to get range...
% flux_in=load('/Volumes/PISTON/PISTON_MISOBOB_2018/TGT/scientific_analysis/NetCDFfiles/PISTON_MISOBOB_2018_flux_10_v4_leg2n3.txt');  
gps=load('/Volumes/Data_1/PISTON_2019/Sally_Ride/flux/Field_Processed/gps/PISTON_2019_gpsnav_1min_all.txt');

% cloudbase=load(['/Volumes/MyPassport/ATOMIC_2020/RHB/Ceilometer/Processed/cloudbase/ATOMIC_2020_cloudbase_' sprintf('%03i',datenum(2020,m,d,0,0,0)-datenum(2020,0,0,0,0,0)) '.txt']);
path_out='/Volumes/Data_1/PISTON_2019/Sally_Ride/Ceilometer/Processed_Mar2021/all/';
file_out=strcat(['PISTON-Ceilometer-15s_RV-Sally-Ride_20190905_R0_thru_20190925.nc']);

cmode = netcdf.getConstant('NETCDF4');
cmode = bitor(cmode,netcdf.getConstant('CLASSIC_MODEL'));
ncid = netcdf.create(strcat(path_out,file_out),cmode);

jdy=data_in.cc.jd;
time_coverage_start =datestr(datenum(datestr('01-Jan-2019 00:00:00'))-1+jdy(1),0); %Time mark the start point of each 10min or hourly interval.
time_coverage_end =datestr(datenum(datestr('01-Jan-2019 00:00:00'))-1+jdy(end),0);
% lat=flux_in(:,34);
% lon=flux_in(:,35);
lat=gps(:,2);
lon=gps(:,3);
%interpolates to 15s
lati=interp1(gps(:,1),lat,jdy);
loni=interp1(gps(:,1),lon,jdy);


%% Define the dimensions (note in NetCDF, first dimension is 0)
[z_dimlen,time_dimlen]=size(data_in.cc.ceil_b);
% trajectory_dimlen=1;  %number of cruise legs
z_dimID= netcdf.defDim(ncid,'z',z_dimlen); 
time_dimID = netcdf.defDim(ncid,'time',time_dimlen);
% dimIDs = [time_dimID, z_dimID];%used to pass the dimids of the dimensions of the NETCDF variables. All the NETCDF variables we are creating share the same dimensions.
dimIDs = [z_dimID, time_dimID];%used to pass the dimids of the dimensions of the NETCDF variables. All the NETCDF variables we are creating share the same dimensions.

%% 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 2019-01-01 00:00:00UTC')
    netcdf.putAtt(ncid, time_ID, 'axis', 'T');
    netcdf.putAtt(ncid, time_ID, 'calendar', 'standard');
    netcdf.putAtt(ncid, time_ID, 'comment', 'Time mark 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, 'axis', 'Y');
    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, 'axis', 'X');
    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.);
z_ID=netcdf.defVar(ncid,'z','double',z_dimID);
    netcdf.putAtt(ncid, z_ID, 'long_name', 'height of each gate in meters');
    netcdf.putAtt(ncid, z_ID, 'standard_name', 'height');    
    netcdf.putAtt(ncid, z_ID, 'units', 'm');
    netcdf.putAtt(ncid, z_ID, 'axis', 'Z');
    netcdf.putAtt(ncid, z_ID, 'positive', 'up');
    netcdf.putAtt(ncid, z_ID, 'valid_min', min(temp.c.range));
    netcdf.putAtt(ncid, z_ID, 'valid_max', max(temp.c.range));
    netcdf.putAtt(ncid, z_ID , 'coverage_content_type', 'coordinate');
    netcdf.defVarFill(ncid, z_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 data collected from the RV Saly Ride by the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Lab (PSL) in collaboration with Oregon State University in the Western Pacific ocean from 2019-09-05 to 2019-09-25');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'summary','The data contained within this file covers about 20 days of volume backscatter and cloud base data in the West Pacific ocean. The data come from a Vaisala single-wavelength backscatter lidar (ceilometer) from NOAA PSL''s deployed on the RV Thomas G Thompson. Codes and heights were automatically computed by the Vaisala software. The averaging period is 15 seconds.') 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'keywords','Earth Science > Spectral/Engineering > Lidar > Lidar Backscatter , Earth Science > Atmosphere > Clouds > Cloud Properties > Cloud Base Height, 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'),'campaign','EUREC4A');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'instrument','Lidar Ceilometer from NOAA PSL');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'platform','RV Sally Ride');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'project','PISTON 2019');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'sea_name','West Pacific ocean');

netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'institution',sprintf('NOAA Physical Sciences Laboratory \n (325 Broadway R/PSD, Boulder, CO 80305)'));
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'naming_authority','gov.noaa.psl');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'id',file_out);

% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lat_min', min(data_in.lat)); %Describes a simple lower latitude limit. 
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lat_max', max(data_in.lat)); %Describes a simple upper latitude limit. 
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lon_min', min(data_in.lon)); %Describes a simple lower longitude limit. 
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'geospatial_lon_max', max(data_in.lon)); %Describes a simple upper longitude limit. 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'time_coverage_start', time_coverage_start); %Describes the time of the first data point in the data set. Use ISO 8601:2004 for date and time. 
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'time_coverage_end', time_coverage_end); % Describes the time of the last data point in the data set. Use ISO 8601:2004 for date and time.
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'time_coverage_duration', 'PT15.S') ; %Describes the duration of the data set. Use ISO 8601:2004 for date and time. PT10.S means a 10-second duration.

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','Ludovic Bariteau'); %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','ludovic.bariteau@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  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'),'product_version','v0');
netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'history','v0: original data');
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'history','v1.0: original data, v1.1: fixed offset in time, v1.2: added cloud base height, v1.2: all daily files merged into one');
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'processing_level','Quality controlled data');
% netcdf.putAtt(ncid,netcdf.getConstant('GLOBAL'),'comment','For reference, the height of the variables are specified in the attributes (coordinates) and can be used as coordinates. If not specified, the variables can be assumed to be at the surface. Bulk estimates of air sea fluxes were computed using the COARE bulk algorithm version 3.5 (Edson et al., 2013: On the exchange of momentum over the open ocean,J. Phys. Oceanogr., 43, 1589-1610). No information on the platform''s sensors (model, SN, calibration, etc...) were available yet and therefore are not included as variables'); 
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 cloud1 at 15-sec resolution) in a given hour or 10-min time period. The percentiles of cloud1 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 cloud1 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).'))

%% Define the platform, grid mapping, instruments, heights, and main variables
%platform
platform_ID=netcdf.defVar(ncid,'Thompson','char',time_dimID);
    netcdf.putAtt(ncid, platform_ID, 'long_name', 'RV Sally Ride');
    netcdf.putAtt(ncid, platform_ID, 'ncei_code', '33P4');
    netcdf.putAtt(ncid, platform_ID, 'call_sign', 'WSAF');
    netcdf.putAtt(ncid, platform_ID, 'imo_code', '9695171');
%grid mapping    
crs_ID=netcdf.defVar(ncid,'crs','double',[]);
    netcdf.putAtt(ncid, crs_ID, 'grid_mapping_name', 'latitude_longitude');
%     netcdf.putAtt(ncid, crs_ID, 'longitude_of_prime_meridian', '0.0d');
%     netcdf.putAtt(ncid, crs_ID, 'semi_major_axis', '6378137.0d');
%     netcdf.putAtt(ncid, crs_ID, 'inverse_flattening', '298.257223563d');
    netcdf.putAtt(ncid, crs_ID, 'epsg_code', 'EPSG:4326');
%instruments    
Ceilo_ID=netcdf.defVar(ncid,'Ceilometer_instrument','char',time_dimID);
    netcdf.putAtt(ncid, Ceilo_ID, 'long_name', 'Vaisala CL31 single-wavelength backscatter lidar (ceilometer)');
    netcdf.putAtt(ncid, Ceilo_ID, 'ncei_name', 'Ceilometer');
    netcdf.putAtt(ncid, Ceilo_ID, 'make_model', 'Vaisala CL31');
%     netcdf.putAtt(ncid, Ceilo_ID, 'serial_number', 'E2310001');
%             netcdf.putAtt(ncid, Ceilo_ID, 'calibration_date', 'unknown');  

%%%%%%
% Define the main variables
vari=data_in.cc.ceil_b;  
    ceil_b_ID=add_mainVar_attributes_v2020(ncid,'ceil_b','double',dimIDs,'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.,'ceilometer_instrument','Thompson','','crs','time lat lon z','time: point longitude: point latitude: point height:point');
%         [nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))] 
vari=data_in.cc.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.,'ceilometer_instrument','Thompson','','crs','time lat lon z','time: point longitude: point latitude: point height:point');
    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=data_in.cc.cloud_base_first; 
    cloud1_ID=add_mainVar_attributes_v2020(ncid,'first_cloud_base','double',time_dimID,'base height of first cloud layer','','m','thematicClassification',[0., 4000.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999.,'ceilometer_instrument','Thompson','','crs','time lat lon z','time: point longitude: point latitude: point height:point');
vari=data_in.cc.cloud_base_second;  
    cloud2_ID=add_mainVar_attributes_v2020(ncid,'second_cloud_base','double',time_dimID,'base height of second cloud layer','','m','thematicClassification',[0., 4000.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999.,'ceilometer_instrument','Thompson','','crs','time lat lon z','time: point longitude: point latitude: point height:point');
vari=data_in.cc.cloud_base_third;  
    cloud3_ID=add_mainVar_attributes_v2020(ncid,'third_cloud_base','double',time_dimID,'base height of third cloud layer','','m','thematicClassification',[0., 4000.],[nanmin(nanmin(vari,[],2)),nanmax(nanmax(vari,[],2))],-9999.,'ceilometer_instrument','Thompson','','crs','time lat lon z','time: point longitude: point latitude: point height:point');

%% 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
data_in.cc.ceil_b(isnan(data_in.cc.ceil_b))=-9999; 
data_in.cc.cloud_base_first(isnan(data_in.cc.cloud_base_first))=-9999; 
data_in.cc.cloud_base_second(isnan(data_in.cc.cloud_base_second))=-9999; 
data_in.cc.cloud_base_third(isnan(data_in.cc.cloud_base_third))=-9999; 
data_in.cc.cloud_code(isnan(data_in.cc.cloud_code))=-9999; 

lati(find(isnan(lati)))=-9999;    
loni(find(isnan(loni)))=-9999;    


%% Then store the dimension variables in
vari=data_in.cc.te;
    netcdf.putVar(ncid,time_ID,vari);
vari=lati;
    netcdf.putVar(ncid,latitude_ID,vari);
vari=loni;
    netcdf.putVar(ncid,longitude_ID,vari);
vari=temp.c.range;
    netcdf.putVar(ncid,z_ID,vari);

%% Then store main variables
vari=data_in.cc.ceil_b;
    netcdf.putVar(ncid,ceil_b_ID,vari);      
vari=data_in.cc.cloud_code;  
    netcdf.putVar(ncid,cloud_code_ID,vari);      
vari=data_in.cc.cloud_base_first;  
    netcdf.putVar(ncid,cloud1_ID,vari);      
vari=data_in.cc.cloud_base_second;  
    netcdf.putVar(ncid,cloud2_ID,vari);      
vari=data_in.cc.cloud_base_third;  
    netcdf.putVar(ncid,cloud3_ID,vari);      

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