% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%   This scripts is for estimate the ratio of GLEAM_Version3.5a/CLM_Output.
%   2021-02-26 by Dayon 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% clear everything
clc
clear all

yearStart = 1986;
yearEnd   = 1999;
nYear     = yearEnd - yearStart + 1;

monthStart = 1;
monthEnd   = 12;
nMonths    = monthEnd - monthStart + 1;

if monthStart < 10
    monthStartStr = ['0', num2str(monthStart)];
else
    monthStartStr = num2str(monthStart);
end
if monthEnd < 10
    monthEndStr = ['0', num2str(monthEnd)];
else
    monthEndStr = num2str(monthEnd);
end

startEndDate = ['yr_', num2str(yearStart),'_',num2str(yearEnd),'_mo_',monthStartStr,'_',monthEndStr];

ET_ratio_min = 0.01;
ET_ratio_max = 100;

res = 0.25;

% longtude and latitude
lon  = 0.5:1:359.5;
lat  = -89.5:1:89.5;
nLon = length(lon);
nLat = length(lat);

% ==================================================================================
% extract CLM5.0 Output Data 
% ==================================================================================
inDir = '/home/dayon/Downloads/GLEAM/ModelOutput';
temp2 = zeros(nLon, nLat, nMonths);
monthDays = [31 28 31 30 31 30 31 31 30 31 30 31];
for iMonth = monthStart : monthEnd
    if iMonth < 10
        iMonthStr = ['0', num2str(iMonth)];
    else
        iMonthStr = num2str(iMonth);
    end

    ET = zeros(nLon, nLat, nYear);
    for iYear = yearStart : yearEnd
        iYearStr  = num2str(iYear);  
        iYearInd  = iYear - yearStart + 1;
        inFile    = ['CLM50SpGs_res1D_CTRL.clm2.h0.',iYearStr,'-',iMonthStr,'.nc'];
        inDirFile = [inDir, '/', inFile];
        temp1     = ncread(inDirFile, 'QVEGE') + ncread(inDirFile, 'QVEGT') ...
                    + ncread(inDirFile, 'QSOIL');
        temp1(temp1 > 10000) = nan;
        % convert units from mm/s to mm/month
        ET(:, :, iYearInd) = monthDays(iMonth) * 3600 * 24 * temp1; 
    end
    temp2(:,:,iMonth) = reshape(mean(ET, 3), nLon, nLat);
end
ET_Model_MonthMean = temp2;

clear temp1 temp2
ET_Model_MonthMean_max = max(max(max(ET_Model_MonthMean)));
ET_Model_MonthMean_min = min(min(min(ET_Model_MonthMean)));

