%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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])



