%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%  Init MData %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Returns 3 cell arrays: Molwc, Mflux, listMeso
%     and 1 vector:      szMolwc, szMflux

% MERL olwc data is stored in cell array Molwc{meso}(i,j)
% 
% Molwc Indices:
% meso - designates mesocosm (1-9)
% i - designates data point (size numData)
% j - designates field (1-7): 
%
% 1     Julian Day
% 2     O2 (mg O/L)
% 3     NH4 (mg N/L)
% 4     NO3 (mg N/L)
% 5     PO4 (mg P/L)
% 6     Si (mg Si/L)
% 7     Temp (deg C)

% Mflux Indices:
% meso - designates mesocosm (1-9)
% i - designates data point (size numData)
% j - designates field (1-7): 
%
% 1     Julian Day
% 2     SOD  (g  O2/m2-d)
% 3     JNH4 (mg  N/m2-d)
% 4     JNO3 (mg  N/m2-d)
% 5     JPO4 (mg  P/m2-d)
% 6     JSi  (mg Si/m2-d)
% 7     JMn  (mg Mn/m2-d)


readMERLraw('../SED_Input/Unformatted/MERLdataTabDelimited.txt')
% MERL data has been read in, and is initially
% stored in 'data(i,j)'
% 'i' designates a datapoint, and 'j' is as above

% strMeso{i} stores strings which identify the mesocosms that
% correspond to each 'i' in data(i,j)
strMeso = textdata(4:end,1);  % The first 3 lines are junk
strMeso = strtrim(strMeso); % remove whitespace
% list1: a cell array which lists each mesocosm name exactly once
list1 = {'C1' 'C5' 'C8' '1X' '2X' '4X' '8X' '16X' '32X'};
%list1 = unique(cMeso); % This line doesn't work: the strings are stored
%out of order.

% Get the set of line indices that correspond to each mesocosm
mesoInd = cell(9,1);
for meso = 1:9
    mesoInd{meso} = strcmp(strMeso,list1{meso});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Make MERL olwc array for SFM %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Store MERL olwc data in Molwc{meso} where the index 'meso'
% designates each mesocosm
Molwc = cell(9,1);
for meso = 1:9
    % Columns 1 and 8-13 are used for forcing
    Molwc{meso} = data(mesoInd{meso},[1,8:13]);
end

% ---------------------------------------------------
% Interpolate through bad data points.
% When bad data is at the edge of a data set,
% it cannot be interpolated and needs to be removed.
% ---------------------------------------------------
for meso = 1:9
    for i = 1:7
        % Copy the days for the field to interpolate
        newDays = Molwc{meso}(:,1);
        % Find the bad data points of this field
        badInd = isnan(Molwc{meso}(:,i));
        % Create a "day" and "field" vector with bad data removed
        unbadDays = Molwc{meso}(:,1);
        unbadDays(badInd) = [];
        unbadField = Molwc{meso}(:,i);
        unbadField(badInd) = [];

        newField = interp1(unbadDays,unbadField,newDays);

        Molwc{meso}(:,i) = newField;
        
        % Now trim any bad data from the endpoints (they're still Nans)
        badInd = isnan(Molwc{meso}(:,i));
        Molwc{meso}(badInd,:) = [];
    end
end


% Count the number of data points for each mesocosm.
szMolwc = zeros(9,1);
for meso = 1:9
    szMolwc(meso) = length(Molwc{meso}(:,1));
end

% --------------------------------------------------------
% ----------- Make MERL spinup array for SFM -------------
% --------------------------------------------------------
% Spinup comes from the 1st year of C1 tank
% Manually remove the only bad spinup datum,
% and set up the input so that the spinup ends on the first day.
MERLspin = [data(1:6,:);data(8:10,:)];
numSpin = 9; %Number of data points / year
% --------------------------------------------------------


% --------------------------------------------------------
% ------------ Make MERL flux array for SFM --------------
% --------------------------------------------------------

% Store MERL flux data in Mflux{meso} where the index 'meso'
% designates each mesocosm
Mflux = cell(9,1);
for meso = 1:9
    % Columns 1-7 are used for calculating the cost
    Mflux{meso} = data(mesoInd{meso},1:7);
end

% ---------------------------------------------------
% Interpolate through bad data points.
% When bad data is at the edge of a data set,
% it cannot be interpolated and needs to be removed.
% ---------------------------------------------------
for meso = 1:9
    for i = 1:7
        % Copy the days for the field to interpolate
        newDays = Mflux{meso}(:,1);
        % Find the bad data points of this field
        badInd = isnan(Mflux{meso}(:,i));
        % Create a "day" and "field" vector with bad data removed
        unbadDays = Mflux{meso}(:,1);
        unbadDays(badInd) = [];
        unbadField = Mflux{meso}(:,i);
        unbadField(badInd) = [];

        newField = interp1(unbadDays,unbadField,newDays);

        Mflux{meso}(:,i) = newField;
        
        % Now trim any bad data from the endpoints (they're still Nans)
        badInd = isnan(Mflux{meso}(:,i));
        Mflux{meso}(badInd,:) = [];
    end
end

baseDate = datenum('Jan 01 1981');
for meso = 1:9
  Mflux{meso}(:,1) = Mflux{meso}(:,1) + baseDate;
end

% Count the number of data points for each mesocosm.
szMflux = zeros(9,1);
for meso = 1:9
    szMflux(meso) = length(Mflux{meso}(:,1));
end

% Redefine listMeso for ease of FORTRAN use
listMeso = {'C01' 'C05' 'C08' '01X' '02X' '04X' '08X' '16X' '32X'};

fprintf(1,'Sub init_Molwc_Mflux done.\n')

% This is a helper script, so clear obselete variables to reduce clutter:
clear data strMeso list1 textdata mesoInd badInd i meso new* unbad*
