clear all
tic;
adc_cell_in;
rng(seed);

%--- set up constants
delta = sqrt(kappa/abs(U)/(L/pi));
alpha = log(qmax/qmin)/L;        % qs(y) = qmax*exp(-alpha*y)
q_thres = qmax*(1-exp(-alpha*delta));
dd = sqrt(2*kappa*dt);
kUdt = ku*U*dt;
i_tsave = 0;

n_iter = round(t_tot/dt);
fsave_step = round(f_save/dt);
tsave_step = round(t_save/dt);
tskip_step = round(t_skip/dt);
n_tsave = floor(n_iter/tsave_step);
n_tavg_save = ceil((n_iter-tskip_step+1)/fsave_step);
nx_c = round(L/bin);

%--- arrays for flux calculation
dyflux = L/(n_flux_pos+1);        % width of y-partition
ypos = (1:1:n_flux_pos)*dyflux;   % y-position to calculate flux
flux = zeros(1,n_flux_pos);

%--- initialize arrays
q_avg = zeros(n_tsave,1);
q_std = zeros(n_tsave,1);
moist_in_t = zeros(n_tsave,1);
rain_x = zeros(nx_c,1);
moist_net_t = zeros(n_tsave,1);
rain_bl = zeros(nx_c,1);

x = rand(1,np)*L;
y = rand(1,np)*L;

%--- saturation profile
%Tmax = 26;                                     % C-C profile
%Tmin = -50;                                    % C-C profile
%T = Tmax - (Tmax - Tmin)*y/L;                  % C-C profile
%qmin = 0.3619*exp(17.67*Tmin/(Tmin+243.3));    % C-C profile
%qmax = 0.3619*exp(17.67*Tmax/(Tmax+243.3));    % C-C profile

qs = qmax*exp(-alpha*y);                        % exp profile
%qs = 0.3619*exp(17.67*T./(T+243.3));           % C-C profile

qs_ypos = qmax*exp(-alpha*ypos);       % for flux calculation

%--- initial condition for (x,y,q)
q = qs;

fid = fopen('xyq.in','w');

if animate
   figure(1)
%  set(gcf,'visible','off');
   scatter(x,y,20,log10(q),'filled');box on;
   xlabel('$x$','interpreter','latex','fontsize',20);
   ylabel('$y$','rotation',0,'interpreter','latex','fontsize',20);
   caxis([-2 0]);
   hc = colorbar('fontsize',24);axis image; 
   set(hc,'ytick',-2:0.5:0)
   title('$\log_{10} q\,$ at time = 0.0', ...
          'interpreter','latex','fontsize',24)
   xlim([0 L]);ylim([0 L]);
%  print(gcf,'-r125','./png/0.png','-dpng');
   pause(0.001)
end

q_avg0 = mean(q);
q_std0 = std(q);
        
for n = 1 : n_iter

tsave = mod(n,tsave_step);
fsave = mod(n,fsave_step);
cal_tavg = (fsave == 0 & n >= tskip_step);
save_field = cal_tavg;
if tsave == 0; i_tsave = i_tsave + 1; end

%---flux calculation: partition index before advection
if cal_tavg
    yn1 = floor(y/dyflux);
end 

%--- advection (sympletic, Euler-Murayama)
f = x + y;
g = x - y;
f = f + kUdt*sin(ku*g);
g = g + kUdt*sin(ku*f);
x = (f + g)*0.5;
y = f - x;
x = x + dd*randn(1,np);
y = y + dd*randn(1,np);

%--- reflecting boundary condition in x
x = abs(x);
ind = x > L;
x(ind) = 2*L - x(ind);

%--- forcing
ind = y <= 0;
if tsave == 0                 % calculate moisture input
   delta_q = qmax - q(ind);
   moist_in_t(i_tsave) = sum(delta_q);
   ind1 = delta_q > q_thres;
   moist_net_t(i_tsave) = sum(delta_q(ind1));
end
q(ind) = qmax;

%--- flux calculation: correct partition index at bottom
if cal_tavg
   yn1(ind) = 0;
