%% THIS SCRIPT IS USED TO CREATE INITIAL CONDITIONS FOR THE WRF TO TEST
% to test Zhong et al 2003 analytic model against LES runs

% we must both create a terrain matrix and a sounding file for each of the
% Zhong cases

close all
clear all
scrz=get(0,'screensize');


dx=50;
dy=50;
xdim=1000;
ydim=10;
xmax=xdim.*dx;
ymax=ydim.*dy;
xc=xdim./2;
yc=ydim./2;

alpha=15;%degrees;
h=500; % basin depth (m)
vfloor= 10*10^3; % width of valley floor;
vflhalf=vfloor./2;
dxh=vflhalf./dx;

hgt=zeros(xdim,1);

for i=1:xdim;
    for j=1;
        if i< (xc-dxh)
            hgt(i,j)=(abs((xc-dxh)-i).*dx).*tand(alpha);
            if hgt(i,j)>h;
                hgt(i,j)=h;
            end
        elseif (i<= xc+dxh) && (i>= xc-dxh)
            hgt(i,j)=0;
        elseif i> xc+dxh
            hgt(i,j)=(abs((xc+dxh)-i).*dx).*tand(alpha);
            if hgt(i,j)>h;
                hgt(i,j)=h;
            end
        end
    end
end

idx=find((round(hgt/10)).*10==250);

%% attempt to construct this terrain with smooth transitions

h_amp=250;% 1/2 amplitude of the elevation change 

wavenum=4*(2*pi)./(xmax/2); % one wave over the depth of the column
xpoints1=1:dx:(xmax/2);
terrain1=h_amp-h_amp.*(tanh(wavenum.*(xpoints1-(xpoints1(idx(1))))));
xpoints2=(xmax/2):dx:xmax;


terrain2=h_amp+h_amp.*(tanh(wavenum.*(xpoints2-(xpoints2(idx(2)-(length(xpoints1)-1))))));
xpoints=[xpoints1 xpoints2(2:end)];
terrain=[terrain1 terrain2(2:end)];

figure(20); 
movegui('northwest')
plot(xpoints,terrain,'r','linewidth',2);
hold on
plot(1:dx:xmax,hgt);
set(gca,'xtick',[0:5000:(xdim.*dx)]);
xtick=get(gca,'xtick');
xlab=xtick./10^3;
set(gca,'xticklabel',{xlab},'xgrid','on','ygrid','on','layer','top');

%%
hgt_mat=repmat(hgt,[1 ydim]);
terrain_mat=repmat(terrain',[1 ydim]);
%%
f=figure(1);
set(f,'color','w','pos',scrz.*.4)
movegui('northwest')
subplot(1,2,1)
%ha=area(1:dx:xmax,hgt_mat(:,5));
ha=area(1:dx:xmax,terrain_mat(:,5));

set(ha,'facecolor',[.7 .7 .7],'linewidth',2);
set(gca,'xtick',[0:5000:(xdim.*dx)]);
xtick=get(gca,'xtick');
xlab=xtick./10^3;
set(gca,'xticklabel',{xlab},'xgrid','on','ygrid','on','layer','top');
ylim([0 1500]);
set(gca,'ytick',0:250:1500);
subplot(1,2,2)

%slope=gradient(squeeze(hgt_mat(:,5)),dx);
slope=gradient(squeeze(terrain_mat(:,5)),dx);
slpang=atand(abs(slope));
subplot(1,2,2);
plot((1:xdim).*dx,slpang,'-k','linewidth',3);
%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]);


%%
save('ZHONG_TERRAIN.mat','terrain_mat')

%% now reproduce Zhong sounding data

d=50; % 50 meter capping inversion depth
thetatop=300;
thetacil=295;
thetabot=293;

dz=10; % ten meter grid
capdepth=(0:dz:(h-d));
capprofile=linspace(thetabot,thetacil,length(capdepth));
cildepth=((h-d):dz:h);
cilprofile=linspace(thetacil,thetatop,length(cildepth));
aloftdepth=(h+dz):dz:6*h;
aloftprofile=ones(1,length(aloftdepth)).*thetatop;
profile=[capprofile(1:end-1) cilprofile aloftprofile];
alt=[capdepth(1:end-1) cildepth aloftdepth]

figure(2);
clf
plot(profile,alt,'linewidth',2);
grid on;
xlim([thetabot-2 thetatop+2]);

v=zeros(size(profile));
v(1:((length(capdepth)-1)+length(cildepth)))=0;
v(length(capdepth)+length(cildepth):end)=15;
figure(3);
clf
plot(v,alt,'linewidth',2);
grid on;
xlim([-2 17]);


mat=[alt' profile' zeros(size(v))' v' zeros(size(v))'];
%%
fname='zhong_input_sounding';
hdr=[1000 profile(1) 0];
dlmwrite(fname,hdr,'delimiter','\t','precision', '%.2f');
dlmwrite(fname,mat,'-append','delimiter','\t','precision', '%.2f');

