% J16b model
% Horsepool Monitoring Station
% Uinta Basin, Eastern Utah
% 29 Jan 2016

%{
     As noted in comments below, for some input parameters, we adopted values
     assigned by Edwards, et al., 'High Winter Ozone Pollution from Carbonyl
     Photolysis in an Oil and Gas Basin,' Nature, 514:351-354 (2014).
%}

clear

%% OBSERVATIONS

%{
     Most data are from 29 Jan 2016 at Horsepool.
     A variable called BLheight (Boundary Layer Height) is not being used in this
     version.  The computation runs OK even if I set it to zero.
     Other variables, M and kOH, are also not being used.  No harm comes when they
     are not assigned values.

%}

load Inputs/zero.csv Inputs.BLheight;
Inputs.BLheight = zero;

% The following input data are loaded from csv files.
load Inputs/J16bCO.csv   Inputs.CO;      %  CO concentrations, ppb
load Inputs/J16bNO.csv   Inputs.NO;      %  NO concentrations, ppb
load Inputs/J16bNO2.csv  Inputs.NO2;     %  NO2 concentrations, ppb
load Inputs/J16bP.csv    Inputs.P;       %  pressure, mbar
load Inputs/J16b_ATB_RH.csv   Inputs.RH;      %  Relative humidity, %
load Inputs/J16b_ATB_T.csv    Inputs.T;       %  Temperature, Kelvin
load Inputs/J16bkdil.csv Inputs.kdil;    %  Dilution constant, hz
load Inputs/J16bTime.csv Inputs.Time;    %  Time of day, hr
Inputs.CO   = J16bCO;
Inputs.NO   = J16bNO;
Inputs.NO2  = J16bNO2;
Inputs.P    = J16bP;
Inputs.RH   = J16b_ATB_RH;
Inputs.T    = J16b_ATB_T;
Inputs.kdil = J16bkdil;
Inputs.Time = J16bTime;

%% METEOROLOGY

%{
   P, T and RH are input above.  SZA (solar zenith angle) is calculated
   by the F0AM routine sun_position from time, date, and location.
   kdil is a physical loss constant controlling dilution of all compounds;
   it has been set to 1.8E-5/sec in daylight hours,
   and 2.0E-6/sec in nighttime hours.
%}

o = ones(size(Inputs.T));
time.year           = 2016*o;
time.month          = 01*o;
time.day            = 29*o;
time.hour           = Inputs.Time;
time.min            = 0*o;
time.sec            = 0*o;
time.UTC            = -7;
location.latitude   = 40.1434;
location.longitude  = -109.4689;
location.altitude   = 1569;
sun = sun_position(time,location); %fields zenith and azimuth

Met = {...
%   names       %values
    'P'          Inputs.P;        % Pressure, mbar
    'T'          Inputs.T;        % Temperature, K
    'RH'         Inputs.RH;       % Relative Humidity, %
    'SZA'        sun.zenith;      % Solar zenith angle, degrees
    'albedo'     0.70;             % Albedo, fraction
    'O3col'      275;             % Ozone column, dobsons units from Giovanni NASA website
    'kdil'       Inputs.kdil;     % Dilution constant, hz, adapted from Edwards et al.
    'jcorr'      0.5;             % Optimizes comparison b/w model and observed NO/NO2
    'BLH'        Inputs.BLheight; % Set to 0, so apparently not used by this version.
    };


%% CHEMICAL CONCENTRATIONS

%{
  Concentrations of NO and NO2 are initialized from observational data entered above.
  Concentrations of organics and aldehydes are initialized below using speciation data
  reported in 2021.  [Lyman et al. High ethylene and propylene in an area dominated by
  oil production, Atmosphere, 12, 1 (2021)]
  Species with HoldMe = 1 will be held constant throughout each step.
  Species with HoldMe = 0 are only initialized at the start of the run, because
       ModelOptions.LinkSteps = 1 (see below).
%}

% Multipliers for input concentrations

org_scalar =     0.65
ald_scalar =     0.65
nox_scalar =     1.00

