%MIT OpenCourseWare: https://ocw.mit.edu
%12.810 Dynamics of the Atmosphere (Spring 2023)
%License: Creative Commons BY-NC-SA 
%For information about citing these materials or our Terms of Use, visit: https://ocw.mit.edu/terms.

% Mountain wave forced by a ridge

% Fast fourier transform is used to determine spectral coefficients for the solution. 

% Topography should either be isolated or periodic such that the domain length is a multiple of the wavelength (otherwise get ringing)

% Adapted from Holton and Hakim matlab lee_wave_1.m routine: 
% -fixed issues with k versus kmin, and issue with use of FFT (need positive and negative wavenumbers or else get aliasing)
% -changed to Gaussian ridge for the pset problem

clear

disp('topographic wave for isolated ridge')


%  Set parameters
u0 = input('give mean zonal wind in m/s  ');
Lz = 15e3;           % depth of domain in m
L = 3e3;             % width of ridge in m as measured by standard deviation
Lx = 70e3;           % Lx is the length of the domain in m 
bv = 1.e-4;          % buoyancy frequency squared s^-2
g = 9.81;            % m s^-2

%%*********************************************
NL = 61;              % number of vertical gridpoints
N = 512;              % number of modes for Fourier transform
kmin = 2*pi/(Lx);   % lowest zonal wavenumber 
uxz = zeros(NL,N);
wxz = zeros(NL,N);
zz = linspace(0,Lz,NL);
xx = linspace(0, Lx, N); 
dz = Lz/(NL-1);
s = [1:N];
xm = Lx/2 ;                      % location of ridge 
hx = 5e2*exp(-(xx-xm).^2/2/L^2); % Gaussian ridge of height 500m

% Fourier transform hx(x) to get hn(s)
hn = fft(hx,N);
hn(1) = 0;
%plot(abs(hn), 'o'); return
%  solve for streamfunction for each fourier mode
w=zeros(NL,N);
u=zeros(NL,N);
psi=zeros(NL,N);
for n = 2:N

    ns = n-1;

    % divide into positive and negative wavenumbers
    if ns<N/2
     k = kmin*ns;
    else
     k = kmin*(ns-N);
    end

    %  vertical wavenumber squared for analytic solution for each zonal wavenumber
    m2 = bv/u0^2-k^2;

    % apply radiation condition where applicable
    if k>0 | m2<0  % m2<0 is an evanescent wave
     m = sqrt(m2);
    else
     m = -sqrt(m2);
    end


    %  vertical velocity for each Fourier mode
    wsurf = i*k*u0*hn(n);
    for j=1:NL
        w(j,n) = wsurf*exp(i*m*zz(j));  %vertical profile of w
    end

    % streamfunction contours (approximate trajectories)
    u(:,n) = -m*w(:,n)/k; 
    psi(:,n) = i*u(:,n)/m;

end


for j = 1:NL
    uxz(j,:) = real(ifft(u(j,:),N));    % disturbance zonal wind in x,z space
    wxz(j,:) = real(ifft(w(j,:),N));    % vertical velocity
    psixz(j,:) = real(ifft(psi(j,:),N)) -u0*j*dz;
    psixz_ref(j,:) = 0*real(ifft(psi(j,:),N)) -u0*j*dz;
end
[x, z] = meshgrid(xx,zz); 

disp('max u prime')
disp(max(max(uxz)))

figure(1)
[cs H] = contour(x/1000,z/1000,wxz);
clabel(cs,H)
title('vertical velocity (m/s)')
xlabel( 'horizontal distance (km)'), ylabel('height (km)')

figure(2)
[cs H] = contour(x/1000,z/1000,psixz);
hold on; contour(x/1000, z/1000, psixz_ref, '--');
set(gcf,'renderer','zbuffer')

hold off
title('streamlines')
xlabel( 'horizontal distance (km)'), ylabel('height (km)')
