%------------------------------------------%
%   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
physicalConstants



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

plotSwitch  = false;
writeSwitch = true;

startYear  = 1996;%1986;
finalYear  = 2014;




randYear = [1997 1998 2004 2009 2002 2013 1996 2001 2006 2014];  



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


%--------------------%
%   Read Thar-file   %
%--------------------%
tharInputFolder = './';
tharInputFile   = 'tharFluxnet.mat';

load([tharInputFolder tharInputFile]);

tharVarName = {'TA_F' 'P_F'   'WS_F' 'SW_IN_F' 'LW_IN_F' 'Q_AIR' 'PA_F'};
tharVarNew  = {'Tair' 'Rainf' 'Wind' 'SWdown'  'LWdown'  'Qair'  'PSurf'};
tharLinInt  = [1 0 1 0 0 1 1];
tharVarNum  = length(tharVarName);



% convert time vector to Matlab time format
dtOld    = 30*m;
dt       = 450; 
dtd      = dt/d;
tharLen  = length(thar.time);

% timeVec   = datevec(thar.time);
% jsbTime   = timeVec(:,1)*1e4 + timeVec(:,2)*1e2   + timeVec(:,3) ...
%           + timeVec(:,4)/24  + timeVec(:,5)/24/60 + timeVec(:,6)/24/60/60;
      
tharYear  = year(thar.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,tharLen);
for iVar = 1:tharVarNum
    
    tharVar = thar.(tharVarName{iVar});
    index   = find(isnan(tharVar));
    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 )
                tharVar(gapStart(ng):gapFinal(ng)) = tharVar(gapFinal(ng)+1);
            elseif (gapFinal(ng)==tharLen)
                tharVar(gapStart(ng):gapFinal(ng)) = tharVar(gapStart(ng)-1);
            else
                weight = (1:gapDurat(ng))/(gapDurat(ng)+1);
                tharVar(gapStart(ng):gapFinal(ng)) ...
                    = tharVar(gapStart(ng)-1).*(1-weight) ...
                    + tharVar(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 = tharVar((gapStart(ng):gapFinal(ng)) - tshift);
                    if ( tharLinInt(iVar) )
                        tharVar(gapStart(ng):gapFinal(ng))...
                           = tharVar(gapStart(ng):gapFinal(ng)) ...
                           + detrend(template);
                    else
                        tharVar(gapStart(ng):gapFinal(ng)) = template;
                    end
                end
            end
        end
    end
    tharVar = linInt(dtOld,dt,tharVar);
    tharVar(end-2:end) = tharVar(end-3);
    thar.(tharVarNew{iVar}) = tharVar;
    tharNew.(tharVarNew{iVar}) = tharVar;
       
end

startNum  = datenum('01-Jan-1996'); 
finalNum  = datenum('31-Dec-2014 23:52:30'); 
nextNum   = finalNum + dtd;
dayNum    = nextNum-startNum;
ndt       = dayNum/dtd;
thar.time = linspace(startNum,finalNum,ndt)';
thar.year = year(thar.time);

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


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

tharNew.time = thar.time;










allYear   = startYear:finalYear;
numYear   = length(allYear);
for iYear = 1:numYear
    currYear = allYear(iYear);


   
    %-------------------%
    %   Read jsb-file   %
    %-------------------%
    jsbInputFolder  = './';
    jsbInputFile    = 'climate_wfdei-gpcc_141-28_3h_1990.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
        %tharVar  = thar.(tharVarName{iVar})';
        %tharVar  = tharVar(tharYear==currYear);

        jsbVar = jsb.var.(jsbVarName{iVar});
        
        currData = thar.(tharVarNew{iVar});
        if(currYear >= 1996)
            actualYear = currYear;
        else
            actualYear = randYear(iYear); 
        end
        
        currData = currData(year(thar.time) == actualYear);
        jsbVarOld = squeeze(jsb.var.(jsbVarName{iVar}));
        jsb.var.(jsbVarName{iVar})(1,1,1:jsbTimeLen) = currData;
      
      
    
  
        jsbVarNew = squeeze(jsb.var.(jsbVarName{iVar}));
        
        if currYear >= 1996
            if(plotSwitch)
                figure
                hold on
                plot(jsbVarNew)
                title([jsbVarName{iVar} '(' num2str(iYear) ')'])
            end
        end
        
        
        
        
        
    end

    %------------------------%
    %   Write new jsb-file   % 
    %------------------------%
    outputFolder    = './Output/';
    outputFile      = ['thar_periodforc_450sec_' num2str(currYear) '.nc'];
    
    if(writeSwitch)
        writeNetcdf([outputFolder outputFile], jsb.glo, jsb.dim, jsb.var)
    end
    if(currYear < 1996)
        disp(['Current year of initialization: ' num2str(randYear(iYear))]);
    end
    

end









