%------------------------------------------%
%   Process forcing-file for JSB-Offline   %
%------------------------------------------%

clc;                % Refresh Command Window
clear;              % clear Workspace
close all;          % close all Figures
format long g       % 14 decimal places
addpath(genpath('../../../Functions')) % load Functions

m = 60; hr = 60*m; d = 24*hr; yr = 365*d;


%-------------------%
%   User settings   %
%-------------------%

plotSwitch = 1;
writeSwitch = true;

startYear  = 1990;
finalYear  = 2003;
qairExp    = true;


randYear = {'y1' 'y2' 'y4' 'y1' 'y2' 'y1' 'y3' 'y2' 'y2' 'y1' 'y3' 'y1' 'y1' 'y2'};   

% minimal wind speed (wind speed ought not drop on zero in numeric models)
eps_wind = 0.01;

T0 = 273.15;


%--------------------%
%   Read trop-file   %
%--------------------%
tropInputFolder = './';
tropInputFile   = 'tropFluxnet.mat';

load([tropInputFolder tropInputFile]);

tropVarName = {'TA_F' 'P_F'   'WS_F' 'SW_IN_F' 'LW_IN_F' 'Q_AIR' 'PA_F'};
tropVarNew  = {'Tair' 'Rainf' 'Wind' 'SWdown'  'LWdown'  'Qair'  'PSurf'};
tropLinInt  = [1 0 1 0 0 1 1];
tropVarNum  = length(tropVarName);



% convert time vector to Matlab time format
dtOld    = 30*m;
dtJsbOld = 3*hr;
dt       = 450; 
dtd      = dt/d;
tropLen  = length(trop.time);

% timeVec   = datevec(trop.time);
% jsbTime   = timeVec(:,1)*1e4 + timeVec(:,2)*1e2   + timeVec(:,3) ...
%           + timeVec(:,4)/24  + timeVec(:,5)/24/60 + timeVec(:,6)/24/60/60;
      
tropYear  = year(trop.time);

% for iVar = 1:diceVarNum
%     x = run.var.(diceVarName{iVar})(1:end-1);
%     x = nanmean(reshape(x,18,[]))';
%     run.var.(diceVarName{iVar}) = x;
% 
% end



