clc
clear all
close all

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

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;
y_index = 0;

for y = start_year:stop_year

    y_index = y_index + 1;
    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

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

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

        read_file = [read_path_base 'Y' y_str '/M' m_str '/GEOSldas_CNCLM40_9km.SCSAT.monthly.' y_str m_str '.nc4'];

        rcsat_mon = ncread(read_file,'SCSAT');

        if (y_index==1)
           rcsat_avg(:,m) = rcsat_mon;
        else
           rcsat_avg(:,m) =  ((y_index-1)/y_index)*rcsat_avg(:,m) + (1/y_index)*rcsat_mon;
        end

    end
end

mat_index = 0;
y_index = 0;

for y = start_year:stop_year

    y_index = y_index + 1;
    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

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

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

        read_file = [read_path_base 'Y' y_str '/M' m_str '/GEOSldas_CNCLM40_9km.SCSAT.monthly.' y_str m_str '.nc4'];

        rcsat_mon = ncread(read_file,'SCSAT');

        rcsat_scaled = rcsat_mon./rcsat_avg(:,m);

        % write scaled data to new file

        save_file = [read_path_base 'Y' y_str '/M' m_str '/GEOSldas_CNCLM40_9km.SCSAT_scaled.monthly.' y_str m_str '.nc4'];

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

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

        ncwrite(save_file,'SCSAT_sc',rcsat_scaled);
        ncwriteatt(save_file,'SCSAT_sc',"long_name","saturated_stomatal_conductance_scaled_by_clim");
        ncwriteatt(save_file,'SCSAT_sc',"units","-");

     end
end
