%% Bottom Gravity Current Test Case Post-Processing Script

% This test case is based on the paper "Dynamics of Two-Dimensional Turbulent Bottom
% Gravity Currents" by Ozgokmen and Chassignet in 2002 published in
% Journal of Physical Oceanography.

% Date: Jan. 2014
% Author: Jing Lin

%% Set up data loading and images saving directory


%Checking variance

   rhosumsq = 0;
   rhosum   = 0;
   usumsq = 0;
   usum   = 0;
   vsumsq = 0;
   vsum   = 0;
   Ntot   = 0;

   for i = 0 : round(Nt/app.PlotIntrvl)

      load(sprintf('%s/%s/%03d', d0, app.savemat, i));

      rhosum = rhosum + rho;
      rhosumsq = rhosumsq + rho.*rho;

      usum = usum + u;
      usumsq = usumsq + u.*u;

      vsum = vsum + v;
      vsumsq = vsumsq + v.*v;

      Ntot = Ntot + 1;

   end;

   rhovar = (rhosumsq - rhosum.*rhosum./Ntot) ./ (Ntot-1);
   uvar = (usumsq - usum.*usum./Ntot) ./ (Ntot-1);
   vvar = (vsumsq - vsum.*vsum./Ntot) ./ (Ntot-1);

   clf; set(figure(1), 'Position', [0, 200, 1200, 1200]);
    
   subplot(3,1,1);
           fld = rhovar;
           fld(1:Nbcs) = NaN;
           xtest = XP*h/1e3;
           ftest = fld(NodeRho(2:app.Ny+1, 2:app.Nx+1));
           ind = find( (xtest<totvarXlim) & (YP<=(b+b1)*0.5) & ~isnan(ftest) );
           rhovarmean = mean(ftest(ind));
           colormap(mk_colmap);
           surf(XP*h/1e3, YP*h/1e3, fld(NodeRho(2:app.Ny+1, 2:app.Nx+1)), 'LineStyle','none');
           axis tight
           view(2); colorbar;
           shading interp    
           TT(1) = xlabel('x(km)');
           TT(2) = ylabel('z(km)');
           TT(3) = zlabel('Density variance');
           TT(4) = title(sprintf(...
               'Density variance (Nx=%g, Ny=%g, dt=%g, time=%3.3gh); meanvar=%3.3g',...
               app.Nx, app.Ny, app.dt*h^2/nu_x, i*app.PlotIntrvl*app.dt*h^2/nu_x/3600, rhovarmean));
           formatplot([], gca, TT, [], font_settings.fontsize);
    
           subplot(3,1,2);
           u_Plot = uvar;
           u_Plot(1:Nbcs) = NaN;   
           xtest = XU*h/1e3;
           ftest = u_Plot(Nodeu(2:app.Ny+1, 2:app.Nx));
           ind = find( (xtest<totvarXlim) & (YU<=(b+b1)*0.5) & ~isnan(ftest) );
           uvarmean = mean(ftest(ind));
           colormap(mk_colmap);
           surf(XU*h/1e3, YU*h/1e3, u_Plot(Nodeu(2:app.Ny+1, 2:app.Nx)), 'LineStyle','none');
           axis tight
           view(2); colorbar;
           shading interp    
           TT(1) = xlabel('x(km)');
           TT(2) = ylabel('z(km)');
           TT(3) = zlabel('u variance');
           TT(4) = title(['u variance; meanvar=',num2str(uvarmean)]);
           formatplot([], gca, TT, [], font_settings.fontsize);
    
           subplot(3,1,3);
           v_Plot = vvar;
            v_Plot(1:Nbcs) = NaN;   
           xtest = XV*h/1e3;
           ftest = v_Plot(Nodev(2:app.Ny, 2:app.Nx+1));
           ind = find( (xtest<totvarXlim) & (YV<=(b+b1)*0.5) & ~isnan(ftest) );
           vvarmean = mean(ftest(ind));
           colormap(mk_colmap);
           surf(XV*h/1e3, YV*h/1e3, v_Plot(Nodev(2:app.Ny, 2:app.Nx+1)), 'LineStyle','none');
           axis tight
           view(2); colorbar;
           shading interp    
           TT(1) = xlabel('x(km)');
           TT(2) = ylabel('z(km)');
           TT(3) = zlabel('w variance');
           TT(4) = title(['w variance; meanvar=',num2str(vvarmean)]);
           formatplot([], gca, TT, [], font_settings.fontsize);

   print('-dpng', sprintf('%s/%s/variance', d0, app.saveimg));

