%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Box modeling for CO2 with advection and entrainment
% made by S. Ono, Sep 10 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set up
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt = 0.01; % time step in days
T = 100; % total time for integration in days
N = floor(T/dt); % total number of integration step
dx = 10000; % lateral length of the box (m), dx = dy
Mair = 28.97*1E-3; % molar weight of air (kg/mol)
g = 9.807; % acceleration of gravity (m/s-2)
ps = 1020*100; % surface pressure in Pa % it was 10000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0 = 0; % initial time
mBL = 400E-6; % initial CO2 mixing ratio for Boundary Layer
mENT = 370E-6; % CO2 mixing ratio for entranment zone
mUP = 380E-6; % CO2 mixing ratio for upwind box
wind = 4;%4; % wind speed, m/s
tau = dx/(wind*60*60*24); % residence time wrt lateral advection (day);
q = 10; % source flux mol CO2/day-1 m-2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HARMONIC FIT FOR BOUNDARY LAYER HEIGHT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ta = t0:dt:dt*(N+1);
wt = 2*pi*ta';
A = [sin(1*wt) cos(1*wt) sin(2*wt) cos(2*wt) sin(3*wt) cos(3*wt) sin(4*wt) cos(4*wt)];
a = [ 75.2235; 62.8391; -29.1925; 8.0102; -7.4486; 6.5066; 1.9185; -11.6844];
a0 = 931.7713;
ph = 100*(A*a+a0); % boundary layer Pressure (Pa) factor 100 is for unit conversion form mbar to Pascal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = t0;
for n=1:N
dp = ph(n+1)-ph(n);
source = Mair*g/(ps-ph(n))*q*dt; % source flux
adv = 1/tau*(mUP-mBL)*dt; % advection term,
ent = -dp/(ps-ph(n))*(mENT - mBL)*(dp<0); % entrainment term, zero if dhdt<0
mBL= mBL + adv + ent + source;
t = t+dt ;
time(n) = t;
MBL(n) = mBL;
ADV(n) = adv;
ENT(n) = ent;
SOURCE(n) = source;
end
figure(1)
subplot(3,1,1)
plot(time(N-300:N), MBL(N-300:N)*1E+6)
xlabel('time (days)')
ylabel('mCO2 (ppm)');
subplot(6,1,3)
plot(time(N-300:N), ADV(N-300:N)/dt)
xlabel('time (days)')
ylabel('Adv, Flux CO2 ');
subplot(6,1,4)
plot(time(N-300:N), ENT(N-300:N)/dt)
xlabel('time (days)')
ylabel('Ent, Clux CO2');
subplot(6,1,5)
plot(time(N-300:N), SOURCE(N-300:N)/dt)
xlabel('time (days)')
ylabel('source flux');
subplot(6,1,6)
plot(time(N-300:N), -ph(N-300:N)/100)
xlabel('time (days)')
ylabel('-P_h(mbar)');
%ylim([-1000,-800])