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

% preamble
read_parameters;
thermo_constants
quickplot=false;
%clear t;

% 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);

hz=0:10:25000;   % height coordinate (m)
pz=NaN+zeros(length(hz),length(aa));
tz=pz; tdz=pz; qz=pz; rhz=pz; uz=pz; vz=pz; latz=pz; lonz=pz; thetaz=pz;

% convert filename timestamps to time variables
for i=1:nsonde;
    balfn=aa(i).name; % simple choice of file for now
    yr(i) = str2double(balfn(1:2))+2000;
    month(i) = str2double(balfn(3:4));
    day (i)  = str2double(balfn(5:6));
    hour(i)  = str2double(balfn(7:8));
    minute(i) = 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_interp2z.mat'])
% MODIFY :: SPdeS, 4 March 2008
% Compress latz,lonz soundings to scalar station of launch to match other years.
for i=1:nsonde
    lat(i)=latz(find(latz(:,i)>-360,1),i); % take first as station lat, lon
    lon(i)=lonz(find(lonz(:,i)>-360,1),i); 
end
% also write lat,lon to netcdf rather than latz,lonz below
% END MODIFY :: SPdeS, 4 March 2008

    % write yearly netcdf
    nc=netcdf([way_proc_data_balloon 'stratus' year 'sonde_z.nc'],'clobber');
    nc('launch')=0; % record dimension
    nc('altitude')=length(hz);
    
    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;

    nc.has_GPS_position_profile=true; % latlon files

% MODIFY :: SPdeS 4 March 2008
%     nc{'lat'}=ncfloat('launch','altitude');
%     nc{'lat'}.long_name='latitude in flight'; % latlon files n/a
%     nc{'lat'}.units='degrees latitude';
%     nc{'lat'}(1:nsonde,:)=latz';
%     
%     nc{'lon'}=ncfloat('launch','altitude');
%     nc{'lon'}.long_name='longitude in flight';
%     nc{'lon'}.units='degrees longitude';
%     nc{'lon'}(1:nsonde,:)=lonz';

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

    nc{'lon'}=ncfloat('launch');
    nc{'lon'}.long_name='longitude';
    nc{'lon'}.units='degrees longitude';
    nc{'lon'}(1:nsonde)=lon;
% END MODIFY :: SPdeS 4 March 2008

    nc{'height'}=ncfloat('altitude');
    nc{'height'}.long_name='altitude'; % above mean sea level';
    nc{'height'}.units='m';
    nc{'height'}(:)=hz';

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

    nc{'pres'}=ncfloat('launch','altitude');
    nc{'pres'}.long_name='pressure';
    nc{'pres'}.units='hPa';
    nc{'pres'}(1:nsonde,:)=pz';

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

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

    nc{'Td'}=ncfloat('launch','altitude');
    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,:)=tdz';

    nc{'wnddir'}=ncfloat('launch','altitude');
    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,:)=wdz';

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

    nc{'u'}=ncfloat('launch','altitude');
    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,:)=uz';

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

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

