% This script solves the steady state equations for the equatorial
% atmospheric response to prescribed sea surface saturation entropy
% anomalies, in the case of vanishing momentum drag. Three maps are
% produced. 
%--------------------------------------------------------------------------
file='Gill_params.csv';
%--------------------------------------------------------------------------
warning('off')
% Read Parameters (see text for definitions):
%T=readtable('Gill_params.csv','ReadRowNames',1);
[T]=main();
alpha=T{1,1};        
Gm=T{1,2};
chi=T{1,3};
chi2=T{1,4};         % (Used only in calculation of hm)
delta=T{1,5};
k=T{1,6};            % Zonal wavenumber ('n' in text)
y0=T{1,7};
dx=T{1,8};
dy=T{1,9};
delx=T{1,10};
dely=T{1,11};
cmp=char(T{1,12});
%--------------------
% Grid specification:
y=-3:dy:3;
x=0:dx:pi;
%
%  Spatial distribution of sea surface saturation entropy
%
s0=exp(-1.5.*(y-y0).^4);
%
% Plotting parameters:
YLIM=[-2 2];            % Truncation of graphs in y direction
%-------------------------------------------------------------------------
close('all')
qy=max(size(y));
qx=max(size(x));
delX=round(delx/(x(2)-x(1)));
delY=round(dely/(y(2)-y(1)));
[~,n]=min(y.^2);
if alpha == 0       % Special case where simple analytic solution exists
    s=s0.*y.^2./(y.^2-1i.*Gm.*k/chi);
else  % Do integrals stated in text integrating poleward on each side of equator
    J=s0.*y.^2./(y.^2-1i.*Gm.*k/chi);
    p=zeros(1,qy);
    for j=n+1:qy
        ym=0.5*(y(j-1)+y(j));
        p(j)=p(j-1)-(J(j)-J(j-1))*ym^(-1i*k*Gm/alpha)*exp(0.5*(chi/alpha)*ym^2);
    end   
    %
    for j=n:-1:1
        ym=0.5*(y(j+1)+y(j));
        p(j)=p(j+1)+(J(j+1)-J(j))*ym^(-1i*k*Gm/alpha)*exp(0.5*(chi/alpha)*ym^2);
    end    
    s=J+y.^(1i*k*Gm/alpha).*exp(-0.5*(chi/alpha).*y.^2).*p;
    %
    %seq=0.5*(s(n-1)+s(n+1)); % First-order estimate of equatorial residual
    %
    %seq1=2.5*s(n-1)-2*s(n-2)+0.5*s(n-3); % second-order estimate of equatorial residual
    %seq2=2.5*s(n+1)-2*s(n+2)+0.5*s(n+3); 
    %
    seq1=(8/3)*s(n-1)-2.5*s(n-2)+s(n-3)-s(n-4)/6; % third-order estimate of equatorial residual
    seq2=(8/3)*s(n+1)-2.5*s(n+2)+s(n+3)-s(n+4)/6;
    %
    seq=0.5*(seq1+seq2);
    %
    s=s-seq;  % Subtract estimate of equatorial residual from all s values
    s(n)=0;   % Set equatorial value of s exactly to zero
end
%
% Calculate u,v, w and hm
%
v=-1i*k*s./(y+1e-6);
%
u=zeros(1,qy);
u(2:end)=(s(2:end)-s(1:end-1))./(dy*(y(2:end)+y(1:end-1))); %  u is calculated at intermediate points in y
u(1:end-1)=(u(1:end-1)+u(2:end));  % Interpolate u back to main y points
%
u(n-4:n-1)=u(n-5);u(n+1:n+4)=u(n+5); % This smoothes over anomalous behavior near equator
u(n)=0.5*(u(n-1)+u(n+1));
%
w=(chi*(s0-s)-alpha*u)/Gm;  
hm=w*(1-Gm)+chi2*s;
%
% Now calculate full fields and make plots
%
ek=1i*exp(1i*2.*x);    % Regardless of value of k, plot over one full wavelength in x
EK=repmat(ek,[qy,1]);
S=real(repmat(s',[1,qx]).*conj(EK));
S0=real(repmat(s0',[1,qx]).*conj(EK));
U=real(repmat(u',[1,qx]).*conj(EK));
V=real(repmat(v',[1,qx]).*conj(EK));
W=real(repmat(w',[1,qx]).*conj(EK));
HM=real(repmat(hm',[1,qx]).*conj(EK));
%
tit2=strcat('\alpha=',num2str(alpha),', G_m=',num2str(Gm),', \chi=',num2str(chi), ...
    ', \chi_2=',num2str(chi2),', \delta=',num2str(delta),', k=',num2str(k),', y_0=',num2str(y0));
figure('units','normalized','position',[0.1 0.4 0.4 0.5])
uc=colormap(cmp);
%uc=uc(60:197,:);       % Good for 'jet'
%uc=uc(120:256,:);       % Good for 'bone' and 'gray'
contourf(x,y,W,20);
cl=get(gca,'CLim');
hold on
contour(x,y,max(S0,0),10,'k','linewidth',1)
contour(x,y,min(S0,0),10,'--k','linewidth',1)
hold off
set(gca,'CLim',cl)
set(gca,'ylim',YLIM)
set(gca,'fontsize',13,'fontweight','bold')
xlabel('x')
ylabel('y','rotation',0);
title([{'s_0 (black contours), w (shading)'};{tit2}])
colormap(uc)
%
figure('units','normalized','position',[0.5 0.4 0.4 0.5])
%
colormap(uc)
contourf(x,y,S,20);
hold on
quiver(x(1:delX:end),y(1:delY:end),U(1:delY:end,1:delX:end), ...
    V(1:delY:end,1:delX:end),'k','linewidth',1.4)
hold off
set(gca,'ylim',YLIM,'xlim',[min(x) max(x)])
set(gca,'fontsize',13,'fontweight','bold')
xlabel('x')
ylabel('y','rotation',0);
title([{'s (shading), surface winds (arrows)'};{tit2}])
%
figure('units','normalized','position',[0.3 0.05 0.35 0.45])
%
colormap(uc)
contourf(x,y,HM,20);
set(gca,'ylim',YLIM,'xlim',[min(x) max(x)])
set(gca,'fontsize',13,'fontweight','bold')
xlabel('x')
ylabel('y','rotation',0);
title([{'h_m (mid-tropospheric moist static energy)'};{tit2}])
%
warning('on')