clc
clear all
close all

storm_name = 'Idai';
data_path = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' storm_name];
data_files = dir(data_path);

lat_lon_file = '/discover/nobackup/projects/smaptc/jkolassa/archive/Idai/Idai_SMAP_4D/diag/Y2019/M01/Idai_SMAP_4D.tavg1_2d_ocn_Nx.20190131_1030z.nc4';
lat_data = double(ncread(lat_lon_file,'lat'));
lon_data = double(ncread(lat_lon_file,'lon'));
frland_data = double(ncread(lat_lon_file,'FRLAND'));

month_vec = string(['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec']);


for f = 3:(length(data_files)-2)

    bt_file = [data_path '/' data_files(f).name];
    bt_data = readmatrix(bt_file);
    if size(bt_data,2)==8 
       bt_data(:,8) = [];
    end

    bt_data_new = [];
    datenum_write = datetime(bt_data(1,2),bt_data(1,3),bt_data(1,4),bt_data(1,5),0,0);

    max_traj = max(bt_data(:,1));

    for t = 1:max_traj
                
        traj_index = (bt_data(:,1)==t);
        bt_data_traj = bt_data(traj_index,:);
        bt_data_traj_new = bt_data_traj;

        tmp= (bt_data_traj(:,5)<0);
        if sum(tmp>0)
           bt_file
           t
        end

        datenum_start = datetime(bt_data_traj(1,2),bt_data_traj(1,3),bt_data_traj(1,4),bt_data_traj(1,5),0,0);

        for i = 2:size(bt_data_traj,1)
 
            datenum_new = datenum_start - (i-1)*3/24;
            datenum_str = char(datenum_new);

            bt_data_traj_new(i,2) = str2double(datenum_str(8:11));
            tmp = string(datenum_str(4:6));
            str_index = strcmp(month_vec,tmp);
            month_index = find(str_index==1);
            bt_data_traj_new(i,3) = month_index;
            bt_data_traj_new(i,4) = str2double(datenum_str(1:2));
            if length(datenum_str) > 11
               bt_data_traj_new(i,5) = str2double(datenum_str(13:14));
            end                      
            if (bt_data_traj_new(i,5)<0)
                datenum_str
            end
            bt_data_traj_new(i,8) = (i-1)*3;
        end % i
        bt_data_new = [bt_data_new; bt_data_traj_new];
    end % t

    land_mask_mat = zeros(size(frland_data));
    land_mask_mat_tr = zeros(size(frland_data,1),size(frland_data,2),max(bt_data_new(:,1)));

    for j = 1:size(bt_data_new,1)

        if(bt_data_new(j,8)<=24)
           lat_diff = sqrt((lat_data-bt_data_new(j,7)).^2);
           lat_index = find(lat_diff==min(lat_diff));
           lon_diff = sqrt((lon_data-bt_data_new(j,6)).^2);
           lon_index = find(lon_diff==min(lon_diff));

           if (frland_data(lon_index,lat_index)<0.5)
              land_mask_mat(lon_index,lat_index) = 1;
              land_mask_mat_tr(lon_index,lat_index,bt_data_new(j,1)) = 1;
           end
        end
    end % j

    f
    sum(land_mask_mat(:))
    ocean_mask_mat = land_mask_mat;

    tmp = char(datenum_start)
    tmp(12) = '_';

    save_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/Idai/land_mask/ocean_mask_2_24h_' tmp '.mat'];
    save(save_file,'ocean_mask_mat');
   
end % f