% ==================================================================================
% extract GLEAM Observation
% ==================================================================================
inDir  = ['/home/dayon/Downloads/GLEAM/res1D_GLEAM'];
inFile = ['E_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
inDirFile = [inDir, '/' , inFile];

% choose the E_Obs Type
ET_Option = 2;
switch ET_Option
case 1
    disp('You choose the Actual evaporation (E) as the Observation ET');
    ET_Obs = ncread(inDirFile, 'E');
case 2
    inFileEt = ['Et_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
    inDirFileEt = [inDir, '/' , inFileEt];
    inFileEi = ['Ei_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
    inDirFileEi = [inDir, '/' , inFileEi];
    inFileEb = ['Eb_1980-2020_GLEAM_v3.5a_MO_res1D.nc'];
    inDirFileEb = [inDir, '/' , inFileEb];
    disp('You choose the sum of Transpiration, Interception loss and Bare-soil evaporation (Et + Ei + Eb) as the Observation ET');
    ET_Obs = ncread(inDirFileEt, 'Et') + ncread(inDirFileEi, 'Ei') + ncread(inDirFileEb, 'Eb');
end


ET_Obs_MonthMean = zeros(nLon, nLat, nMonths);
for iMonth = monthStart : monthEnd
    ipMonStart = (yearStart - 1980)*12 + iMonth;
    ipMonEnd   = (yearEnd   - 1980)*12 + iMonth;
    ipArray    = [ipMonStart : 12 : ipMonEnd];
    temp3      = ET_Obs(:, :, ipArray);
    ET_Obs_MonthMean(:, :, iMonth) = mean(temp3, 3);
end
ET_Obs_MonthMean_min = min(min(min(ET_Obs_MonthMean)));
ET_Obs_MonthMean_max = max(max(max(ET_Obs_MonthMean)));

% ==================================================================================
% Ratio
% ==================================================================================
ET_ratio = ET_Obs_MonthMean./ET_Model_MonthMean;

ET_ratio(isnan(ET_ratio)) = 1.0;

ET_ratio(ET_ratio <= 0) = 1.0;

ET_ratio(ET_ratio >= ET_ratio_max) = ET_ratio_max;

ET_ratio(ET_ratio <= ET_ratio_min) = ET_ratio_min;

ET_ratio_max = max(max(max(ET_ratio)));
ET_ratio_min = min(min(min(ET_ratio)));

% ==================================================================================
% Make netcdf
% ==================================================================================
outDir  = '/home/dayon/Downloads/GLEAM/ET_ratio';

switch ET_Option
case 1 
    outFile = ['Global_ETactual_ratio_GLEAM_', startEndDate, '.nc'];
case 2
    outFile = ['Global_ETthree_ratio_GLEAM_', startEndDate, '.nc'];
end
outDirFile = [outDir, '/', outFile];
outid   = netcdf.create(outDirFile, 'CLOBBER');

% define dimension
latdimID = netcdf.defDim(outid, 'lat', nLat);
londimID = netcdf.defDim(outid, 'lon', nLon);
timedimID = netcdf.defDim(outid, 'time', nMonths);

% define variable
latid  = netcdf.defVar(outid, 'lat', 'double', latdimID);
lonid  = netcdf.defVar(outid, 'lon', 'double', londimID);
timeid = netcdf.defVar(outid, 'time', 'double', timedimID);
ETid   = netcdf.defVar(outid, 'ET_ratio', 'double', [londimID, latdimID, timedimID]);
QVEGEid = netcdf.defVar(outid, 'QVEGE_ratio', 'double', [londimID, latdimID, timedimID]); 
QVEGTid = netcdf.defVar(outid, 'QVEGT_ratio', 'double', [londimID, latdimID, timedimID]); 
QSOILid = netcdf.defVar(outid, 'QSOIL_ratio', 'double', [londimID, latdimID, timedimID]); 
netcdf.endDef(outid);

netcdf.putVar(outid, latid, lat);
netcdf.putVar(outid, lonid, lon);
timeVar = 0:1:11;
netcdf.putVar(outid, timeid, timeVar);
netcdf.putVar(outid, ETid, ET_ratio);
netcdf.putVar(outid, QVEGEid, ET_ratio);
netcdf.putVar(outid, QVEGTid, ET_ratio);
netcdf.putVar(outid, QSOILid, ET_ratio);
netcdf.reDef(outid);

netcdf.putAtt(outid,latid,'long_name','latitude');
netcdf.putAtt(outid,latid,'units','degrees_north');
netcdf.putAtt(outid,latid,'standard_name','latitude');

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

netcdf.putAtt(outid,timeid,'long_name','month');
netcdf.putAtt(outid,timeid,'units','months since 1980-01-01 00:00:00');
netcdf.putAtt(outid,timeid,'standard_name','month');

undef = -9999.9;
netcdf.putAtt(outid,ETid,'long_name','ET_ratio');
netcdf.putAtt(outid,ETid,'units','unitless');
netcdf.putAtt(outid,ETid,'_FillValue',undef);

netcdf.putAtt(outid,QVEGEid,'long_name','QVEGE_ratio');
netcdf.putAtt(outid,QVEGEid,'units','unitless');
netcdf.putAtt(outid,QVEGEid,'_FillValue',undef);

netcdf.putAtt(outid,QVEGTid,'long_name','QVEGT_ratio');
netcdf.putAtt(outid,QVEGTid,'units','unitless');
netcdf.putAtt(outid,QVEGTid,'_FillValue',undef);

netcdf.putAtt(outid,QSOILid,'long_name','QSOIL_ratio');
netcdf.putAtt(outid,QSOILid,'units','unitless');
netcdf.putAtt(outid,QSOILid,'_FillValue',undef);

netcdf.close(outid);