% gfd_ds_read_grid.m
% reads a dropsonde netcdf file from BADC
% interpolates onto a regular height grid using interp1
% then interpolates onto a distance grid 
% Creates an output mat file to be used with plotting code
% Original code by Ian Renfrew
% Modified by Stephen Outten 2007, and again by Ian Renfrew in Feb 2008

%------------------------------------------
% set up which dataset to use, etc.
%------------------------------------------
clear all;
iflight = 268;

if (iflight == 268)
    dir_name = 'F:\data\gfd\FAAM\B268\';
    dir_name_out = dir_name;
    input_leg = 4; % input_leg = 1, 2, 3, 4=all ds
    z_min = 0; z_max = 4500;    % min and max for vert coord. interpolation
elseif (iflight == 276)
    dir_name = 'C:\data\gfd\FAAM\B276\';
    dir_name_out = dir_name;
    input_leg = 14; 
end; 

flip_slice = [1 2]; % reverses the direction of these dropsonde legs
%------------------------------------------
% set input files
%------------------------------------------
if input_leg==1  
	infile{1} = [dir_name '210207_d1.nc'];
	infile{2} = [dir_name '210207_d2.nc'];
	infile{3} = [dir_name '210207_d4.nc'];  % Note: ds3 failed
    lat_dropsonde=[62 62 62];
    lon_dropsonde=[-40.0118 -40.735 -42.008];
    ofilename=([dir_name_out '210207_leg1.mat']);
elseif input_leg==2
    infile{1} = [dir_name '210207_d2.nc']; 
    infile{2} = [dir_name '210207_d5.nc']; 
	infile{3} = [dir_name '210207_d6.nc'];
	infile{4} = [dir_name '210207_d7.nc'];
	infile{5} = [dir_name '210207_d8.nc'];    
    infile{6} = [dir_name '210207_d10.nc']; 
    infile{7} = [dir_name '210207_d11.nc'];
    lat_dropsonde=[62       61.5   61.14   60.8    60.485  59.7385 59.651];     % launch location
    lon_dropsonde=[-40.735 -40.9  -41.188 -41.35  -41.5   -41.8725 -41.586 ];    % 6th lat/lon is average of 10-11
    ofilename=([dir_name_out '210207_leg2.mat']);
elseif input_leg==3
    infile{1} = [dir_name '210207_d9.nc'];
	infile{2} = [dir_name '210207_d10.nc'];
	infile{3} = [dir_name '210207_d11.nc'];
    load([dir_name '210207_ds12.mat']);    % synthetic dropsonde
    lat_dropsonde=[59.972   59.826   59.651  59.355];
    lon_dropsonde=[-42.995 -42.287  -41.586 -40.499];
    ofilename=([dir_name_out '210207_leg3.mat']);
elseif input_leg==4
	infile{1} = [dir_name '210207_d1.nc'];
	infile{2} = [dir_name '210207_d2.nc'];
	infile{3} = [dir_name '210207_d4.nc'];  
    infile{4} = [dir_name '210207_d5.nc'];
	infile{5} = [dir_name '210207_d6.nc'];    
    infile{6} = [dir_name '210207_d7.nc']; 
    infile{7} = [dir_name '210207_d8.nc'];
    infile{8} = [dir_name '210207_d9.nc'];
	infile{9} = [dir_name '210207_d10.nc'];    
    infile{10} = [dir_name '210207_d11.nc'];
    load([dir_name '210207_ds12_ds13.mat']);    % synthetic dropsondes
    %lat_dropsonde=[62 62 62 61.5   61.14   60.8    60.485 59.972   59.826   59.651];
    %lon_dropsonde=[-40.0118 -40.735 -42.008 -40.9  -41.188 -41.35  -41.5  -42.995 -42.287  -41.586];
    lat_dropsonde=[62 62 62 61.5   61.14   60.8    60.485 59.972   59.826   59.651 59.355 60.16];
    lon_dropsonde=[-40.0118 -40.735 -42.008 -40.9  -41.188 -41.35  -41.5  -42.995 -42.287  -41.586 -40.499 -40.29];
    ofilename=([dir_name_out '210207_leg4.mat']);
elseif input_leg==14
    infile{1} = [dir_name '050307_d1.nc'];
	infile{2} = [dir_name '050307_d2.nc'];
	infile{3} = [dir_name '050307_d3.nc'];
  	infile{4} = [dir_name '050307_d4.nc'];
   	infile{5} = [dir_name '050307_d5.nc'];
   	infile{6} = [dir_name '050307_d6.nc'];
   	infile{7} = [dir_name '050307_d7.nc'];
%    	infile{8} = [dir_name '050307_d8.nc'];
    lat_dropsonde=[65.77 65.57 65.37 65.19 64.99 64.80 64.61]; % 64.41];
    lon_dropsonde=[-33.97 -33.56 -33.20 -32.86 -32.51 -32.17 -31.83]; % -31.50];
    ofilename=([ dir_name_out '050307_leg1.mat']);
end

