POP_lat_transects_extract_IndoPac.m

function [] = POP_lat_transects_extract_IndoPac(input_filename)

% extract latitude transects from POP timestep file, and output to netCDF file


grid_res = 0.1;


output_file_form = 'pop.h.nday5.JC.transects.IndoPac_1deg';

lat_transects_vec = [-34.7; ((-34):1:70)'];
lon_bounds_array = [[20; 22; 27; 29; (30:0.4:32)'; 33; (35:(-0.25):34)'; (35.2:1.2:38.8)'; (40:(-0.2):38)'; (38.8:0.8:50)'; (43:(-1):37)'; (36.8:(-0.4):32)'; (121:(-0.5):117)'; (118:3.5:132)'; (136:1:140)'; (140*ones(4,1)); 136; (134.5:1.5:140.5)'; 154; 155; 162; 178; 177; 179; 180; 184; 182; 32; 28] (360 + [-71; ((-71):0.1:(-69.6))'; -70; -70; ((-71):(-1.5):(-75.5))'; ((-76):(-0.5):(-80))'; -80; -79; ((-80.5):0.7:(-77))'; ((-77)*ones(4,1)); -78; ((-83):(-1):(-87))'; -90.5; -92; -93; -99.5; -101.5; -104; ((-105)*ones(3,1)); -105.5; ((-106):(-1):(-112))'; -112.5; -116.5; -117; -118; ((-120):(-0.6):(-123))'; ((-123.5)*ones(6,1)); -122.5; -122; ((-122):(-1):(-128))'; ((-130):(-1.5):(-134.5))'; -138; -145; -165; -162; -160; ((-161)*ones(3,1)); -164; -84; -84])];

depth_range = [0 6000];


input_dayfile = input_filename(length(input_filename) + ((-12):1:(-3)));

longitude = ncread(input_filename,'TLONG');
land_ind = find(longitude < -1e-5);
longitude(land_ind) = max(max(max(longitude)));
latitude = ncread(input_filename,'TLAT');
latitude(land_ind) = max(max(max(latitude)));
level = 0.01*ncread(input_filename,'z_t');
time = ncread(input_filename,'time');
time_bound = ncread(input_filename,'time_bound');


longitude_vel = ncread(input_filename,'ULONG');
latitude_vel = ncread(input_filename,'ULAT');
% lev_wvel = 0.01*ncread(input_filename,'z_w_bot');
dz = 0.01*ncread(input_filename,'dz');

size_array = [size(longitude) length(level) 1];


start_lon_vec = NaN([1 length(lat_transects_vec)]);
grid_dir_array = NaN([size_array(1) length(lat_transects_vec)]);
dist_from_start_transect_min_array = NaN([size_array(1) length(lat_transects_vec)]);
dist_transect_array = NaN([size_array(1) length(lat_transects_vec)]);
lon_transect_array = NaN([size_array(1) length(lat_transects_vec)]);
lon_transect_min_array = NaN([size_array(1) length(lat_transects_vec)]);
lon_transect_max_array = NaN([size_array(1) length(lat_transects_vec)]);
dx_transect_array = NaN([size_array(1) length(lat_transects_vec)]);
angle_transect_array = NaN([size_array(1) length(lat_transects_vec)]);
dz_transect_array = NaN([size_array(1) length(lat_transects_vec) size_array(3)]);
transect_mask_array = NaN([size_array(1) length(lat_transects_vec) size_array(3)]);
vel_cross_transect_array = NaN([size_array(1) length(lat_transects_vec) size_array(3:4)]);
temp_transect_array = NaN([size_array(1) length(lat_transects_vec) size_array(3:4)]);
velT_cross_transect_array = NaN([size_array(1) length(lat_transects_vec) size_array(3:4)]);
for transect_ind = 1:length(lat_transects_vec)

    lat_transect = lat_transects_vec(transect_ind);
    lon_bounds = lon_bounds_array(transect_ind,:);
    
    size_array = [size(longitude) length(level) 1];
    
    if ((mod(lon_bounds(2) - lon_bounds(1) - 1e-5,360) + 1e-5 > mod(max(max(max(longitude))) - lon_bounds(1) - 1e-5,360) + 1e-5) || (mod(lon_bounds(2) - lon_bounds(1) - 1e-5,360) + 1e-5 > mod(lon_bounds(2) - min(min(min(longitude))) - 1e-5,360) + 1e-5))
        in_region_grid_ind = find(((longitude - (mod(lon_bounds(1) - min(min(min(longitude))),360) + min(min(min(longitude)))) >= 0) | (longitude - (mod(lon_bounds(2) - min(min(min(longitude))),360) + min(min(min(longitude)))) <= 0)) & (latitude - lat_transect >= (-2*grid_res)) & (latitude - lat_transect <= (2*grid_res)));
    else
        in_region_grid_ind = find(((longitude - (mod(lon_bounds(1) - min(min(min(longitude))),360) + min(min(min(longitude)))) >= 0) & (longitude - (mod(lon_bounds(2) - min(min(min(longitude))),360) + min(min(min(longitude)))) <= 0)) & (latitude - lat_transect >= (-2*grid_res)) & (latitude - lat_transect <= (2*grid_res)));
    end
    in_region_i_ind = mod(in_region_grid_ind - 1,size_array(1)) + 1;
    in_region_j_ind = ceil(in_region_grid_ind/size_array(1));    
%     if max(diff(unique(in_region_i_ind))) > 1.5
%         curr_unique = unique(in_region_i_ind);
%         [~,gap_ind] = max(diff(curr_unique));
%         region_i_ind = [(curr_unique(gap_ind + 1):1:size_array(1))'; (1:1:curr_unique(gap_ind))'];
%     else
        region_i_ind = unique(in_region_i_ind);
