%   ES210_MATLAB_VB10b1.m

%   Version:  08 November 2024
%   Added new End User Input to address Teaching with MATLAB Workshop
%   comment -> allow for a variable display of temperature profiles. 

%   Version:  09 October 2024
%   Minor corrections for clarity

%   Version:  02 December 2020
%   Minor correction to title in figure(102) - there is no diffusion here.

%   Version:  03 Dec. 2016
%   M-file for calculation of temperature advection (one dimensional motion):
%     - User defined time-step and duration
%     - East/West Domain:  3000 km, dx = 4 km
%     - Incorporates Courant-Friedrichs-Levy (CFL) stability condition

%%   Header

clear variables;
close all;

format shortg;
format compact;
clc;

%%   Burger's Equation Model (i.e., without diffusion)

%   Advection Scheme selection
Scheme = input('User Input - Select Finite-Difference Scheme; 0 = Leapfrog, 1 = Adams-Bashforth:  ');

%   Experimental parameters:  Time
Hours = input('User Input - Time Step; default is 2 hours:  ');   %  User Input
    if isempty(Hours) == 1
        Hours = 2;
        dt = Hours*3600;   %  Convert input from hours into seconds
        disp('Time Step unspecified, use 2 hours');
    else
        dt = Hours*3600;   %  Convert input from hours into seconds
    end
Weeks = input('User Input - Duration (weeks); default is 8:  ');   %  User Input
    if isempty(Weeks) == 1
        Weeks = 8;
        disp('Duration unspecified, use 8 weeks');
    end
Duration = Weeks*7*24*60*60;   %  Duration of Advection Event (sec.)
Nt = floor(Duration/dt) + 1;   %  Number of time steps
disp('Total number of Time Steps:'),Nt 
LimitFigures = input('User Input - View Temperature Profiles, input an integer < Time Steps; default is a figure every 25th time step: ');   %  User Input
    if isempty(LimitFigures) == 1
        LimitFigures = 25;
        disp('Figure Display unspecified, use default');
    end
    
%   Experimental parameters:  Horizontal spatial dimension
Lx = 3000000;   %  Basin length horizontal scale (m) --> 3000 km
dx = 4000;   %  Horizontal Grid Spacing (m)
% disp('Horizontal Spatial Resolution (m): '),2*dx
Nx = floor(Lx/dx) + 1;   %  Number of horizontal grid points
x = linspace(0,Lx,Nx);   %  Calculate the grid points

%   Experimental parameters:  Initial Temperature Distribution
T_model = linspace(22,22,Nx);   %  Base model temperature profile
T_model(1,20:45) = 28;   %  Create the western tropical Pacific Ocean warm pool initial distribution
T_model = repmat(T_model,Nt,1);   %  Duplicate the initial temperature profile for all times

figure(1);
plot(x,T_model(1,:),'k.-');
axis([0, Lx, 20, 30]);
xlabel('Spatial dimension; Eastward (m)')
ylabel('Sea Surface Temperature ( ^{o}C)')
title('Initial Temperature Distribution:  western tropical Pacific warm pool');

u = 0.25;   %  Advection velocity (m/s); Eastward

%  Courant-Friedrichs-Lewy stability condition.
CFL = u*dt/dx;
disp('CFL Stability Condition; CFL(Adams-Bashforth) < 0.72, Leapfrog CFL < 1.0: '),CFL

%%   Integrate forward in time for all grid points
time = zeros(1,Nt);   %   Preallocate for speed

for n = 3:Nt
    time(n) = (n-1)*dt;
    for i = 2:Nx-1
    %   Student Input:  Recursive Formula for the Time-Differencing Scheme HERE
        if Scheme == 0
            T_model(n+1,i) = T_model(n-1,i) - CFL*(T_model(n,i+1) - T_model(n,i-1));   %  Leapfrog and Centered-in-Space Methods (2nd Order)
        end
        if Scheme == 1
            T_model(n+1,i) = T_model(n,i) - (CFL/24)*(23*(T_model(n,i+1) - T_model(n,i-1)) - 16*(T_model(n-1,i+1) - T_model(n-1,i-1)) + 5*(T_model(n-2,i+1) - T_model(n-2,i-1)));   %  Adams-Bashforth Scheme (3rd Order)
        end
        
    %   Student Input:  Upstream and Downstream Advective Boundary Conditions HERE
        T_model(n+1,2) = T_model(n,2) + CFL*(T_model(n,2) - T_model(n,1));   %  Open Boundary Condition (upstream) with Forward differencing Scheme (1st Order) Advection
        T_model(n+1,Nx) = T_model(n,Nx) - CFL*(T_model(n,Nx) - T_model(n,Nx-1));   %  Open Boundary Condition (downstream) with Forward differencing Scheme (1st Order) Advection
    end
    
    %   Display the temperature profile every 25th time step during the simulation.
    if mod(n,LimitFigures) == 0
        figure(n);
        plot(x,T_model(1,:),'k.-');
        plot(x,T_model(n+1,:),'r.-');
        axis([0, Lx, 20, 30]);
        xlabel('Spatial dimension; Eastward (m)')
        ylabel('Sea Surface Temperature ( ^{o}C)')
        if Scheme == 0
            text(0.1*Lx,29,'Leapfrog');
        end
        if Scheme == 1
            text(0.1*Lx,29,'Adams-Bashforth');
        end
        title(['Temperature Advection:  western tropical Pacific warm pool; Elapsed: ',num2str(time(n)/(24*3600)),' days']);
    end
end

%%   Three-dimensional display of results
[Time,X] = meshgrid(x/1000,time/(24*3600));   %  Convert time into days and streamwise distance to km for display purposes

figure(101);
mesh(Time,X,T_model(1:end-1,:))
C = colorbar; C.Label.String = 'Water Temperature ( ^{o}C)';   %  Label the colorbar
view(-15,65);
xlabel('Horizontal Distance (km)')
ylabel('Time (days)')
zlabel('Temperature ( ^{o}C)')
title(['Model Advection (without Diffusion) of western tropical Pacific warm pool; u = ',num2str(10*u),' cm s^{-1} - along Equator'])

[Days,Longitude] = meshgrid(170 + x/(111*1000),time/(24*3600));   %  Convert axis vectors into arrays and convert units:  Longitude (degrees, East) and hours 

figure(102);
mesh(Days,Longitude,T_model(1:end-1,:))
C = colorbar; C.Label.String = 'Water Temperature ( ^{o}C)';   %  Label the colorbar
view(-0,90);   %  View the 3D plot from specified azimuth and elevation angles
axis ij;   % Reverse the coordinate system so that the  _y_ values increase from top to bottom.
axis tight;   %  Set the axis limits to equal the range of the data 
xlabel('Longitude, East ( ^o)')
ylabel('Elapsed Time (days)')   %  Left axis label in one row
zlabel('Temperature ( ^{o}C)')
title(['Model Advection (without Diffusion) of western tropical Pacific warm pool; U = +',num2str(100*u),' cm s^{-1} - along Equator'])