%--------------------------------------------------
% read in netCDF data, loop over no. of ds files 
%--------------------------------------------------
for counter = 1:length(infile)
	time_c = nc_varget(infile{counter},'time');
	date = datestr(datenum('01-Jan-1970') + floor(time_c(1)/86400));
	time = (time_c./86400) + datenum('01-Jan-1970') - datenum(date);
	hour=floor(time*24);
	min=floor((time*24-hour)*60); % this produces an approximate time stamp
	% sec = ((time*24-hour)*60-min)*60;
    for i = 1:length(min)
        if min(i)<10; 
            min_str(i,:)=['0' num2str(min(i))]; 
        else 
            min_str(i,:)=num2str(min(i)); 
        end;
    end
    time_str = [num2str(hour) min_str];
	data(:,1) = time_c(:);
	data(:,2) = str2num([datestr(date,'dd') datestr(date,'mm') datestr(date,'yy')]);
	data(:,3) = str2num(time_str);
	data(:,4) = nc_varget(infile{counter},'pres');
	data(:,5) = nc_varget(infile{counter},'tdry');
	data(:,6) = nc_varget(infile{counter},'dp');
	data(:,7) = nc_varget(infile{counter},'rh');
	data(:,8) = nc_varget(infile{counter},'u_wind');
	data(:,9) = nc_varget(infile{counter},'v_wind');
	data(:,10) = nc_varget(infile{counter},'dz');
	data(:,11) = nc_varget(infile{counter},'wspd');
	data(:,12) = nc_varget(infile{counter},'wdir');
	data(:,13) = nc_varget(infile{counter},'mr');
	data(:,14) = nc_varget(infile{counter},'theta');
	data(:,15) = nc_varget(infile{counter},'lat');
    data(:,16) = nc_varget(infile{counter},'lon');
    data(:,17) = nc_varget(infile{counter},'alt');
    data(:,18) = nc_varget(infile{counter},'theta_e');
    disp(['leg ds number =' num2str(counter) ' initial lat= ' num2str(data(10,15)) 'and lon= ' num2str(data(10,16))]);  
	data = flipud(data);
	[m n] = size(data);
    data = filterminus999(data);  
    data = data(2:m,:);
    ds{counter} = data;     % store each dropsonde in cell array ds
    clear data min_str
end
clear counter date hour i infile j m min n nc time time_c time_str 

%------------------------------------------------------------------------
% Extracts data to variables ds_p_var (etc); loop over no. of dropsondes
% This is done separately for thermodynamic variables and winds
% Note: need to convert mixing ratio (r, in g/kg) to spec. humidity q (in g/kg)
%--------------------------------------------------------------------------
for i=1:length(ds);
    count=1; % extract data to p, t, rh, q, theta also removes points if any of these are NaN
    for j=1:length(ds{i})
        if isnan(ds{i}(j,4))|isnan(ds{i}(j,5))|isnan(ds{i}(j,7))|isnan(ds{i}(j,13))|isnan(ds{i}(j,14))
        else
            ds_p_var{i}(count)=ds{i}(j,4);
            ds_t{i}(count)=ds{i}(j,5);
            ds_rh{i}(count)=ds{i}(j,7);
            ds_q{i}(count)=((ds{i}(j,13)/1000)./(1 + ds{i}(j,13)/1000))*1000;
            ds_theta{i}(count)=ds{i}(j,14);
            ds_thetae{i}(count)=ds{i}(j,18);
            ds_z{i}(count)=ds{i}(j,17);
            count=count+1;
        end
    end
end

for i=1:length(ds);
    count=1; % extract data to p, wd, & wsh, also removes points if any of these are NaN
    for j=1:length(ds{i})
        if isnan(ds{i}(j,4))|isnan(ds{i}(j,8))|isnan(ds{i}(j,9))|isnan(ds{i}(j,12))|isnan(ds{i}(j,11))|isnan(ds{i}(j,17))
        else
            ds_p_wind{i}(count)=ds{i}(j,4);
            ds_u{i}(count)=ds{i}(j,8);
            ds_v{i}(count)=ds{i}(j,9);
            ds_wd{i}(count)=ds{i}(j,12);
            ds_ws{i}(count)=ds{i}(j,11);
            ds_z_wind{i}(count)=ds{i}(j,17);
            count=count+1;
        end
    end
end
    
%---------------------------------------------------------------
% linear interpolation of profile onto regular HEIGHT (z) grid
%---------------------------------------------------------------
zi = z_min:10:z_max;         % vertical coordinate for grid - 10 m resolution
for i=1:length(lat_dropsonde)
	vert_t(:,i)=interp1(ds_z{i},ds_t{i},zi)';  
    vert_rh(:,i)=interp1(ds_z{i},ds_rh{i},zi)';  
    vert_q(:,i)=interp1(ds_z{i},ds_q{i},zi)';  
    vert_theta(:,i)=interp1(ds_z{i},ds_theta{i},zi)';  
    vert_thetae(:,i)=interp1(ds_z{i},ds_thetae{i},zi)';
	vert_p(:,i)=interp1(ds_z{i},ds_p_var{i},zi)';      
    vert_u(:,i)=interp1(ds_z_wind{i},ds_u{i},zi)';  
    vert_v(:,i)=interp1(ds_z_wind{i},ds_v{i},zi)';  
    vert_ws(:,i)=interp1(ds_z_wind{i},ds_ws{i},zi)';  
    vert_wd(:,i)=interp1(ds_z_wind{i},ds_wd{i},zi)';  