InitConc = {...
    % names           conc(ppb)              HoldMe

    %Inorganics
    'O3'                  56                   0;
    'CO'               Inputs.CO               1;

    %NOy
    'NO'                Inputs.NO*nox_scalar   0;
    'NO2'               Inputs.NO2*nox_scalar  0;
    'NOx'               {'NO2','NO'}        []; %family conservation

    %Alkanes
    'CH4'               4920*org_scalar                 1; %methane
    'C2H6'              123*org_scalar                  1; %ethane
    'C3H8'              63*org_scalar                   1; %propane
    'NC4H10'            25*org_scalar                   1; %n-butane
    'IC4H10'            15*org_scalar                   1; %isobutane
    'NC5H12'            10*org_scalar                   1; %n-pentane
    'IC5H12'            11*org_scalar                   1; %isopentane
    'NC6H14'            4*org_scalar                    1; %n-hexane
    'M2PE'              3*org_scalar                    1; %2-methylpentane
    'M3PE'              2*org_scalar                    1; %3-methylpentane
    'M22C4'             0.3*org_scalar                  1; %2,2-dimethylbutane
    'M23C4'             1.6*org_scalar                  1; %2,3-dimethylbutane
    'NC7H16'            1.9*org_scalar                  1; %n-heptane
    'M2HEX'             0.8*org_scalar                  1; %2-methylhexane
    'M3HEX'             1.2*org_scalar                  1; %3-methylhexane
    'NC8H18'            0.8*org_scalar                  1; %n-octane
    'NC9H20'            0.2*org_scalar                  1; %n-nonane
    'NC10H22'           0.2*org_scalar                  1; %n-decane
    'CHEX'              1.8*org_scalar                  1; %cyclohexane

    %Alkenes & alkyne
    'C2H4'              1.1*org_scalar                  1; %ethylene
    'C3H6'              0.1*org_scalar                  1; %propylene
    'BUT1ENE'           0*org_scalar                    1; %1-butene
    'CBUT2ENE'          0*org_scalar                    1; %cis-2-butene
    'TBUT2ENE'          0*org_scalar                    1; %trans-2-butene
    'PENT1ENE'          0*org_scalar                    1; %1-pentene
    'CPENT2ENE'         0*org_scalar                    1; %cis-2-pentene
    'TPENT2ENE'         0*org_scalar                    1; %trans-2-pentene
    'HEX1ENE'           0*org_scalar                    1; %1-hexene
    'C2H2'              1.8*org_scalar                  1; %acetylene

    %Aromatics
    'BENZENE'           1.2*org_scalar                  1; %benzene
    'TOLUENE'           1.2*org_scalar                  1; %toluene
    'OXYL'              0.1*org_scalar                  1; %o-xylene
    'MXYL'              0.2*org_scalar                  1; %m-xylene
    'PXYL'              0.1*org_scalar                  1; %p-xylene
    'EBENZ'             0.1*org_scalar                  1; %ethylbenzene
    'PBENZ'             0*org_scalar                    1; %propylbenzene
    'IPBENZ'            0*org_scalar                    1; %isopropylbenzene
    'TM123B'            0.1*org_scalar                  1; %1,2,3-trimethylbenzene
    'TM124B'            0*org_scalar                    1; %1,2,4-trimethylbenzene
    'TM135B'            0*org_scalar                    1; %1,3,5-trimethylbenzene
    'OETHTOL'           0*org_scalar                    1; %1-ethyl-2-methylbenzene
    'METHTOL'           0*org_scalar                    1; %1-ethyl-3-methylbenzene
    'PETHTOL'           0*org_scalar                    1; %1-ethyl-4-methylbenzene
    'STYRENE'           0*org_scalar                    1; %styrene

    %Alcohols
    'CH3OH'             10*org_scalar                   1; % methanol 10
    'C2H5OH'            0.3*org_scalar                  1; % ethanol 0.3
    'IPROPOL'           2.3*org_scalar                  1; % isopropanol 2.3

    %Carbonyls
    'HCHO'              6.5*ald_scalar                  1; %formaldehyde
    'CH3CHO'            2.9*ald_scalar                  1; %acetaldehyde
    'C2H5CHO'           0*ald_scalar                    1; %propionaldehyde
    'C3H7CHO'           1.3*ald_scalar                  1; %butyraldehyde
    'C4H9CHO'           0*ald_scalar                    1; %valeraldehyde
    'ACR'               1.6*ald_scalar                  1; %acrolein
    'MACR'              0.6*ald_scalar                  1; %methacrolein
    'C4ALDB'            0*ald_scalar                    1; %crotonaldehyde
    'BENZAL'            4.4*ald_scalar                  1; %benzaldehyde
    'CH3COCH3'          3*ald_scalar                    1; %acetone
    'MEK'               0.7*ald_scalar                  1; %methyl ethyl ketone
    'CYHEXONE'          0*ald_scalar                    1; %cyclohexanone
    };

