% Model 1: test of convective parameterization
%
% 1 single O/A column with fixed Fo and Fa
%
clear
Param

Nt = 10000;
Dt = DAY;
Time = [0:1:Nt-1]*Dt;

% Oceanic & Atm. cooling
Fa = 0; %Wm-2
Fo = 0; %Wm-2
Te = TEMI;

Ta(1) = 300; %initial Atm. temp. (K)
To(1) = 280; %initial (oceanic) surf. temp. (K)
for n = 1:Nt-1
   [Fs(n),Feva(n)] = calc_Fs(Te,Ta(n),To(n)); %net surface heat flux (W m-2)
   Ft(n) = calc_Ft(Te,Ta(n),To(n)); %net top-of-atm heat flux (W m-2)
   Tend_atm = (Fs(n) + Ft(n) - Fa)/HCA;
   Tend_oce = -(Fs(n) + Fo)/HCO;
   Ta(n+1) = Ta(n) + Dt*Tend_atm;
   To(n+1) = To(n) + Dt*Tend_oce;
end
[Fs(Nt),Feva(Nt)] = calc_Fs(Te,Ta(Nt),To(Nt));
Ft(Nt) = calc_Ft(Te,Ta(Nt),To(Nt));

% Final temperature
  To_end = To(end),
  Ta_end = Ta(end),
  DT_end = To_end-Ta_end,
  DTcrit,
  
figure(1), clf, 
subplot(2,1,1)
plot(Time/MONTH,Ta,'b',Time/MONTH,To,'g');
%axis([0 1 100 300]);
subplot(2,1,2)
plot(Time/MONTH,Fs,'c',Time/MONTH,Feva,'b',Time/MONTH,Ft,'m',Time/MONTH,Fs+Ft,'k');
%axis([0 160 -300 300]);

figure(2), clf,
plot(Time/MONTH,To-Ta-DTcrit_conv,'g');