%     end
    
    start_marker_grid_ind = find((abs((mod(longitude - lon_bounds(1) + 180,360) - 180) - 0) < (0.51*grid_res)) & (latitude - lat_transect >= (-3*grid_res)) & (latitude - lat_transect <= (3*grid_res)));
    start_marker_i_ind = mod(start_marker_grid_ind - 1,size_array(1)) + 1;
    end_marker_grid_ind = find((abs((mod(longitude - lon_bounds(2) + 180,360) - 180) - (-0)) < (0.51*grid_res)) & (latitude - lat_transect >= (-3*grid_res)) & (latitude - lat_transect <= (3*grid_res)));
    end_marker_i_ind = mod(end_marker_grid_ind - 1,size_array(1)) + 1;
    if max(diff(unique(start_marker_i_ind))) > 1.5
        curr_unique = unique(start_marker_i_ind);
        [~,gap_ind] = max(diff(curr_unique));
        start_marker_i_ind = [(curr_unique(gap_ind + 1):1:size_array(1))'; (1:1:curr_unique(gap_ind))'];
    else
        start_marker_i_ind = unique(start_marker_i_ind);
    end
    if max(diff(unique(end_marker_i_ind))) > 1.5
        curr_unique = unique(end_marker_i_ind);
        [~,gap_ind] = max(diff(curr_unique));
        end_marker_i_ind = [(curr_unique(gap_ind + 1):1:size_array(1))'; (1:1:curr_unique(gap_ind))'];
    else
        end_marker_i_ind = unique(end_marker_i_ind);
    end
    
    curr_unique = unique([region_i_ind; start_marker_i_ind; end_marker_i_ind]);
    if diff(lon_bounds) > 359
        in_lon_range_ind = [start_marker_i_ind; ((start_marker_i_ind(length(start_marker_i_ind)) + 1):1:size_array(1))'; (1:1:(end_marker_i_ind(1) - 1))'; end_marker_i_ind];
    elseif max(diff(curr_unique)) > 2
        [~,gap_ind] = max(diff(curr_unique));
        if ((ismember(curr_unique(gap_ind + 1),start_marker_i_ind) == 1) || (isempty(find(ismember([1 size_array(1)],curr_unique) == 1,1)) == 0))
            in_lon_range_ind = [(curr_unique(gap_ind + 1):1:size_array(1))'; (1:1:curr_unique(gap_ind))'];
        else
            in_lon_range_ind = (curr_unique(1):1:curr_unique(length(curr_unique)))';
        end
    else
        in_lon_range_ind = (curr_unique(1):1:curr_unique(length(curr_unique)))';
    end
    
    
    % expand the index range to include a buffer
    n_i_buffer_ind = 2;
    n_j_buffer_ind = 2;
    in_lon_range_ind = [(mod(in_lon_range_ind(1) + ((-n_i_buffer_ind):1:(-1))' - 1,size_array(1)) + 1); in_lon_range_ind; (mod(in_lon_range_ind(length(in_lon_range_ind)) + (1:1:n_i_buffer_ind)' - 1,size_array(1)) + 1)]; %#ok<AGROW>
    in_lon_range_ind = in_lon_range_ind((in_lon_range_ind >= 1) & (in_lon_range_ind <= size_array(1)));
    in_lat_range_ind = [(min(in_region_j_ind) + ((-n_j_buffer_ind):1:(-1))'); unique(in_region_j_ind); (max(in_region_j_ind) + (1:1:n_j_buffer_ind)')];
    in_lat_range_ind = in_lat_range_ind((in_lat_range_ind >= 1) & (in_lat_range_ind <= size_array(2)));
    in_lat_range_ind = (min(in_lat_range_ind):1:max(in_lat_range_ind))';
    
    in_lev_range_ind = find((level - depth_range(1) >= -1e-5) & (level - depth_range(2) <= 1e-5));
    in_lev_range_ind = (max([1 (min(in_lev_range_ind) - 1)]):1:min([length(level) (max(in_lev_range_ind) + 2)]))';


    lon_in_range = longitude(in_lon_range_ind,in_lat_range_ind);
    lat_in_range = latitude(in_lon_range_ind,in_lat_range_ind);
    depth_in_range = double(level(in_lev_range_ind));
    
    lon_in_range_vel = longitude_vel(in_lon_range_ind,in_lat_range_ind);
    lat_in_range_vel = latitude_vel(in_lon_range_ind,in_lat_range_ind);
    if ismember(1,in_lev_range_ind) == 1
        depth_in_range_wvel_top = [0; double(cumsum(dz(in_lev_range_ind(1:(length(in_lev_range_ind) - 1)))))];
        depth_in_range_wvel = double(cumsum(dz(in_lev_range_ind)));
    else
        depth_in_range_wvel_top = double(sum(dz(1:(min(in_lev_range_ind) - 1)))) + [0; double(cumsum(dz(in_lev_range_ind(1:(length(in_lev_range_ind) - 1)))))];
        depth_in_range_wvel = double(sum(dz(1:(min(in_lev_range_ind) - 1))) + cumsum(dz(in_lev_range_ind)));
    end
    
    
    size_array = [length(in_lon_range_ind) length(in_lat_range_ind) length(in_lev_range_ind) 1];
    
    n_depth_subsets = ceil(prod(size_array([1 2]))/(1e6));
    size_depth_subset = ceil(size_array(3)/n_depth_subsets);
    
    
    bottom_depth_in_range = NaN([size_array(1:2) 1 1]);
    uvel_in_range = NaN(size_array);
    vvel_in_range = NaN(size_array);
    temp_in_range = NaN(size_array);
    uet_in_range = NaN(size_array);
    vnt_in_range = NaN(size_array);
    SSH_in_range = NaN([size_array(1:2) 1 size_array(4)]);
    hte = NaN([size_array(1:2) 1 1]);
    htn = NaN([size_array(1:2) 1 1]);
    dxu = NaN([size_array(1:2) 1 1]);
    dyu = NaN([size_array(1:2) 1 1]);
    tarea = NaN([size_array(1:2) 1 1]);
    angle = NaN([size_array(1:2) 1 1]);
    anglet = NaN([size_array(1:2) 1 1]);
    
    for depth_subset_ind = 1:ceil(size_array(3)/size_depth_subset)
        curr_subset_ind = ((((depth_subset_ind - 1)*size_depth_subset) + 1):1:min([(depth_subset_ind*size_depth_subset) size_array(3)]))';
        curr_in_lev_range_ind = in_lev_range_ind(curr_subset_ind);
        
        if ((length(find(ismember([1 size(longitude,1)],in_lon_range_ind) == 1)) > 1) && ((in_lon_range_ind(1) ~= 1) || (in_lon_range_ind(length(in_lon_range_ind)) ~= size(longitude,1))))
            if length(find(diff(in_lon_range_ind) >= 2)) >= 2
                disp('Error: multiple gaps')
            end
            large_gap_ind = find(diff(in_lon_range_ind) >= 2);
            if isempty(large_gap_ind) == 0
                lon_range_ind_span_1 = max(in_lon_range_ind) - in_lon_range_ind(large_gap_ind + 1) + 1;
                lon_range_ind_span_2 = in_lon_range_ind(large_gap_ind) - min(in_lon_range_ind) + 1;
                in_lon_range_ind = in_lon_range_ind([((large_gap_ind + 1):1:length(lon_in_range))'; (1:1:large_gap_ind)']);
                lon_in_range = lon_in_range([((large_gap_ind + 1):1:length(lon_in_range))'; (1:1:large_gap_ind)'],in_lat_range_ind);
                lon_in_range_vel = lon_in_range_vel([((large_gap_ind + 1):1:length(lon_in_range_vel))'; (1:1:large_gap_ind)'],in_lat_range_ind);
            end
    
            if isempty(find(diff(in_lon_range_ind) < 0,1)) == 0
                lon_range_ind_span_1 = max(in_lon_range_ind) - in_lon_range_ind(1) + 1;
                lon_range_ind_span_2 = in_lon_range_ind(length(in_lon_range_ind)) - min(in_lon_range_ind) + 1;
            end
            lat_range_ind_span = max(in_lat_range_ind) - min(in_lat_range_ind) + 1;
            depth_range_ind_span = max(curr_in_lev_range_ind) - min(curr_in_lev_range_ind) + 1;
    
            start_vec = [in_lon_range_ind(1) min(in_lat_range_ind) min(curr_in_lev_range_ind) 1];
            count_vec = [lon_range_ind_span_1 lat_range_ind_span depth_range_ind_span 1];
            
            curr_uvel_in_range = 0.01*double(ncread(input_filename,'UVEL',start_vec,count_vec));
            curr_vvel_in_range = 0.01*double(ncread(input_filename,'VVEL',start_vec,count_vec));
            curr_temp_in_range = double(ncread(input_filename,'TEMP',start_vec,count_vec));
            curr_uet_in_range = double(ncread(input_filename,'UET',start_vec,count_vec));
            curr_vnt_in_range = double(ncread(input_filename,'VNT',start_vec,count_vec));
            if depth_subset_ind == 1
                curr_bottom_depth_in_range = 0.01*double(ncread(input_filename,'HT',start_vec(1:2),count_vec(1:2)));
                curr_SSH_in_range = reshape(0.01*double(ncread(input_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))),[count_vec(1:2) 1 count_vec(4)]);
                curr_hte_in_range = 0.01*double(ncread(input_filename,'HTE',start_vec(1:2),count_vec(1:2)));
                curr_htn_in_range = 0.01*double(ncread(input_filename,'HTN',start_vec(1:2),count_vec(1:2)));
                curr_dxu_in_range = 0.01*double(ncread(input_filename,'DXU',start_vec(1:2),count_vec(1:2)));
                curr_dyu_in_range = 0.01*double(ncread(input_filename,'DYU',start_vec(1:2),count_vec(1:2)));
                curr_tarea_in_range = (1e-4)*double(ncread(input_filename,'TAREA',start_vec(1:2),count_vec(1:2)));
                curr_angle_in_range = double(ncread(input_filename,'ANGLE',start_vec(1:2),count_vec(1:2)));
                curr_anglet_in_range = double(ncread(input_filename,'ANGLET',start_vec(1:2),count_vec(1:2)));
            end
    
            start_vec = [min(in_lon_range_ind) min(in_lat_range_ind) min(curr_in_lev_range_ind) 1];
            count_vec = [lon_range_ind_span_2 lat_range_ind_span depth_range_ind_span 1];
    
            curr_uvel_in_range_added = 0.01*double(ncread(input_filename,'UVEL',start_vec,count_vec));
            curr_vvel_in_range_added = 0.01*double(ncread(input_filename,'VVEL',start_vec,count_vec));
            curr_temp_in_range_added = double(ncread(input_filename,'TEMP',start_vec,count_vec));
            curr_uet_in_range_added = double(ncread(input_filename,'UET',start_vec,count_vec));
            curr_vnt_in_range_added = double(ncread(input_filename,'VNT',start_vec,count_vec));
            if depth_subset_ind == 1
                curr_bottom_depth_in_range_added = 0.01*double(ncread(input_filename,'HT',start_vec(1:2),count_vec(1:2)));
                curr_SSH_in_range_added = reshape(0.01*double(ncread(input_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))),[count_vec(1:2) 1 count_vec(4)]);
                curr_hte_in_range_added = 0.01*double(ncread(input_filename,'HTE',start_vec(1:2),count_vec(1:2)));
                curr_htn_in_range_added = 0.01*double(ncread(input_filename,'HTN',start_vec(1:2),count_vec(1:2)));
                curr_dxu_in_range_added = 0.01*double(ncread(input_filename,'DXU',start_vec(1:2),count_vec(1:2)));
                curr_dyu_in_range_added = 0.01*double(ncread(input_filename,'DYU',start_vec(1:2),count_vec(1:2)));
                curr_tarea_in_range_added = (1e-4)*double(ncread(input_filename,'TAREA',start_vec(1:2),count_vec(1:2)));
                curr_angle_in_range_added = double(ncread(input_filename,'ANGLE',start_vec(1:2),count_vec(1:2)));
                curr_anglet_in_range_added = double(ncread(input_filename,'ANGLET',start_vec(1:2),count_vec(1:2)));
            end
    
            curr_uvel_in_range = [curr_uvel_in_range; curr_uvel_in_range_added]; %#ok<AGROW>
            curr_uvel_in_range(curr_uvel_in_range > 1e15) = 0;
            clear curr_uvel_in_range_added
            uvel_in_range(:,:,curr_subset_ind,1) = curr_uvel_in_range(:,:,:,1);
            curr_vvel_in_range = [curr_vvel_in_range; curr_vvel_in_range_added]; %#ok<AGROW>
            curr_vvel_in_range(curr_vvel_in_range > 1e15) = 0;
            clear curr_vvel_in_range_added
            vvel_in_range(:,:,curr_subset_ind,1) = curr_vvel_in_range(:,:,:,1);
            curr_temp_in_range = [curr_temp_in_range; curr_temp_in_range_added]; %#ok<AGROW>
            curr_temp_in_range(curr_temp_in_range > 1e15) = 0;
            clear curr_temp_in_range_added
            temp_in_range(:,:,curr_subset_ind,1) = curr_temp_in_range(:,:,:,1);
            curr_uet_in_range = [curr_uet_in_range; curr_uet_in_range_added]; %#ok<AGROW>
            curr_uet_in_range(curr_uet_in_range > 1e15) = 0;
            clear curr_uet_in_range_added
            uet_in_range(:,:,curr_subset_ind,1) = curr_uet_in_range(:,:,:,1);
            curr_vnt_in_range = [curr_vnt_in_range; curr_vnt_in_range_added]; %#ok<AGROW>
            curr_vnt_in_range(curr_vnt_in_range > 1e15) = 0;
            clear curr_vnt_in_range_added
            vnt_in_range(:,:,curr_subset_ind,1) = curr_vnt_in_range(:,:,:,1);
            if depth_subset_ind == 1
                curr_bottom_depth_in_range = [curr_bottom_depth_in_range; curr_bottom_depth_in_range_added]; %#ok<AGROW>
                curr_bottom_depth_in_range(curr_bottom_depth_in_range > 1e15) = 0;
                clear curr_bottom_depth_in_range_added
                bottom_depth_in_range(:,:,1,1) = curr_bottom_depth_in_range;
                curr_SSH_in_range = [curr_SSH_in_range; curr_SSH_in_range_added]; %#ok<AGROW>
                curr_SSH_in_range((curr_SSH_in_range > 1e15) | (isnan(curr_SSH_in_range) == 1)) = 0;
                clear curr_SSH_in_range_added
                SSH_in_range(:,:,1,1) = curr_SSH_in_range(:,:,:,1);
                curr_hte_in_range = [curr_hte_in_range; curr_hte_in_range_added]; %#ok<AGROW>
                curr_hte_in_range(curr_hte_in_range > 1e15) = 0;
                clear curr_hte_in_range_added
                hte(:,:,1,1) = curr_hte_in_range;
                curr_htn_in_range = [curr_htn_in_range; curr_htn_in_range_added]; %#ok<AGROW>
                curr_htn_in_range(curr_htn_in_range > 1e15) = 0;
                clear curr_htn_in_range_added
                htn(:,:,1,1) = curr_htn_in_range;
                curr_dxu_in_range = [curr_dxu_in_range; curr_dxu_in_range_added]; %#ok<AGROW>
                curr_dxu_in_range(curr_dxu_in_range > 1e15) = 0;
                clear curr_dxu_in_range_added
                dxu(:,:,1,1) = curr_dxu_in_range;
                curr_dyu_in_range = [curr_dyu_in_range; curr_dyu_in_range_added]; %#ok<AGROW>
                curr_dyu_in_range(curr_dyu_in_range > 1e15) = 0;
                clear curr_dyu_in_range_added
                dyu(:,:,1,1) = curr_dyu_in_range;
                curr_tarea_in_range = [curr_tarea_in_range; curr_tarea_in_range_added]; %#ok<AGROW>
                curr_tarea_in_range(curr_tarea_in_range > 1e15) = 0;
                clear curr_tarea_in_range_added
                tarea(:,:,1,1) = curr_tarea_in_range;
                curr_angle_in_range = [curr_angle_in_range; curr_angle_in_range_added]; %#ok<AGROW>
                curr_angle_in_range((curr_angle_in_range > 1e15) | (isnan(curr_angle_in_range) == 1)) = 0;
                clear curr_angle_in_range_added
                angle(:,:,1,1) = curr_angle_in_range;
                curr_anglet_in_range = [curr_anglet_in_range; curr_anglet_in_range_added]; %#ok<AGROW>
                curr_anglet_in_range((curr_anglet_in_range > 1e15) | (isnan(curr_anglet_in_range) == 1)) = 0;
                clear curr_anglet_in_range_added
                anglet(:,:,1,1) = curr_anglet_in_range;
            end
        else
            lon_range_ind_span = max(in_lon_range_ind) - min(in_lon_range_ind) + 1;
            lat_range_ind_span = max(in_lat_range_ind) - min(in_lat_range_ind) + 1;
            depth_range_ind_span = max(curr_in_lev_range_ind) - min(curr_in_lev_range_ind) + 1;
    
            start_vec = [min(in_lon_range_ind) min(in_lat_range_ind) min(curr_in_lev_range_ind) 1];
            count_vec = [lon_range_ind_span lat_range_ind_span depth_range_ind_span 1];
    
            curr_uvel_in_range = 0.01*double(ncread(input_filename,'UVEL',start_vec,count_vec));
            curr_vvel_in_range = 0.01*double(ncread(input_filename,'VVEL',start_vec,count_vec));
            curr_temp_in_range = double(ncread(input_filename,'TEMP',start_vec,count_vec));
            curr_uet_in_range = double(ncread(input_filename,'UET',start_vec,count_vec));
            curr_vnt_in_range = double(ncread(input_filename,'VNT',start_vec,count_vec));
            curr_uvel_in_range(curr_uvel_in_range > 1e15) = 0;
            curr_vvel_in_range(curr_vvel_in_range > 1e15) = 0;
            curr_temp_in_range(curr_temp_in_range > 1e15) = 0;
            curr_uet_in_range(curr_uet_in_range > 1e15) = 0;
            curr_vnt_in_range(curr_vnt_in_range > 1e15) = 0;
            uvel_in_range(:,:,curr_subset_ind,1) = curr_uvel_in_range(:,:,:,1);
            vvel_in_range(:,:,curr_subset_ind,1) = curr_vvel_in_range(:,:,:,1);
            temp_in_range(:,:,curr_subset_ind,1) = curr_temp_in_range(:,:,:,1);
            uet_in_range(:,:,curr_subset_ind,1) = curr_uet_in_range(:,:,:,1);
            vnt_in_range(:,:,curr_subset_ind,1) = curr_vnt_in_range(:,:,:,1);
            
            if depth_subset_ind == 1
                curr_bottom_depth_in_range = 0.01*double(ncread(input_filename,'HT',start_vec(1:2),count_vec(1:2)));
                curr_SSH_in_range = reshape(0.01*double(ncread(input_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))),[count_vec(1:2) 1 count_vec(4)]);
                curr_hte_in_range = 0.01*double(ncread(input_filename,'HTE',start_vec(1:2),count_vec(1:2)));
                curr_htn_in_range = 0.01*double(ncread(input_filename,'HTN',start_vec(1:2),count_vec(1:2)));
                curr_dxu_in_range = 0.01*double(ncread(input_filename,'DXU',start_vec(1:2),count_vec(1:2)));
                curr_dyu_in_range = 0.01*double(ncread(input_filename,'DYU',start_vec(1:2),count_vec(1:2)));
                curr_tarea_in_range = (1e-4)*double(ncread(input_filename,'TAREA',start_vec(1:2),count_vec(1:2)));
                curr_angle_in_range = double(ncread(input_filename,'ANGLE',start_vec(1:2),count_vec(1:2)));
                curr_anglet_in_range = double(ncread(input_filename,'ANGLET',start_vec(1:2),count_vec(1:2)));
                curr_bottom_depth_in_range(isnan(curr_bottom_depth_in_range) == 1) = 0;
                bottom_depth_in_range(:,:,1,1) = curr_bottom_depth_in_range;
                curr_SSH_in_range((curr_SSH_in_range < -999) | (isnan(curr_SSH_in_range) == 1)) = 0;
                SSH_in_range(:,:,1,1) = curr_SSH_in_range(:,:,:,1);
                curr_hte_in_range(isnan(curr_hte_in_range) == 1) = 0;
                hte(:,:,1,1) = curr_hte_in_range;
                curr_htn_in_range(isnan(curr_htn_in_range) == 1) = 0;
                htn(:,:,1,1) = curr_htn_in_range;
                curr_dxu_in_range(isnan(curr_dxu_in_range) == 1) = 0;
                dxu(:,:,1,1) = curr_dxu_in_range;
                curr_dyu_in_range(isnan(curr_dyu_in_range) == 1) = 0;
                dyu(:,:,1,1) = curr_dyu_in_range;
                curr_tarea_in_range(isnan(curr_tarea_in_range) == 1) = 0;
                tarea(:,:,1,1) = curr_tarea_in_range;
                curr_angle_in_range((curr_angle_in_range > 1e15) | (isnan(curr_angle_in_range) == 1)) = 0;
                angle(:,:,1,1) = curr_angle_in_range;
                curr_anglet_in_range((curr_anglet_in_range > 1e15) | (isnan(curr_anglet_in_range) == 1)) = 0;
                anglet(:,:,1,1) = curr_anglet_in_range;
            end
        end
    end
    
    angle_imax = NaN(size(angle));
    angle_imax(:,2:size_array(2),1,1) = angle(:,2:size_array(2),1,1) - (diff(angle,1,2)/2);
    angle_jmax = NaN(size(angle));
    angle_jmax(2:size_array(1),:,1,1) = angle(2:size_array(1),:,1,1) - (diff(angle,1,1)/2);
    
    
    clear curr*
    
    
    % adjust longitude values to account for branch cut
    diff_lon = diff(lon_in_range,1,1);
    diff_lon(abs(diff_lon) > 100) = mod(diff_lon(abs(diff_lon) > 100) + 180,360) - 180;
    lon_in_range = repmat(lon_in_range(1,:) + (360*round((lon_bounds(1) - lon_in_range(1,:))/360)),[size(lon_in_range,1) 1]) + [zeros([1 size(lon_in_range,2)]); cumsum(diff_lon,1,'forward')];
    diff_lon = diff(lon_in_range_vel,1,1);
    diff_lon(abs(diff_lon) > 100) = mod(diff_lon(abs(diff_lon) > 100) + 180,360) - 180;
    lon_in_range_vel = repmat(lon_in_range_vel(1,:) + (360*round((lon_bounds(1) - lon_in_range_vel(1,:))/360)),[size(lon_in_range_vel,1) 1]) + [zeros([1 size(lon_in_range_vel,2)]); cumsum(diff_lon,1,'forward')];
    
    
    if ismember(1,in_lev_range_ind) == 1
        depth_in_range_wvel_withtop = [0; depth_in_range_wvel];
    else
        depth_in_range_wvel_withtop = [sum(dz(1:(min(in_lev_range_ind) - 1))); depth_in_range_wvel];
    end
    
    depth_upper_array = depth_range(1)*ones([size_array(1:2) 1]);
    depth_lower_array = depth_range(2)*ones([size_array(1:2) 1]);
    
    depth_upper_bound_interp_ind = interp1(depth_in_range_wvel_withtop,1:1:length(depth_in_range_wvel_withtop),reshape(depth_upper_array,[numel(depth_upper_array) 1]));
    depth_upper_array = reshape(depth_in_range_wvel_withtop(max([round(depth_upper_bound_interp_ind) ones([length(depth_upper_bound_interp_ind) 1])],[],2)),size(depth_upper_array));
    depth_lower_bound_interp_ind = interp1(depth_in_range_wvel_withtop,1:1:length(depth_in_range_wvel_withtop),reshape(depth_lower_array,[numel(depth_lower_array) 1]));
    depth_lower_array = reshape(depth_in_range_wvel_withtop(min([round(depth_lower_bound_interp_ind) length(depth_in_range_wvel_withtop)*ones([length(depth_upper_bound_interp_ind) 1])],[],2)),size(depth_lower_array));
    
    
    bottom_depth_minus_cell_top = repmat(bottom_depth_in_range,[1 1 size_array(3)]) - repmat(reshape(depth_in_range_wvel_withtop(1:size_array(3)),[1 1 size_array(3)]),[size_array(1:2) 1]);
    dzt = repmat(reshape(dz(in_lev_range_ind),[1 1 size_array(3)]),[size_array(1:2) 1]);
    dzt = reshape(min([reshape(dzt,[prod(size_array(1:3)) 1]) reshape(bottom_depth_minus_cell_top,[prod(size_array(1:3)) 1])],[],2),size_array(1:3));
    dzt(dzt < 0) = 0;
    % if ismember(1,in_lev_range_ind) == 1
    %     dzt(:,:,1,:) = dzt(:,:,1,:) + SSH_in_range;
    % end
    dzu = zeros(size_array(1:3));
    dzu(1:(size_array(1) - 1),:,:) = reshape(min([reshape(dzt(1:(size_array(1) - 1),:,:),[prod(size_array(1:3) - [1 0 0]) 1]) reshape(dzt(2:size_array(1),:,:),[prod(size_array(1:3) - [1 0 0]) 1])],[],2),size_array(1:3) - [1 0 0]);
    dzv = zeros(size_array(1:3));
    dzv(:,1:(size_array(2) - 1),:) = reshape(min([reshape(dzt(:,1:(size_array(2) - 1),:),[prod(size_array(1:3) - [0 1 0]) 1]) reshape(dzt(:,2:size_array(2),:),[prod(size_array(1:3) - [0 1 0]) 1])],[],2),size_array(1:3) - [0 1 0]);
    dzu_ugrid = zeros(size(dzt));
    dzu_ugrid(:,1:(size_array(2) - 1),:) = reshape(min([reshape(dzu(:,1:(size_array(2) - 1),:),[prod(size_array(1:3) - [0 1 0]) 1]) reshape(dzu(:,2:size_array(2),:),[prod(size_array(1:3) - [0 1 0]) 1])],[],2),size_array(1:3) - [0 1 0]);
    dzv_ugrid = zeros(size(dzt));
    dzv_ugrid(1:(size_array(1) - 1),:,:) = reshape(min([reshape(dzv(1:(size_array(1) - 1),:,:),[prod(size_array(1:3) - [1 0 0]) 1]) reshape(dzv(2:size_array(1),:,:),[prod(size_array(1:3) - [1 0 0]) 1])],[],2),size_array(1:3) - [1 0 0]);
    
    
    % compute (weighted) interpolated flux velocities
    uvel_interp = NaN(size_array);
    uvel_in_range(isnan(uvel_in_range)) = 0;
    uvel_interp(:,2:size_array(2),:,:) = ((0.5*repmat(repmat(dyu(:,1:(size_array(2) - 1)),[1 1 size_array(3)]).*dzu_ugrid(:,1:(size_array(2) - 1),:),[1 1 1 size_array(4)]).*uvel_in_range(:,1:(size_array(2) - 1),:,:)) + (0.5*repmat(repmat(dyu(:,2:size_array(2)),[1 1 size_array(3)]).*dzu_ugrid(:,2:size_array(2),:),[1 1 1 size_array(4)]).*uvel_in_range(:,2:size_array(2),:,:)))./repmat(repmat(hte(:,2:size_array(2)),[1 1 size_array(3)]).*dzu(:,2:size_array(2),:),[1 1 1 size_array(4)]);
    uvel_interp(isnan(uvel_interp) == 1) = 0;
    vvel_interp = NaN(size_array);
    vvel_in_range(isnan(vvel_in_range)) = 0;
    vvel_interp(2:size_array(1),:,:,:) = ((0.5*repmat(repmat(dxu(1:(size_array(1) - 1),:),[1 1 size_array(3)]).*dzv_ugrid(1:(size_array(1) - 1),:,:),[1 1 1 size_array(4)]).*vvel_in_range(1:(size_array(1) - 1),:,:,:)) + (0.5*repmat(repmat(dxu(2:size_array(1),:),[1 1 size_array(3)]).*dzv_ugrid(2:size_array(1),:,:),[1 1 1 size_array(4)]).*vvel_in_range(2:size_array(1),:,:,:)))./repmat(repmat(htn(2:size_array(1),:),[1 1 size_array(3)]).*dzv(2:size_array(1),:,:),[1 1 1 size_array(4)]);
    vvel_interp(isnan(vvel_interp) == 1) = 0;
    
    % compute tracer fluxes in correct units
    uT_imax = (repmat(repmat(tarea,[1 1 size_array(3)]).*dzt,[1 1 1 size_array(4)]).*uet_in_range)./repmat(repmat(hte,[1 1 size_array(3)]).*dzu,[1 1 1 size_array(4)]);
    uT_imax((isnan(uT_imax) == 1) | (isinf(uT_imax) == 1)) = 0;
    vT_jmax = (repmat(repmat(tarea,[1 1 size_array(3)]).*dzt,[1 1 1 size_array(4)]).*vnt_in_range)./repmat(repmat(htn,[1 1 size_array(3)]).*dzv,[1 1 1 size_array(4)]);
    vT_jmax((isnan(vT_jmax) == 1) | (isinf(vT_jmax) == 1)) = 0;
    
    clear uet_in_range vnt_in_range
    
    
    % interpolate from T-grid values
    temp_interp_x = NaN(size_array);
    temp_interp_x(1:(size_array(1) - 1),:,:,:) = temp_in_range(2:size_array(1),:,:,:) - (diff(temp_in_range,1,1)/2);
    temp_interp_x(isnan(temp_interp_x) == 1) = 0;
    temp_interp_y = NaN(size_array);
    temp_interp_y(:,1:(size_array(2) - 1),:,:,:) = temp_in_range(:,2:size_array(2),:,:) - (diff(temp_in_range,1,2)/2);
    temp_interp_y(isnan(temp_interp_y) == 1) = 0;
    
    
    % transect masks
    
    in_region_ind_1 = find((lon_in_range - (lon_bounds(1) + (360*round((min(min(lon_in_range)) - lon_bounds(1))/360))) >= 0) & (lon_in_range - (lon_bounds(2) + (360*round((max(max(lon_in_range)) - lon_bounds(2))/360))) < 0) & (lat_in_range - lat_transect >= (-2*grid_res)) & (lat_in_range - lat_transect < 0));
    in_region_mask_2D = zeros(size_array(1:2));
    in_region_mask_2D(in_region_ind_1) = -0.5; %#ok<FNDSB>
    in_region_ind_2 = find((lon_in_range - (lon_bounds(1) + (360*round((min(min(lon_in_range)) - lon_bounds(1))/360))) >= 0) & (lon_in_range - (lon_bounds(2) + (360*round((max(max(lon_in_range)) - lon_bounds(2))/360))) < 0) & (lat_in_range - lat_transect >= 0) & (lat_in_range - lat_transect < (2*grid_res)));
    in_region_mask_2D(in_region_ind_2) = 0.5; %#ok<FNDSB>
    
    depth_upper_cell_top = reshape(max([reshape(repmat(depth_upper_array,[1 1 size_array(3)]),[prod(size_array(1:3)) 1]) reshape(repmat(reshape(depth_in_range_wvel_withtop(1:size_array(3)),[1 1 size_array(3)]),[size_array(1:2) 1]),[prod(size_array(1:3)) 1])],[],2),size_array(1:3));
    depth_lower_cell_bottom = reshape(min([reshape(repmat(reshape(min([reshape(depth_lower_array,[prod(size_array([1 2])) 1]) reshape(bottom_depth_in_range,[prod(size_array([1 2])) 1])],[],2),[prod(size_array(1:2)) 1]),[size_array(3) 1]),[prod(size_array(1:3)) 1]) reshape(repmat(reshape(depth_in_range_wvel_withtop(2:(size_array(3) + 1)),[1 1 size_array(3)]),[size_array(1:2) 1]),[prod(size_array(1:3)) 1])],[],2),size_array(1:3));
    in_region_mask = repmat(in_region_mask_2D,[1 1 size_array(3)]);
    in_region_mask(depth_lower_cell_bottom - depth_upper_cell_top < 1e-5) = 0;
    
    transect_mask_x = zeros(size_array(1:3));
    transect_mask_x(1:(size_array(1) - 1),:,:) = diff(in_region_mask,1,1);
    transect_mask_x(abs(transect_mask_x) < 0.9) = 0;
    transect_mask_y = zeros(size_array(1:3));
    transect_mask_y(:,1:(size_array(2) - 1),:) = diff(in_region_mask,1,2);
    transect_mask_y(abs(transect_mask_y) < 0.9) = 0;
    
    lon_grid_imax = NaN(size_array(1:2));
    lat_grid_imax = NaN(size_array(1:2));
    lon_grid_imax(1:(size_array(1) - 1),:) = lon_in_range(2:size_array(1),:) - (diff(lon_in_range,1,1)/2);
    lat_grid_imax(1:(size_array(1) - 1),:) = lat_in_range(2:size_array(1),:) - (diff(lat_in_range,1,1)/2);
    lon_grid_jmax = NaN(size_array(1:2));
    lat_grid_jmax = NaN(size_array(1:2));
    lon_grid_jmax(:,1:(size_array(2) - 1)) = lon_in_range(:,2:size_array(2)) - (diff(lon_in_range,1,2)/2);
    lat_grid_jmax(:,1:(size_array(2) - 1)) = lat_in_range(:,2:size_array(2)) - (diff(lat_in_range,1,2)/2);
    lon_grid_vel_imin = NaN(size_array(1:2));
    lon_grid_vel_imin(2:size_array(1),:) = lon_in_range_vel(1:(size_array(1) - 1),:);
    lon_grid_vel_jmin = NaN(size_array(1:2));
    lon_grid_vel_jmin(:,2:size_array(2)) = lon_in_range_vel(:,1:(size_array(2) - 1));
    lon_grid_vel_min_iface = reshape(min([reshape(lon_grid_vel_jmin,[prod(size_array(1:2)) 1]) reshape(lon_in_range_vel,[prod(size_array(1:2)) 1])],[],2),size_array(1:2));
    lon_grid_vel_min_jface = reshape(min([reshape(lon_grid_vel_imin,[prod(size_array(1:2)) 1]) reshape(lon_in_range_vel,[prod(size_array(1:2)) 1])],[],2),size_array(1:2));
    lon_grid_vel_max_iface = reshape(max([reshape(lon_grid_vel_jmin,[prod(size_array(1:2)) 1]) reshape(lon_in_range_vel,[prod(size_array(1:2)) 1])],[],2),size_array(1:2));
    lon_grid_vel_max_jface = reshape(max([reshape(lon_grid_vel_imin,[prod(size_array(1:2)) 1]) reshape(lon_in_range_vel,[prod(size_array(1:2)) 1])],[],2),size_array(1:2));
    lat_grid_vel_jmin = NaN(size_array(1:2));
    lat_grid_vel_jmin(:,2:size_array(2)) = lat_in_range_vel(:,1:(size_array(2) - 1));
    
    % find interfaces that are along transect
    along_transect_ind = find(abs([sum(transect_mask_x,3) sum(transect_mask_y,3)]) > 0.5);
    
    % determine start (westernmost) face of transect
    curr_array = [reshape(lon_grid_vel_min_iface,[prod(size_array(1:2)) 1]); reshape(lon_grid_vel_min_jface,[prod(size_array(1:2)) 1])];
    lon_transect_min = curr_array(along_transect_ind);
    [~,start_transect_ind] = min(lon_transect_min);
    start_transect_ind = along_transect_ind(start_transect_ind);
    if start_transect_ind <= prod(size_array(1:2))
        start_lon = lon_grid_vel_jmin(start_transect_ind);
    else
        start_lon = lon_grid_vel_imin(mod(start_transect_ind - 1,prod(size_array(1:2))) + 1);
    end
    
    % recalculate to ensure consistency regardless of exact endpoint locations on land
    dist_from_start_bound = NaN([size_array(1:2) 2]);
    dist_from_start_bound(:,:,1) = 111.1*abs((cosd((0.5*lat_transect) + (0.5*lat_grid_imax)).*(lon_grid_imax - start_lon)) + (1i*(lat_grid_imax - lat_transect)));
    dist_from_start_bound(:,:,2) = 111.1*abs((cosd((0.5*lat_transect) + (0.5*lat_grid_jmax)).*(lon_grid_jmax - start_lon)) + (1i*(lat_grid_jmax - lat_transect)));
    dist_from_start_imax_jmin = 111.1*abs((cosd((0.5*lat_transect) + (0.5*lat_grid_vel_jmin)).*(lon_in_range_vel - start_lon)) + (1i*(lat_grid_vel_jmin - lat_transect)));
    dist_from_start_imin_jmax = 111.1*abs((cosd((0.5*lat_transect) + (0.5*lat_in_range_vel)).*(lon_grid_vel_imin - start_lon)) + (1i*(lat_in_range_vel - lat_transect)));
    % sort transect faces in order
    [dist_transect,sorted_along_transect_ind] = sort(dist_from_start_bound(along_transect_ind),'ascend');
    sorted_transect_ind = along_transect_ind(sorted_along_transect_ind);
    
    
    start_lon_vec(transect_ind) = start_lon;
    grid_dir_array(1:length(sorted_transect_ind),transect_ind) = ceil(sorted_transect_ind/prod(size_array(1:2)));
    dist_transect_array(1:length(sorted_transect_ind),transect_ind) = dist_transect;
    curr_array = [reshape(dist_from_start_imax_jmin,[prod(size_array(1:2)) 1]); reshape(dist_from_start_imin_jmax,[prod(size_array(1:2)) 1])];
    dist_from_start_transect_min_array(1:length(sorted_transect_ind),transect_ind) = curr_array(sorted_transect_ind);
    curr_array = [reshape(lon_grid_imax,[prod(size_array(1:2)) 1]); reshape(lon_grid_jmax,[prod(size_array(1:2)) 1])];
    lon_transect_array(1:length(sorted_transect_ind),transect_ind) = curr_array(sorted_transect_ind);
    curr_array = [reshape(lon_grid_vel_min_iface,[prod(size_array(1:2)) 1]); reshape(lon_grid_vel_min_jface,[prod(size_array(1:2)) 1])];
    lon_transect_min_array(1:length(sorted_transect_ind),transect_ind) = curr_array(sorted_transect_ind);
    curr_array = [reshape(lon_grid_vel_max_iface,[prod(size_array(1:2)) 1]); reshape(lon_grid_vel_max_jface,[prod(size_array(1:2)) 1])];
    lon_transect_max_array(1:length(sorted_transect_ind),transect_ind) = curr_array(sorted_transect_ind);
    curr_array = [reshape(hte,[prod(size_array(1:2)) 1]); reshape(htn,[prod(size_array(1:2)) 1])];
    dx_transect_array(1:length(sorted_transect_ind),transect_ind) = curr_array(sorted_transect_ind);
    curr_array = [reshape(angle_imax,[prod(size_array(1:2)) 1]); reshape(angle_jmax,[prod(size_array(1:2)) 1])];
    angle_transect_array(1:length(sorted_transect_ind),transect_ind) = curr_array(sorted_transect_ind);
    curr_array = [reshape(dzu,[prod(size_array(1:2)) size_array(3)]); reshape(dzv,[prod(size_array(1:2)) size_array(3)])];
    dz_transect_array(1:length(sorted_transect_ind),transect_ind,:) = reshape(curr_array(sorted_transect_ind,:),[length(sorted_transect_ind) 1 size_array(3)]);
    curr_array = [reshape(transect_mask_x,[prod(size_array(1:2)) size_array(3)]); reshape(transect_mask_y,[prod(size_array(1:2)) size_array(3)])];
    transect_mask_array(1:length(sorted_transect_ind),transect_ind,:) = reshape(curr_array(sorted_transect_ind,:),[length(sorted_transect_ind) 1 size_array(3)]);
    curr_array = [reshape(uvel_interp,[prod(size_array(1:2)) size_array(3:4)]); reshape(vvel_interp,[prod(size_array(1:2)) size_array(3:4)])];
    vel_cross_transect_array(1:length(sorted_transect_ind),transect_ind,:,:) = reshape(curr_array(sorted_transect_ind,:,:),[length(sorted_transect_ind) 1 size_array(3:4)]);
    curr_array = [reshape(temp_interp_x,[prod(size_array(1:2)) size_array(3:4)]); reshape(temp_interp_y,[prod(size_array(1:2)) size_array(3:4)])];
    temp_transect_array(1:length(sorted_transect_ind),transect_ind,:,:) = reshape(curr_array(sorted_transect_ind,:,:),[length(sorted_transect_ind) 1 size_array(3:4)]);
    curr_array = [reshape(uT_imax,[prod(size_array(1:2)) size_array(3:4)]); reshape(vT_jmax,[prod(size_array(1:2)) size_array(3:4)])];
    velT_cross_transect_array(1:length(sorted_transect_ind),transect_ind,:,:) = reshape(curr_array(sorted_transect_ind,:,:),[length(sorted_transect_ind) 1 size_array(3:4)]);
    
