if ~app.NoPlot 

    clf; set(figure(1), 'Position', [0, 800, 1200, 600]);

    subplot(2,1,1);

    S_Plot = rho*DelS;

    S_Plot(1:Nbcs) = NaN;   

    colormap(mk_colmap);

    surf(XP*h/1e3, YP*h/1e3, S_Plot(NodeRho(2:app.Ny+1, 2:app.Nx+1)), 'LineStyle','none');

    axis tight

    view(2); caxis([0,DelS]); colorbar;

    shading interp    

    TT(1) = xlabel('x(km)');

    TT(2) = ylabel('z(km)');

    TT(3) = zlabel('Salinity(psu)');

    TT(4) = title(sprintf(...

        'Salinity Variation (unit: psu) (Nx=%g, Ny=%g, dt=%g, time=%3.3gh)',...

        app.Nx, app.Ny, app.dt*h^2/nu_x, T*h^2/nu_x/3600));

    formatplot([], gca, TT);

    subplot(2,1,2);

    P_Plot = P*(nu_x/h)^2;

    P_Plot(1:Nbcs) = NaN;   

    colormap(mk_colmap);

    surf(XP*h/1e3, YP*h/1e3, P_Plot(NodeP(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('Normalized Pressure(Pa/(kg/m^3))');

    TT(4) = title('Normalized Pressure(Pa/(kg/m^3)) (Without hydrostatic part)');

    formatplot([], gca, TT);

    drawnow;

    set(figure(1), 'PaperPosition', [0, 0, 12, 6]);

    print('-dpng', '-r400', sprintf('%s/%s/%03d', d0, app.saveimg, k/app.PlotIntrvl));

end;

%% Save data

save(sprintf('%s/%s/%03d', d0, app.savemat, k/app.PlotIntrvl),'T','P','u','v','rho', ... 

     'Pi','ui','vi','rhoi','Yi');

%% Display current progress

display(sprintf('%d time steps completed, %d time steps remaining.', k, Nt-k)); 

toc; 

display(sprintf('Current time: %04d/%02d/%02d %02d:%02d:%02d\n', round(clock)));