% script to make terrain for WRF LES runs
% this script constructs the topography, then writes it out as an ordered
% list to be read in my module_initialize_les.F

% 

close all
clear all


% define the grid size (note that this MUST match with the namelist setup
% for WRF!!!!!!!)

xdim=1500;
ydim=100;
dx=30;
dy=30;

plane=300; % elevation of the plane
%basin=500; % elevation of the basin floor

hm=500; % height of the confining ridges;
mtnleft=xdim./4;
icm=xdim./2;
mtnright=xdim.*(3/4);
xa1=30; % inner mountain width parameter
xa2=100; % outer mountain width
hgt=zeros(xdim,ydim);

for i=1:xdim-1;
    for j=1:ydim;
        if i< mtnleft
            hgt(i,j)=plane+(hm-plane)./(1+((i-mtnleft)./xa2).^2);
        elseif i>=mtnleft & i< icm
            hgt(i,j)=hm./(1+((i-mtnleft)./xa1).^2);
        elseif i>= icm & i<=mtnright
            hgt(i,j)=hm./(1+((i-mtnright)./xa1).^2);
        else
            hgt(i,j)=plane+(hm-plane)./(1+((i-mtnright)./xa2).^2);
        end
    end
end
hgt(end,:)=hgt(1,:);

f=figure(1);
set(f,'color','w')
set(gcf,'units','inches');
set(gcf,'pos',[13 4 18 3])
subplot(1,2,1)
%plot((1:xdim).*dx,hgt(:,10));
h=area((1:xdim).*dx,hgt(:,10));
set(h,'FaceColor',[.75 .75 .75])
set(h,'linewidth',2);
set(gca,'xtick',[0:5000:(xdim.*dx)]);
xtick=get(gca,'xtick');
xlab=xtick./10^3;
set(gca,'xticklabel',{xlab});
xlabel('X Dist [km]');
ylabel('Altitude [m]');
ylim([0 hm+50]);
xlim([0 (xdim.*dx)]);


% compute the slopes to make sure we get something reasonable for the SLV

slope=gradient(squeeze(hgt(:,10)),dx);
slpang=atand(abs(slope));
subplot(1,2,2);
plot((1:xdim).*dx,slpang,'-k','linewidth',2);
%area((1:xdim).*dx,slpang);
set(gca,'xtick',[0:5000:(xdim.*dx)]);
xtick=get(gca,'xtick');
xlab=xtick./10^3;
set(gca,'xticklabel',{xlab});
xlabel('X Dist [km]');
ylabel('slope angle [deg]');
%ylim([0 hm+50]);
xlim([0 (xdim.*dx)]);
%ylim([0 hm+50]);

%terrain_mat=hgt(:);

save('TERRAIN.mat','hgt')
return
%%
%fname='TERRAIN';
%hdr=[xdim ydim dx dy];
%dlmwrite(fname,hdr,'delimiter','\t','precision', '%.0f');
%dlmwrite(fname,terrain_mat,'-append','delimiter','\t','precision', '%.2f');

%% NOW ADD THE SOUNDING TO THE TOPOGRAPHY