end

good_ind = find(isnan(lon_transect_array) == 0);
max_good_i_ind = max(mod(good_ind - 1,size(lon_transect_array,1)) + 1);
grid_dir_array = grid_dir_array(1:max_good_i_ind,:);
dist_from_start_transect_min_array = dist_from_start_transect_min_array(1:max_good_i_ind,:);
dist_transect_array = dist_transect_array(1:max_good_i_ind,:);
lon_transect_array = lon_transect_array(1:max_good_i_ind,:);
lon_transect_min_array = lon_transect_min_array(1:max_good_i_ind,:);
lon_transect_max_array = lon_transect_max_array(1:max_good_i_ind,:);
dx_transect_array = dx_transect_array(1:max_good_i_ind,:);
angle_transect_array = angle_transect_array(1:max_good_i_ind,:);
dz_transect_array = dz_transect_array(1:max_good_i_ind,:,:);
transect_mask_array = transect_mask_array(1:max_good_i_ind,:,:);
vel_cross_transect_array = vel_cross_transect_array(1:max_good_i_ind,:,:,:);
temp_transect_array = temp_transect_array(1:max_good_i_ind,:,:,:);
velT_cross_transect_array = velT_cross_transect_array(1:max_good_i_ind,:,:,:);

size_array = size(vel_cross_transect_array);

