%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Produce a MERLolwc*.inp file %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

spinupOnly = 0;

numSpinYrs = 25;

init_Molwc_Mflux %_Spline

% init_Molwc creates
% listMeso{meso} - stores a list of strings identifying the mesocosms
%
% Molwc{meso}(i,j) - stores all the MERL olwc data
% Indices:
% meso - designates mesocosm (1-9)
% i - designates data point (size szMolwc(meso))
% 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{meso}(i,j) - stores all the MERL flux data
% Indices:
% meso - designates mesocosm (1-9)
% i - designates data point (size szMflux(meso))
% 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)
%
% MERLspin(i,j) - stores MERL spinup data
% i - designates data point (size 9)
% j- designates field as Molwc

% ---------------------------------------------------
% For each mesocosm, open a MERLolwc*.inp forcing
% file, and print the MERL olwc data
% ---------------------------------------------------

for meso = 1:9
    
    % Open the input file
    if(spinupOnly == 1)
      filenameSTR = sprintf('../SED_Input/MERLolwc_spinOnly.inp')
    else
      filenameSTR = sprintf('../SED_Input/Spin_25yr/MERLolwc%s_spin%d.inp', listMeso{meso},numSpinYrs);
    end
    
    inpFile = fopen(filenameSTR,'w');

    % Initialize/Calculate some values for printing
    refDate = datenum('01-Jan-1981');
    % PON dep. fluxes by mesocosm,[mg N/m2-day]:
    %JPON = [47.56, 45, 45, 73.93, 64.05, 80.68, 72.63, 146.75, 108.65];
    % Dom's values: 
    JPON = [45, 45, 45, 70, 55, 90, 100, 115, 130];
    
    %Redfield: CNP  106:16:1 (moles)
    %               1273.06 : 224.16 : 30.97 (mass)
    JPOC = JPON(meso) * (1273.01)/(224.16);
    JPOP = JPON(meso) * (30.97)/(224.16);

    JSORP = 0;
    JSORSI = 7.2;
    sal = 35;
    pH = 0;
    refDate = datenum('01-Jan-1981'); %+ Molwc{meso}(1,1); % Days are in MATLAB days
    
    % Write the input file header
    fprintf(inpFile, '        day       jpoc       jpon       jpop      jsorp     jsorsi        sal');
    fprintf(inpFile, '        po4        nh4       no23         si         do       temp         pH\n');
    
    
    % --------- 25 year spinup -------------

    for year = numSpinYrs:-1:1  % Each iteration prints a year of forcing
    for spinIndex = 1:numSpin  % Each iteration prints a line
        fprintf(inpFile,'%11.0f',MERLspin(spinIndex,1) + refDate ...
            - round(365.25 * year)); % Model expects input days to start at 0.
        fprintf(inpFile,'%11.0f',JPOC(1)); %
        fprintf(inpFile,'%11.1f',JPON(1)); %
        fprintf(inpFile,'%11.2f',JPOP(1)); %
        fprintf(inpFile,'%11.1f',JSORP ); %
        fprintf(inpFile,'%11.1f',JSORSI); %
        fprintf(inpFile,'%11.0f',sal   ); %
        fprintf(inpFile,'%11.1e',MERLspin(spinIndex,11)/1000); % olwcPO4 - convert to mg/L
        fprintf(inpFile,'%11.1e',MERLspin(spinIndex,09)/1000); % olwcNH4
        fprintf(inpFile,'%11.1e',MERLspin(spinIndex,10)/1000); % olwcNO23
        fprintf(inpFile,'%11.1e',MERLspin(spinIndex,12)/1000); % olwcSi
        fprintf(inpFile,'%11.1e',MERLspin(spinIndex,08)); % olwcO2 - don't convert. Already mg/L
        fprintf(inpFile,'%11.1f',MERLspin(spinIndex,13)); % Tempereature
        fprintf(inpFile,'%11.1f',pH); % pH
        fprintf(inpFile,'\n');
    end  % Each iteration prints a dataIndex
    end  % Each iteration prints a year
    
    if(spinupOnly == 0)
        % ---------- MERL run -------------------------------------------------

        % Each iteration prints a line.  nuMolwc(meso) is the number of lines
        % to print per experiment
        for line = 1:szMolwc(meso)
            fprintf(inpFile,'%11.0f',Molwc{meso}(line,1)); % Model only cares about the difference between days
            fprintf(inpFile,'%11.0f',JPOC      ); %
            fprintf(inpFile,'%11.1f',JPON(meso)); %
            fprintf(inpFile,'%11.2f',JPOP      ); %
            fprintf(inpFile,'%11.1f',JSORP     ); %
            fprintf(inpFile,'%11.1f',JSORSI    ); %
            fprintf(inpFile,'%11.0f',sal       ); %
            fprintf(inpFile,'%11.4e',Molwc{meso}(line,5)/1000); % olwcPO4 - convert from ug/L to mg/L
            fprintf(inpFile,'%11.4e',Molwc{meso}(line,3)/1000); % olwcNH4
            fprintf(inpFile,'%11.4e',Molwc{meso}(line,4)/1000); % olwcNO23
            fprintf(inpFile,'%11.4e',Molwc{meso}(line,6)/1000); % olwcSi
            fprintf(inpFile,'%11.4e',Molwc{meso}(line,2)); % olwcO2 - don't convert. Already mg/L
            fprintf(inpFile,'%11.1f',Molwc{meso}(line,7)); % Tempereature
            fprintf(inpFile,'%11.1f',pH); % pH
            fprintf(inpFile,'\n');
        end  % Each iteration prints a line
    end % check to see if we're doing a "spinup only run"
    
end % Loop over each MERL experiment
fclose(inpFile);


% ---------------------------------------------------
% For each mesocosm, open a MERLflux*.inp data
% file, and print the MERL flux data
% ---------------------------------------------------


for meso = 1:9
    % Open the input file
    filenameSTR = sprintf('../SED_Input/Sediment_Fluxes/MERLflux%s.inp', listMeso{meso});
    inpFile = fopen(filenameSTR,'w');
    fprintf(inpFile, '        day        SOD       JNH4       JNO3       JPO4        JSi\n');
    % Each iteration prints a line.  szMflux(meso) is the number of lines
    % to print per experiment
    for line = 1:szMflux(meso)
        fprintf(inpFile,'%11.0f',Mflux{meso}(line,01)); % Days are in MATLAB days
        fprintf(inpFile,'%11.4f',Mflux{meso}(line,02)); % measured SOD
        fprintf(inpFile,'%11.4f',Mflux{meso}(line,03)); % measured JNH4
        fprintf(inpFile,'%11.4f',Mflux{meso}(line,04)); % measured JNO3
        fprintf(inpFile,'%11.4f',Mflux{meso}(line,05)); % measured JPO4
        fprintf(inpFile,'%11.4f',Mflux{meso}(line,06)); % measured JSi
        fprintf(inpFile,'\n');
    end  % Each iteration prints a line
end % Loop over each MERL experiment
fclose(inpFile);

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

% To reduce clutter, clear the obselete variables.
%clear J* refDate filenameSTR inpFile line meso pH sal