end

%----------------------------------------------------------------
% if B268, leg2, then average ds 10 and 11 to get end sounding
%----------------------------------------------------------------
if (input_leg == 2)
    vert_t(:,6) = mean([vert_t(:,6) vert_t(:,7)],2);  
    vert_rh(:,6) = mean([vert_rh(:,6) vert_rh(:,7)],2);  
    vert_q(:,6) = mean([vert_q(:,6) vert_q(:,7)],2);   
    vert_theta(:,6) = mean([vert_theta(:,6) vert_theta(:,7)],2);  
    vert_thetae(:,6) = mean([vert_thetae(:,6) vert_thetae(:,7)],2); 
    vert_p(:,6) = mean([vert_p(:,6) vert_p(:,7)],2);     
    vert_u(:,6) = mean([vert_u(:,6) vert_u(:,7)],2);  
    vert_v(:,6) = mean([vert_v(:,6) vert_v(:,7)],2);  
    vert_ws(:,6) = mean([vert_ws(:,6) vert_ws(:,7)],2);  
    vert_wd(:,6) = mean([vert_wd(:,6) vert_wd(:,7)],2);
    % snip off 7th column
    vert_t = vert_t(:,1:6);
    vert_rh = vert_rh(:,1:6);
    vert_q = vert_q(:,1:6);
    vert_theta = vert_theta(:,1:6);
    vert_thetae = vert_thetae(:,1:6);
    vert_p = vert_t(:,1:6);
    vert_u = vert_u(:,1:6);
    vert_v = vert_v(:,1:6);
    vert_ws = vert_ws(:,1:6);
    vert_wd = vert_wd(:,1:6);
    % snip off last lat/lon
    lat_dropsonde = lat_dropsonde(1:6);
    lon_dropsonde = lon_dropsonde(1:6);
end;

%-----------------------------------------------------------------
% second linear 1D interpolation onto a regular DISTANCE grid
%-----------------------------------------------------------------
lat_ds_rad=lat_dropsonde*pi/180;   % convert lat lon to radians
lon_ds_rad=lon_dropsonde*pi/180;
Re=6371;    
for i=2:length(lat_dropsonde)       % convert to distance along leg
    leg1_distance(i)=Re*(acos(sin(lat_ds_rad(1))*sin(lat_ds_rad(i))+cos(lat_ds_rad(1))*cos(lat_ds_rad(i))*cos(lon_ds_rad(1)-lon_ds_rad(i))));
end
leg1_distance(1)=0;
if find(flip_slice == input_leg)
    leg1_distance=fliplr(leg1_distance);  % reverses slice
end
slice_length=max(leg1_distance);
distance = [0:10:slice_length];             % horizontal coordinate of slice (10 km grid)
slice_t=interp1(leg1_distance,vert_t',distance,'linear')';
slice_rh=interp1(leg1_distance,vert_rh',distance,'linear')';
slice_q=interp1(leg1_distance,vert_q',distance,'linear')';
slice_theta=interp1(leg1_distance,vert_theta',distance,'linear')';
slice_thetae=interp1(leg1_distance,vert_thetae',distance,'linear')';
slice_p=interp1(leg1_distance,vert_p',distance,'linear')';
slice_wd = interp1(leg1_distance,vert_wd',distance,'linear')'; 
slice_ws = interp1(leg1_distance,vert_ws',distance,'linear')'; 
slice_u = interp1(leg1_distance,vert_u',distance,'linear')'; 
slice_v = interp1(leg1_distance,vert_v',distance,'linear')'; 

%-------------------------------------------------------
% calculate u_10N from "flight-level" dropsonde wind speed
% and also from flight-level profiles (filtered)
%-------------------------------------------------------
if (input_leg == 4) 
    q_ds_40m = vert_q(find(zi==40),:);
    t_ds_40m = vert_t(find(zi==40),:)+273.14;
    theta_ds_40m = vert_theta(find(zi==40),:);
    thetae_ds_40m = vert_thetae(find(zi==40),:);
    ws_ds_40m = vert_ws(find(zi==40),:);   
    ws_ds_50m = vert_ws(find(zi==50),:);    
    ws_ds_600m = vert_ws(find(zi==600),:);    
    ustar_B268 = 1.1245;    % mean from flight-level data
    %ws_ds_10N = ws_ds_40m - (ustar_B268/0.4)*log(40/10);
    ws_ds_10N = ws_ds_50m - (ustar_B268/0.4)*log(50/10);
end;

%----------------------------------------
% additional diagnostics 
%----------------------------------------
%gfd_ds_diagnostics
    
%-------------------------------------
% save output file in matlab format
%-------------------------------------
save (ofilename)