end

%--- reflecting boundary condition in y
y(ind) = -y(ind);
ind = y > L;
y(ind) = 2*L - y(ind);
q(ind) = qmin;               % hit upper boundary: q -> qmin

%--- flux calculation
if cal_tavg
   yn1(ind) = n_flux_pos;    % correct partition index at top
   yn2 = floor(y/dyflux);    % partition index after advection

   for j1 = 1 : n_flux_pos
       ind = yn1 < j1 & yn2 >= j1;           % +ve flux
       flux(j1) = flux(j1) + sum(min(q(ind),qs_ypos(j1)));
       ind = yn1 >= j1 & yn2 < j1;           % -ve flux
       flux(j1) = flux(j1) - sum(q(ind));
   end
end

%--- condensation
qs = qmax*exp(-alpha*y);                    % exp profile
%T = Tmax - (Tmax - Tmin)*y/L;              % C-C profile
%qs = 0.3619*exp(17.67*T./(T+243.3));       % C-C profile
if cal_tavg                   % precipitation calculation
   ind = q > qs;
   moist_cond = q(ind) - qs(ind);
   x_cond = ceil(x(ind)/bin);
   ind1 = y(ind) < delta;
   for j = 1 : nx_c
      ind = (x_cond == j);
      rain_x(j) = rain_x(j) + sum(moist_cond(ind));
      ind = ind & ind1;
      rain_bl(j) = rain_bl(j) + sum(moist_cond(ind));
   end
end
q = min(q,qs);

%--- compute mean and s.d.
if tsave == 0
   q_avg(i_tsave) = mean(q);
   q_std(i_tsave) = std(q);
end

%--- output field data
if save_field
   fwrite(fid,x,'single');
   fwrite(fid,y,'single');
   fwrite(fid,q,'single');
end
    
%--- animation
if animate
   figure(1)
%  set(gcf,'visible','off');
   scatter(x,y,20,log10(q),'filled');box on
   xlabel('$x$','interpreter','latex','fontsize',20);
   ylabel('$y$','rotation',0,'interpreter','latex','fontsize',20);
   caxis([-2 0]);
   hc = colorbar('fontsize',24);axis image; 
   set(hc,'ytick',-2:0.5:0)
   title(['$\log_{10} q\,$ at time = ' num2str(n*dt,'%.2f')], ...
          'interpreter','latex','fontsize',24)
   xlim([0 L]);ylim([0 L]);
%  print(gcf,'-r125',['./png/' num2str(n) '.png'],'-dpng');
   pause(0.001)
end

end
fclose(fid);

%----- Save data
time = [0 tsave_step : tsave_step : n_iter]'*dt;
out = [time [q_avg0;q_avg]];
save q_avg.out out -ascii

out = [time [q_std0;q_std]];
save q_std.out out -ascii

tsave_skip = ceil(t_skip/t_save);
mean_q = mean(q_avg(tsave_skip:end));
save mean_q.out mean_q -ascii

flux = flux*L*L/np/dt/n_tavg_save;
out = [flux' ypos'];
save flux.out out -ascii

time = time(2:end);
moist_in_t = moist_in_t*L*L/np/dt;
out = [time moist_in_t];
save moist_in_t.out out -ascii

moist_in = mean(moist_in_t(tsave_skip:end));
save evap_s.out moist_in -ascii

moist_net_t = moist_net_t*L*L/np/dt;
out = [time moist_net_t];
save moist_net_t.out out -ascii

moist_net = mean(moist_net_t(tsave_skip:end));
save evap_net.out moist_net -ascii

rain_pos = (1:1:nx_c)'*bin+bin/2;
rain_x = rain_x*nx_c*L/np/dt/n_tavg_save;
out = [rain_pos rain_x]; 
save precip_x.out out -ascii

rain_bl = rain_bl*nx_c*L/np/dt/n_tavg_save;
out = [rain_pos rain_bl]; 
save precip_bl_x.out out -ascii

runtime = toc;
save runtime.txt runtime -ascii;

exit