%% CHEMISTRY
% The chemistry mechanism is MCM v3.3.1.

ChemFiles = {...
   'MCMv331_K(Met)';
   'MCMv331_J(Met,2)'; %Jmethod flag of 0 specifies default MCM parameterization, 2 uses albedo
   'mcm_subset';
   };

%% DILUTION CONCENTRATIONS
%{
  Dilution rates for compounds are determined by the background concentrations
  and the value of kdil in Met.  Background concentrations default to zero except for
  ozone.
%}
BkgdConc = {...
    % names     conc
    'O3'        50;  %same as used by Edwards
    };

%% OPTIONS
%{
"Verbose" can be set from 0 - 3; affecting level of run-time output
"EndPointsOnly" is set to 1 because we only want the last point of each step.
"LinkSteps" Can be set to 0 or 1. With 0, concentrations are initialized at start of each step.
                                  With 1, initialize at first step only.
"IntTime" is the integration time for each step, equal to the spacing of the data (60 minutes).
"TimeStamp" is set to the hour-of-day for observations.
"SavePath" give the filename only (in this example); the default save directory is the UWCMv3\Runs folder.
"FixNOx" forces total NOx to be reset to constrained values at the beginning of every step.
%}

ModelOptions.Verbose        = 2;
ModelOptions.EndPointsOnly  = 1;
ModelOptions.LinkSteps      = 1;
ModelOptions.IntTime        = 3600; % 3600 seconds/hour
ModelOptions.TimeStamp      = Inputs.Time;
ModelOptions.SavePath       = 'J16bOutput';
% ModelOptions.FixNOx       = 1; % if you use this, disable family conservation above.


%% INPUT REPLICATION AND INTERPOLATION

% INTERPOLATION
% The InputInterp function allows for interpolation to a finer time resolution.

if 0
    dt = 1800; % time spacing, seconds

    Time_interp = (0:dt:(86400-dt))'/3600; %interpolation timebase, fractional hours (to match Inputs.Time)
    circularFlag = 1; % time wraps around at midnight
    [Met,InitConc,BkgdConc] = ...
        InputInterp(Met,InitConc,BkgdConc,Inputs.Time,Time_interp,circularFlag);
    ModelOptions.TimeStamp = Time_interp;
    ModelOptions.IntTime = dt;
end

% REPLICATION
% Replicate the run over 4 days, to achieve three days of "spin-up."
%  "repIndex" will be used to separate the days with SplitRun later.
nRep = 4; %number of days to repeat
[Met,InitConc,BkgdConc,repIndex] = InputReplicate(Met,InitConc,BkgdConc,nRep);
ModelOptions.TimeStamp = repmat(ModelOptions.TimeStamp,nRep,1);

%% MODEL RUN

% Now we run the model.  It takes either 2 or 4 minutes on two different  MacBook Pros.

S = F0AM_ModelCore(Met,InitConc,ChemFiles,BkgdConc,ModelOptions);

%% OUTPUT OF RESULTS

% This splits the run into four separate days
SplitRun(S,'custom',repIndex)

% Write out temperature, relative humidity, the three scalar multipliers, solar zenith angle, and concentrations

writematrix(Inputs.T,  'J16b_ATB_T.csv')
writematrix(Inputs.RH, 'J16b_ATB_RH.csv')
writematrix(org_scalar,'J16b_ATB_Org_scalar.csv')
writematrix(nox_scalar,'J16b_ATB_Nox_scalar.csv')
writematrix(ald_scalar,'J16b_ATB_Ald_scalar.csv')
writecell(Met(4,2),    'J16b_ATB_SZA.csv')
writestruct(S.Conc,    'J16b_ATB_conc.xml')

% The S structure is also saved in J16bOutput