%% Detecting and filling gaps
maxlingapfill=2*hr;    % Longest time period for linear gap filling in s
gapfilled=zeros(1,tropLen);
for iVar = 1:tropVarNum
    
    tropVar = trop.(tropVarName{iVar});
    index   = find(isnan(tropVar));
    nIdx    = length(index);
    
    if (nIdx>1)  % Only if there are gaps
        % Detecting gaps
        gapIndex = find(index(2:nIdx)-index(1:nIdx-1)>1);
        gapStart = index(gapIndex+1);
        gapFinal = index(gapIndex);
        gapStart = [index(1) gapStart'];
        gapFinal = [gapFinal' index(nIdx)];
        gapDurat = gapFinal-gapStart+1;
        
        for ng=1:length(gapStart)
            gapfilled(gapStart(ng):gapFinal(ng)) = 1;
            if ( gapStart(ng) == 1 )
                tropVar(gapStart(ng):gapFinal(ng)) = tropVar(gapFinal(ng)+1);
            elseif (gapFinal(ng)==tropLen)
                tropVar(gapStart(ng):gapFinal(ng)) = tropVar(gapStart(ng)-1);
            else
                weight = (1:gapDurat(ng))/(gapDurat(ng)+1);
                tropVar(gapStart(ng):gapFinal(ng)) ...
                    = tropVar(gapStart(ng)-1).*(1-weight) ...
                    + tropVar(gapFinal(ng)+1).*weight;
            end
            if ( gapDurat(ng)*dt > maxlingapfill )
                if ( gapDurat(ng)*dt >= 10*d )
                    shiftstep = yr;
                else
                    shiftstep = d;
                end
                tshift = fix(gapDurat(ng)*dt/(shiftstep)+1)*shiftstep/dt;
                if ( gapStart(ng) - tshift > 0 )
                    template = tropVar((gapStart(ng):gapFinal(ng)) - tshift);
                    if ( tropLinInt(iVar) )
                        tropVar(gapStart(ng):gapFinal(ng))...
                           = tropVar(gapStart(ng):gapFinal(ng)) ...
                           + detrend(template);
                    else
                        tropVar(gapStart(ng):gapFinal(ng)) = template;
                    end
                end
            end
        end
    end
    tropVar = linInt(dtOld,dt,tropVar);
    trop.(tropVarNew{iVar}) = tropVar;
    tropNew.(tropVarNew{iVar}) = tropVar;
       
%        figure(iVar)
%        hold on
%        plot(currVar)
%        plot(trop.(tropVarName{iVar}))
%        title(tropVarNew{iVar})
end

startNum  = datenum('01-Sep-2001'); 
finalNum  = datenum('31-Oct-2003 23:52:30'); 
nextNum   = finalNum + dtd;
dayNum    = nextNum-startNum;
ndt       = dayNum/dtd;
trop.time = linspace(startNum,finalNum,ndt)';
trop.year = year(trop.time);

% Remove calmes
trop.Wind(trop.Wind<eps_wind) = eps_wind;
trop.Tair = trop.Tair - T0;


tropNew.Wind(tropNew.Wind<eps_wind) = eps_wind;
tropNew.Tair = tropNew.Tair - T0;

tropNew.time = trop.time;

%trop = tropNew;
%save('tropFluxProcessed.mat','trop')


% Create two complete years for periodic forcing
finalPer   = datenum('31-Aug-2003 23:52:30');
[~,idxPer] = min(abs(trop.time-finalPer));

for iVar = 1:tropVarNum
    tropVar   = trop.(tropVarNew{iVar});
    tropVar   = tropVar(1:idxPer);
    tropTime  = trop.time(1:idxPer);
    tropOld   = tropVar;
    trop2003  = tropVar(year(tropTime)==2003);
    tropRest  = tropVar(year(tropTime)~=2003);
    tropVar   = [trop2003;tropRest];
    trop.y1.(tropVarNew{iVar}) = tropVar(1:end/2);
    trop.y2.(tropVarNew{iVar}) = tropVar(end/2+1:end);
    % leap year
    trop.y3.(tropVarNew{iVar}) = [tropVar(1:1/dtd);tropVar(1:end/2)];
    trop.y4.(tropVarNew{iVar}) = [tropVar(end/2+1:end/2+1/dtd);tropVar(end/2+1:end)];
    
    

end




allYear   = startYear:finalYear;
numYear   = length(allYear);
for iYear = 1:numYear
    
    currYear = allYear(iYear);
    
    
    %-------------------%
    %   Read jsb-file   %
    %-------------------%
    jsbInputFolder  = '~/Desktop/Work/JSB_OFFLINE/Forcing/Input/';
    jsbInputFile    = ['climate_wfdei-gpcc_164-50_3h_' num2str(currYear) '.nc'];
    
    [jsb.glo, jsb.dim, jsb.var] = readNetcdf([jsbInputFolder jsbInputFile]);
   
    jsbVarName   = jsb.var.Name(jsb.dim.Num+1:end);
    jsbVarNum    = length(jsbVarName);
    
    startNum   = datenum(num2str(currYear),'YYYY');
    nextNum    = datenum(num2str(currYear+1),'YYYY');
    finalNum   = nextNum - dtd;
    dayNum     = nextNum-startNum;
    ndt        = dayNum/dtd;
    time    = linspace(startNum,finalNum,ndt)';
    timeVec    = datevec(time);
    
    jsbTime   = timeVec(:,1)*1e4 + timeVec(:,2)*1e2   + timeVec(:,3) ...
          + timeVec(:,4)/24  + timeVec(:,5)/24/60 + timeVec(:,6)/24/60/60;
      
    jsb.var.time = jsbTime;
    jsbTimeLen   = length(jsb.var.time);
    jsb.dim.time = jsbTimeLen;
    jsb.dim.Len(strcmp(jsb.dim.Name,'time')) = jsbTimeLen;
    
    
    for iVar = 1:jsbVarNum
        %tropVar  = trop.(tropVarName{iVar})';
        %tropVar  = tropVar(tropYear==currYear);

        jsbVar = jsb.var.(jsbVarName{iVar});
        
        
        jsb.var.(jsbVarName{iVar})(1,1,1:jsbTimeLen) = trop.(randYear{iYear}).(tropVarNew{iVar});
        
        jsbVarOld = squeeze(jsb.var.(jsbVarName{iVar}));
        tropVar = trop.(tropVarNew{iVar})(trop.year==currYear);
        
        
        
        if currYear == 2001
            
            startNum = trop.time(1);
            nextNum  = datenum(num2str(currYear+1),'YYYY');
            finalNum = nextNum - dtd;
            
            [~,startJsb]  = min(abs(time - startNum));
            [~,finalJsb]  = min(abs(time - finalNum));
            
            if(qairExp && iVar == 6)
                n       = 28; % number of days
                y1      = tropVar(1:n/dtd)'-[1:n/dtd]*0.00001;
                y2      = tropVar(n/dtd+1:2*n/dtd)'-[n/dtd:-1:1]*0.00001;
                tropVar = [y1';y2';tropVar(2*n/dtd+1:end)];
                tropVar(tropVar<0) = 0;
                jsb.var.(jsbVarName{iVar})(1,1,startJsb:finalJsb) = tropVar';
            end
            jsb.var.(jsbVarName{iVar})(1,1,startJsb:finalJsb) = tropVar';
     
               
            
        elseif currYear == 2002
            
            jsb.var.(jsbVarName{iVar})(1,1,:) = tropVar';
            
        elseif currYear == 2003
            
            startNum = datenum(num2str(currYear),'YYYY');
            finalNum = trop.time(end);
            
            [~,startJsb]  = min(abs(time - startNum));
            [~,finalJsb]  = min(abs(time - finalNum));
            
            % Remove last 3 Nan's from Interpolation 
            tropVar(end-2:end) = tropVar(end-3);
            jsb.var.(jsbVarName{iVar})(1,1,startJsb:finalJsb) = tropVar';
            
        end
        
        
        
        
        jsbVarNew = squeeze(jsb.var.(jsbVarName{iVar}));
        
        if currYear >= 2001
            if(plotSwitch)
                figure
                hold on
                plot(jsbVarOld)
                plot(jsbVarNew)
                title([jsbVarName{iVar} '(' num2str(iYear) ')'])
            end
        end
        
        
        
        
        
    end

    %------------------------%
    %   Write new jsb-file   % 
    %------------------------%
    outputFolder    = './Output/';
    outputFile      = ['trop_periodforc_450sec_qair0_' num2str(currYear) '.nc'];
    
    if(writeSwitch)
        writeNetcdf([outputFolder outputFile], jsb.glo, jsb.dim, jsb.var)
    end

    

end









