% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% Aim: This scripts is for deal GLEAM DATA to agree with CLM Output NC file format and coordinate.
%      The GLEAM is 0.25 x 0.25 deg, the CLM-Output is 0.5 x 0.5 deg.
% Date: 2021-04-10
% Author: Dayang Wang, Sun Yat-sen University
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 


% clear everything
clc; clear all

disp(['----------  Script Executes Start at : ',datestr(datetime),'---------- ']);
% read nc data and Var
varName = 'Ei';
varNameLong = 'Interception loss';
varUnits = 'mm.month-1';

InPath = 'H:\TempData\GLEAM_Monthly_v3.5a\monthly_update';
Infile = [varName,'_1980-2020_GLEAM_v3.5a_MO.nc'];
InPathFile = [InPath, '\', Infile];

varValue = ncread(InPathFile, varName);

% pack the 4 grids with 0.25 deg to 1 cell with 0.5 deg
DistrLat  = 2*ones(1, 360);
DistrLon  = 2*ones(1, 720);
DistrTim  = 1*ones(1, 492);
E_Cell = mat2cell(varValue, DistrLat, DistrLon, DistrTim);

% sum the 16 grids to 1 grid, except nan grids
varValueR05D = nan(size(E_Cell));
for iLat = 1:360
    for iLon = 1:720
        parfor iT = 1:492
            tempGrid = E_Cell{iLat, iLon, iT};
            tempGridAlign = tempGrid(:);
            % check the number of NaN
            NumNan = sum(isnan(tempGridAlign));
            if  NumNan < 2 
                ECellNonan = tempGridAlign(~isnan(tempGridAlign));
                ECellSum   = sum(ECellNonan)/(4-NumNan);   % the denominator is the (4-NumNan)
            else
                ECellSum = nan;
            end
            varValueR05D(iLat, iLon, iT) = ECellSum;
        end
    end
end

% shift the lon format from [-180, 180] to [0, 360] and shift fomat to CLM output
varValueR05D = circshift(varValueR05D, [0, 360, 0]);
for iT = 1:492
    tmp = varValueR05D(:,:,iT);
    tmp = flip(tmp);
    tmp = tmp';
    TransF(:,:,iT) = tmp;
end
varValueR05D = TransF;

% ==================================================================================
% make nc file
% ==================================================================================

outDir  = 'H:\TempData\GLEAM_Monthly_v3.5a\monthly_res05D';
outFile = [varName, '_1980-2020_GLEAM_v3.5a_MO_res05D.nc'];
outDirFile = [outDir, '\', outFile];
outid   = netcdf.create(outDirFile, 'CLASSIC_MODEL');

% define dimensions
latdimID  = netcdf.defDim(outid, 'lat',  360);
londimID  = netcdf.defDim(outid, 'lon',  720);
timedimID = netcdf.defDim(outid, 'time', 492);

% define variables
lat_id  = netcdf.defVar(outid, 'lat', 'NC_DOUBLE', latdimID);
lon_id  = netcdf.defVar(outid, 'lon', 'NC_DOUBLE', londimID);
time_id = netcdf.defVar(outid, 'time','NC_DOUBLE', timedimID);
E_id    = netcdf.defVar(outid, varName,'NC_DOUBLE', [londimID, latdimID, timedimID]);

netcdf.endDef(outid);

% time
timeVar  = [1:492]';

% recalculate the lat and lon gridcell centre point
latrR05D = [-89.75 : 0.5 : 89.75]';
lonrR05D = [0.25 : 0.5 : 359.75]';

% put variable into nc
netcdf.putVar(outid, lat_id, latrR05D);
netcdf.putVar(outid, lon_id, lonrR05D);
netcdf.putVar(outid, time_id, timeVar);
netcdf.putVar(outid, E_id, varValueR05D);

netcdf.reDef(outid);

% put Attribute into nc
netcdf.putAtt(outid, lat_id, 'standard_name', 'latitude');
netcdf.putAtt(outid, lat_id, 'units', 'degrees_north');

netcdf.putAtt(outid, lon_id, 'standard_name', 'longitude');
netcdf.putAtt(outid, lon_id, 'units', 'degrees_east');

netcdf.putAtt(outid, time_id, 'standard_name', 'time');
netcdf.putAtt(outid, time_id, 'units', 'months since 1980-01');

netcdf.putAtt(outid, E_id, 'long_name', varNameLong);
netcdf.putAtt(outid, E_id, 'units', varUnits);
netcdf.putAtt(outid, E_id, '_FillValue', -999.0);
netcdf.putAtt(outid, E_id, 'missing_value', -999.0);

% define global attribution
netcdf.putAtt(outid, netcdf.getConstant('NC_GLOBAL'), 'Frequency', 'mon');
netcdf.putAtt(outid, netcdf.getConstant('NC_GLOBAL'), 'Resolution', '0.5 deg');
netcdf.putAtt(outid, netcdf.getConstant('NC_GLOBAL'), 'Create_time', datestr(datetime));
netcdf.putAtt(outid, netcdf.getConstant('NC_GLOBAL'), 'E_mail', 'wangdy58@mail2.sysu.edu.cn');

netcdf.close(outid)

disp(['----------   Script Finishes at: ',datestr(datetime),'---------- ']);