% Read mat file of interpolated soundings and output netcdf.
% (STRATUS 2007, SPdeS)

% preamble
read_parameters;
thermo_constants
latlon=false;

% construct title string, e.g. 'STRATUS 2007'
tit = cruise;
lt = length(cruise);
tit(lt+1) = ' ';
tit((lt+2):lt+5) = year; % This year is a string read in from the parameter file

% locate raw sounding files
aa = dir([way_raw_data_balloon '*_LATLON.txt']);
nsonde=length(aa);
% load yymmddhhMM dates from filenames
fst=str2mat(aa.name);

pint=(1030:-1:10)'; % pressure coordinate for interp. (hPa)
hint=NaN+zeros(length(pint),nsonde);
timeint=hint;
tint=hint;
tdint=hint;
thetaint=hint;
qint=hint;
rhint=hint;
uint=hint;
vint=hint;
wsint=hint;
wdint=hint;
latint=hint;
lonint=hint;

% convert filename timestamps to time variables
for son=1:nsonde;
    balfn=aa(son).name; % simple choice of file for now
    yr(son) = str2double(balfn(1:2))+2000;
    month(son) = str2double(balfn(3:4));
    day (son)  = str2double(balfn(5:6));
    hour(son)  = str2double(balfn(7:8));
    minute(son) = str2double(balfn(9:10));
end
yday=datenum(yr,month,day,hour,minute,0)-datenum(yr,0,0,0,0,0);

load([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2p.mat'])

for son=1:nsonde
    latlaunch(son)=latint(find(isfinite(latint(:,son)),1,'first'),son);
    lonlaunch(son)=lonint(find(isfinite(lonint(:,son)),1,'first'),son);
end

    % write yearly netcdf
    if latlon
        nc=netcdf([way_proc_data_balloon 'stratus' year 'sonde_p_latlon.nc'],'clobber');
    else
        nc=netcdf([way_proc_data_balloon 'stratus' year 'sonde_p.nc'],'clobber');
    end
        nc('launch')=0; % record dimension
    nc('pressure')=length(pint);
    
    nc{'year'}=ncint('launch');
    nc{'year'}.long_name='Year of the cruise';
    nc{'year'}.units='Gregorian year';
    nc{'year'}(1:nsonde)=yr;

    nc{'yday'}=ncfloat('launch');
    nc{'yday'}='year day of sonde';
    nc{'yday'}.units='days';
    nc{'yday'}(1:nsonde)=yday;
    
    nc{'month'}=ncshort('launch');
    nc{'month'}.long_name='calendar month number';
    nc{'month'}(1:nsonde)=month;

    nc{'day'}=ncshort('launch');
    nc{'day'}.long_name='day of month';
    nc{'day'}(1:nsonde)=day;

    nc{'hour'}=ncshort('launch');
    nc{'hour'}.long_name='hour UTC of launch';
    nc{'hour'}(1:nsonde)=hour;

    nc{'minute'}=ncshort('launch');
    nc{'minute'}.long_name='minute of hour of launch';
    nc{'minute'}(1:nsonde)=minute;

    if latlon
        nc.has_GPS_position_profile=true; % latlon files

        nc{'lat'}=ncfloat('launch','pressure');
        nc{'lat'}.long_name='latitude in flight'; % latlon files n/a
        nc{'lat'}.units='degrees latitude';
        nc{'lat'}(1:nsonde,:)=latint';

        nc{'lon'}=ncfloat('launch','pressure');
        nc{'lon'}.long_name='longitude in flight';
        nc{'lon'}.units='degrees longitude';
        nc{'lon'}(1:nsonde,:)=lonint';
    else
        nc.has_GPS_position_profile=false;

        nc{'lat'}=ncfloat('launch');
        nc{'lat'}.long_name='latitude of launch'; % latlon files n/a
        nc{'lat'}.units='degrees latitude';
        nc{'lat'}(1:nsonde)=latlaunch;
        
        nc{'lon'}=ncfloat('launch');
        nc{'lon'}.long_name='longitude of launch';
        nc{'lon'}.units='degrees longitude';
        nc{'lon'}(1:nsonde)=lonlaunch;
    end
    nc{'pres'}=ncfloat('pressure');
    nc{'pres'}.long_name='pressure'; % above mean sea level';
    nc{'pres'}.units='m';
    nc{'pres'}(:)=pint';

    % dependent variables...
%     nc{'time'}=ncfloat('launch','pressure');
%     nc{'time'}.long_name='time of flight since launch';
%     nc{'time'}.units='seconds';
%     nc{'time'}(1:nsonde,:)=timez';

    nc{'height'}=ncfloat('launch','pressure');
    nc{'height'}.long_name='altitude';
    nc{'height'}.units='m';
    nc{'height'}(1:nsonde,:)=hint';

    nc{'T'}=ncfloat('launch','pressure');
    nc{'T'}.long_name='temperature';
    nc{'T'}.units='degree C';
    nc{'T'}.source='observed';
    nc{'T'}(1:nsonde,:)=tint';

    nc{'RH'}=ncfloat('launch','pressure');
    nc{'RH'}.long_name='relative humidity';
    nc{'RH'}.units='percent';
    nc{'RH'}.source='observed';
    nc{'RH'}(1:nsonde,:)=rhint';

    nc{'Td'}=ncfloat('launch','pressure');
    nc{'Td'}.long_name='dew point temperature';
    nc{'Td'}.units='degree C';
    nc{'Td'}.source='Vaisala calculated from T and RH';
    nc{'Td'}(1:nsonde,:)=tdint';

    nc{'wnddir'}=ncfloat('launch','pressure');
    nc{'wnddir'}.long_name='wind direction';
    nc{'wnddir'}.units='degrees heading';
    nc{'wnddir'}.source='from interpolated u,v calculated from Vaisala/GPS';
    nc{'wnddir'}.comment='meteorological convention: 0 degree is northerly, 90 degre is easterly';
    nc{'wnddir'}.quality='poor due to interpolation';
    nc{'wnddir'}(1:nsonde,:)=wdint';

    nc{'wndspd'}=ncfloat('launch','pressure');
    nc{'wndspd'}.long_name='wind speed';
    nc{'wndspd'}.units='m/s';
    nc{'wndspd'}.source='Vaisala/GPS';
    nc{'wndspd'}(1:nsonde,:)=wsint';

    nc{'u'}=ncfloat('launch','pressure');
    nc{'u'}.long_name='wind zonal velocity';
    nc{'u'}.units='m/s';
    nc{'u'}.source='calculated from Vaisala wndspd and wnddir';
    nc{'u'}(1:nsonde,:)=uint';

    nc{'v'}=ncfloat('launch','pressure');
    nc{'v'}.long_name='wind meridional velocity';
    nc{'v'}.units='m/s';
    nc{'v'}.source='calculated from wndspd and wnddir';
    nc{'v'}(1:nsonde,:)=vint';

    nc{'theta'}=ncfloat('launch','pressure');
    nc{'theta'}.long_name='potential temperature';
    nc{'theta'}.units='K';
    nc{'theta'}.source='interpolated T and pres';
    nc{'theta'}(1:nsonde,:)=thetaint';
        
    nc{'q'}=ncfloat('launch','pressure');
    nc{'q'}.long_name='wind meridional velocity';
    nc{'q'}.units='g/kg';
    nc{'q'}.source='calculated from T and RH';
    nc{'q'}(1:nsonde,:)=qint';
    nc=close(nc);