output_filename = [output_file_form,'.',input_dayfile,'.nc'];


nccreate(output_filename,'i_ind','Dimensions',{'i_ind',size_array(1)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'i_ind',(1:1:size_array(1))')
ncwriteatt(output_filename,'i_ind','long_name','grid cell number along transect')
ncwriteatt(output_filename,'i_ind','units','integer')
nccreate(output_filename,'lat_transect','Dimensions',{'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'lat_transect',lat_transects_vec)
ncwriteatt(output_filename,'lat_transect','long_name','latitude')
ncwriteatt(output_filename,'lat_transect','units','degrees N')
nccreate(output_filename,'z_t','Dimensions',{'z_t',size_array(3)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'z_t',depth_in_range)
ncwriteatt(output_filename,'z_t','long_name','depth from surface to midpoint of layer')
ncwriteatt(output_filename,'z_t','units','meters')
ncwriteatt(output_filename,'z_t','positive','down')
nccreate(output_filename,'z_w','Dimensions',{'z_w',size_array(3)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'z_w',depth_in_range_wvel_top)
ncwriteatt(output_filename,'z_w','long_name','depth from surface to top of layer')
ncwriteatt(output_filename,'z_w','units','meters')
ncwriteatt(output_filename,'z_w','positive','down')
nccreate(output_filename,'z_w_bot','Dimensions',{'z_w_bot',size_array(3)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'z_w_bot',depth_in_range_wvel)
ncwriteatt(output_filename,'z_w_bot','long_name','depth from surface to bottom of layer')
ncwriteatt(output_filename,'z_w_bot','units','meters')
ncwriteatt(output_filename,'z_w_bot','positive','down')
nccreate(output_filename,'time','Dimensions',{'time',Inf},'Datatype','double','Format','64bit')
ncwrite(output_filename,'time',time)
ncwriteatt(output_filename,'time','long_name','time')
ncwriteatt(output_filename,'time','units','days since 0000-01-01 00:00:00')
ncwriteatt(output_filename,'time','bounds','time_bound')
ncwriteatt(output_filename,'time','calendar','noleap')
nccreate(output_filename,'time_bound','Dimensions',{'d2',2,'time',Inf},'Datatype','double','Format','64bit')
ncwrite(output_filename,'time_bound',time_bound)
ncwriteatt(output_filename,'time_bound','long_name','boundaries for time-averaging interval')
ncwriteatt(output_filename,'time_bound','units','days since 0000-01-01 00:00:00')
nccreate(output_filename,'lon_start','Dimensions',{'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'lon_start',start_lon_vec)
ncwriteatt(output_filename,'lon_start','long_name','longitude at start of transect')
ncwriteatt(output_filename,'lon_start','units','degrees E')
ncwriteatt(output_filename,'lon_start','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'lon_start','missing_value',-9.96e+36)
nccreate(output_filename,'dist_from_start_transect_min','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'dist_from_start_transect_min',dist_from_start_transect_min_array)
ncwriteatt(output_filename,'dist_from_start_transect_min','long_name','Cartesian distance from start of transect to start of cell face')
ncwriteatt(output_filename,'dist_from_start_transect_min','units','kilometers')
ncwriteatt(output_filename,'dist_from_start_transect_min','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'dist_from_start_transect_min','missing_value',-9.96e+36)
nccreate(output_filename,'dist_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'dist_transect',dist_transect_array)
ncwriteatt(output_filename,'dist_transect','long_name','Cartesian distance from start of transect to midpoint of cell face')
ncwriteatt(output_filename,'dist_transect','units','kilometers')
ncwriteatt(output_filename,'dist_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'dist_transect','missing_value',-9.96e+36)
nccreate(output_filename,'grid_dir','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'grid_dir',grid_dir_array)
ncwriteatt(output_filename,'grid_dir','long_name','grid face indicator (1 = east face, 2 = north face)')
ncwriteatt(output_filename,'grid_dir','units','integer')
ncwriteatt(output_filename,'grid_dir','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'grid_dir','missing_value',-9.96e+36)
nccreate(output_filename,'lon_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'lon_transect',lon_transect_array)
ncwriteatt(output_filename,'lon_transect','long_name','longitude at center of grid face')
ncwriteatt(output_filename,'lon_transect','units','degrees E')
ncwriteatt(output_filename,'lon_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'lon_transect','missing_value',-9.96e+36)
nccreate(output_filename,'lon_transect_min','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'lon_transect_min',lon_transect_min_array)
ncwriteatt(output_filename,'lon_transect_min','long_name','longitude at start of grid face')
ncwriteatt(output_filename,'lon_transect_min','units','degrees E')
ncwriteatt(output_filename,'lon_transect_min','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'lon_transect_min','missing_value',-9.96e+36)
nccreate(output_filename,'lon_transect_max','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'lon_transect_max',lon_transect_max_array)
ncwriteatt(output_filename,'lon_transect_max','long_name','longitude at end of grid face')
ncwriteatt(output_filename,'lon_transect_max','units','degrees E')
ncwriteatt(output_filename,'lon_transect_max','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'lon_transect_max','missing_value',-9.96e+36)
nccreate(output_filename,'dx_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'dx_transect',dx_transect_array)
ncwriteatt(output_filename,'dx_transect','long_name','cell width along grid face')
ncwriteatt(output_filename,'dx_transect','units','meters')
ncwriteatt(output_filename,'dx_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'dx_transect','missing_value',-9.96e+36)
nccreate(output_filename,'angle_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'angle_transect',angle_transect_array)
ncwriteatt(output_filename,'angle_transect','long_name','angle grid makes with latitude line')
ncwriteatt(output_filename,'angle_transect','units','radians')
ncwriteatt(output_filename,'angle_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'angle_transect','missing_value',-9.96e+36)
nccreate(output_filename,'dz_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2),'z_t',size_array(3)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'dz_transect',dz_transect_array)
ncwriteatt(output_filename,'dz_transect','long_name','thickness of grid face')
ncwriteatt(output_filename,'dz_transect','units','meters')
ncwriteatt(output_filename,'dz_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'dz_transect','missing_value',-9.96e+36)
nccreate(output_filename,'transect_mask','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2),'z_t',size_array(3)},'Datatype','double','Format','64bit')
ncwrite(output_filename,'transect_mask',transect_mask_array)
ncwriteatt(output_filename,'transect_mask','long_name','direction of northward flow across grid face')
ncwriteatt(output_filename,'transect_mask','units','integer')
ncwriteatt(output_filename,'transect_mask','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'transect_mask','missing_value',-9.96e+36)
nccreate(output_filename,'vel_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2),'z_t',size_array(3),'time',Inf},'Datatype','double','Format','64bit')
ncwrite(output_filename,'vel_transect',vel_cross_transect_array)
ncwriteatt(output_filename,'vel_transect','long_name','velocity across grid face')
ncwriteatt(output_filename,'vel_transect','units','meters/second')
ncwriteatt(output_filename,'vel_transect','cell_methods','time: mean')
ncwriteatt(output_filename,'vel_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'vel_transect','missing_value',-9.96e+36)
nccreate(output_filename,'temp_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2),'z_t',size_array(3),'time',Inf},'Datatype','double','Format','64bit')
ncwrite(output_filename,'temp_transect',temp_transect_array)
ncwriteatt(output_filename,'temp_transect','long_name','potential temperature along grid face')
ncwriteatt(output_filename,'temp_transect','units','degrees Celsius')
ncwriteatt(output_filename,'temp_transect','cell_methods','time: mean')
ncwriteatt(output_filename,'temp_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'temp_transect','missing_value',-9.96e+36)
nccreate(output_filename,'velT_transect','Dimensions',{'i_ind',size_array(1),'lat_transect',size_array(2),'z_t',size_array(3),'time',Inf},'Datatype','double','Format','64bit')
ncwrite(output_filename,'velT_transect',velT_cross_transect_array)
ncwriteatt(output_filename,'velT_transect','long_name','flux of potential temperature across grid face')
ncwriteatt(output_filename,'velT_transect','units','(meters/second)*(degrees Celsius)')
ncwriteatt(output_filename,'velT_transect','cell_methods','time: mean')
ncwriteatt(output_filename,'velT_transect','_FillValue',-9.96e+36)
ncwriteatt(output_filename,'velT_transect','missing_value',-9.96e+36)