%% ----------------------------------------
%% EPS131: Intro to Phys. Oceanog.
%% Section-07: ex. 1
%% Sverdrup relation and zonal flow
%% Laure Zanna
%% ----------------------------------------
clear all; close all;

%parameters
tau_0=0.15;  rho_0=1025;
L_x=5000*10^3 ; L_y=3330*10^3 ; delta_x=100*10^3;
R=6371*10^3; omega=7.3*10^-5; theta_0=35*pi/180;
beta=(1/R)*2*omega*cos(theta_0); H=400;

% domain: latitude in kms
y=0:L_y/20:L_y;

% -----------------
% wind stress
% -----------------
tau_x=-tau_0*cos(pi*y/L_y);
figure(1); subplot(1,2,1); plot(tau_x,y); title('\tau [N/m^2]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 

% -----------------------
% curl of the wind stress
% -----------------------
curl_tau=-(tau_0*pi/(L_y))*sin(pi*y/L_y);
figure(1); subplot(1,2,2); plot(curl_tau,y); title('\nabla \times \tau [N/m^3]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 

% ---------------------------------------------------
% V_int= vertically integrated N-S transport in the interior
% ---------------------------------------------------
V_int=-(tau_0*pi/(rho_0*beta*L_y))*sin(pi*y/L_y);
figure(2); subplot(2,2,1);plot(V_int,y); title('V [m^2/s]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 

% ----------------------------------------------------------
% U_int= vertically integrated E-W transport in the interior
% ----------------------------------------------------------
x=0:L_x/100:L_x;
[X,Y]=meshgrid(x,y);
U=(tau_0*(pi^2)/(rho_0*beta*L_y^2))*cos(pi*Y/L_y).*(X-L_x);
figure(2); subplot(2,2,2); contourf(x,y,U); title('U [m^2/s]');colorbar;
set(gca,'XLim',[0 L_x],'YLim',[0 L_y])
set(gca,'XTick',[0 L_x],'YTick',[0 L_y/2 L_y]);
set(gca,'XTicklabel',{'60W','10W'},'YTicklabel',{'20N','35N','50N'}); 

% ------------------------------------------------------------------
% streamlines of the interior flow : U_int=-dpsi/dy, V_int=dpsi/dy
% ------------------------------------------------------------------
psi=-(tau_0*pi/(rho_0*beta*L_y))*sin(pi*Y/L_y).*(X-L_x);
figure(2); subplot(2,2,3); contourf(x,y,psi/10^6,20); title('\psi [Sv]'); colorbar;
set(gca,'XLim',[0 L_x],'YLim',[0 L_y])
set(gca,'XTick',[0 L_x],'YTick',[0 L_y/2 L_y]);
set(gca,'XTicklabel',{'60W','10W'},'YTicklabel',{'20N','35N','50N'}); 

% ---------------------------------------------------------
% V_wb: Vertically integrated transport in the western boundary
% ---------------------------------------------------------
V_wb=(tau_0*pi*(L_x-delta_x)/(delta_x*rho_0*beta*L_y))*sin(pi*y/L_y);
figure(2); subplot(2,2,4); plot(V_wb,y); title('V=vdz in the w.b. [m^2/s]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 

% ----------------------------------------------
% Total transport [Sv] in the interior: V_int * width
% ----------------------------------------------
V_total=-(tau_0*pi*(L_x-delta_x)/(rho_0*beta*L_y))*sin(pi*y/L_y);
figure(3); subplot(1,3,1); plot(V_total/10^6,y); title('Total interior [Sv]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 

% -----------------------------------
% velocity in the interior: v_int [m/s]
% -----------------------------------
v=-(tau_0*pi/(rho_0*beta*L_y*H))*sin(pi*y/L_y);
figure(3); subplot(1,3,2);plot(v,y); title('v interior[m/s]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 

% ---------------------------------------------
% velocity in the western boundary: v_wb [m/s]
% --------------------------------------------
v_wb=Vtotal_wb/H
figure(3); subplot(1,3,3);plot(v_wb,y); title('v in the w.b [m/s]')
set(gca,'YLim',[0 L_y])
set(gca,'YTick',[0 L_y/2 L_y]);
set(gca,'YTicklabel',{'20N','35N','50N'}); 
