clc
clear all 
close all

glob_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/';

glob_rst_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CNCLM40_9km.catchcnclm40_internal_rst.19810101_0000';
glob_tile_id = ncread(glob_rst_file,'TILE_ID');

glob_coord_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/rc_out/GEOSldas_CNCLM40_9km.ldas_tilecoord.bin';
[glob_tile_coord ] = read_tilecoord(glob_coord_file);

reg_coord_file = '/discover/nobackup/ylim5/EXP_ldas/Rst_v14_240517/SMAP_Nature_v9.1/output/SMAP_EASEv2_M09_GLOBAL/rc_out/SMAP_Nature_v9.1.ldas_tilecoord.bin';
[reg_tile_coord ] = read_tilecoord(reg_coord_file);

reg_index = zeros(size(glob_tile_coord.tile_id));

for i = 1:length(glob_tile_coord.tile_id)

    comp_index = (reg_tile_coord.tile_id==glob_tile_coord.tile_id(i));

    if sum(comp_index)~=0
       reg_index(i) = 1;
    end
end

reg_index = reg_index==1;

% read data

start_year = 1981;
stop_year = 2021; 

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];

mat_index = 0;

for y = start_year:stop_year

    y_str = num2str(y,'%04i');

    % define month length based on leap year status
    if sum(y==leapYears)>0
       monthLength = monthLength_ly;
    else
       monthLength = monthLength_nly;
    end

    for m = 1:12

        mat_index = mat_index + 1;

        m_str = num2str(m,'%02i');

        glob_read_path = [glob_path 'Y' y_str '/M' m_str '/'];
        glob_read_file = [glob_read_path 'GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        glob_scuns_data = ncread(glob_read_file,'SCUNS');
        reg_scuns_data = glob_scuns_data(reg_index);

        % write regional data to new file

        save_file = [glob_read_path 'GEOSldas_CNCLM40_9km.SCUNS.monthly.' y_str m_str '.nc4'];

        cmd = ['rm -rf ' save_file];
        status = system(cmd);

        nccreate(save_file,"SCUNS","Dimensions",{"tile",size(reg_scuns_data,1),"time",1}, ...
                   "FillValue",1e15);

        ncwrite(save_file,'SCUNS',reg_scuns_data);
        ncwriteatt(save_file,'SCUNS',"long_name","unstressed_stomatal_conductance");
        ncwriteatt(save_file,'SCUNS',"units","m s-1");

     end
end
