%% CREATE WRF SOUNDING FOR WRF IDEALIZED CASES

% Folowing Strang and Fernando (2001) JFM, we need to define a bulk
% richardson number, Rb= ((delta buoancy).* D)./ (delta U)
% for a stratified sheared interface with an upper mixed layer. The lower
% layer may be nuetrally stratified of positively stratified

%%%% DEFINITIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% buoyancy = g*(Ttheta - Theta_mean)./Theta_mean
% delta buoyancy b= g(Theta_upper - Theta_lower)./Theta_mean
% delta u= U_upper - U_lower, where U_lower = 0 --->  delta u == U_upper
% D = depth of the upper mixed layer
% H= total depth
% ds= shear layer depth
% db = buoyancy jump depth
% d = displacement of the center of the two tranisiton layers


H=3000; % meters
D=4.*H./5; % half way
U_upper= 10;% m/s
U_lower= 0;
T_upper = 305;
T_lower = 295;
T_mean = 300;
g=9.81;
ds=200; % shear layer depth

delta_b= g.*(T_upper-T_lower)./T_mean;
rib=(delta_b.*D)./U_upper.^2; % bulk Richardson following Strang and Fernando
disp(strcat({'The bulk Richardson number is: '},num2str(rib)));

ris=delta_b.*ds./(U_upper.^2); % shear layer Richardson number;
disp(strcat({'The shear layer Richardson number is: '},num2str(ris)));


%% now construct a profile for model initialization

hgts=0:10:H;
hanom=hgts-(H/5);%mean(hgts);
wavenum=(2*pi)./max(hgts); % one wave over the depth of the column
t_amp=5; % 10 K inversion layer (note total amp = 2 times amp)
t_base= 300; % 300 K mean temperature
u_amp= 5; % 10 m/s velocity jump
u_base=5; % no wind within the CAP;
d_offset=0; % 10 meter offset of layer thickness;
temp=t_base+t_amp.*(tanh(4*wavenum.*hanom));
dthetadz=gradient(temp,10);
u=zeros(1,length(temp));
u(hgts>(mean(hgts)+125))=u_amp*2;

%u=u_base+u_amp.*(tanh(4*wavenum.*(hanom + d_offset)));

figure(1); clf
subplot(1,2,1)
plot(temp,hgts);
xlim([min(temp)-2 max(temp)+2]);
subplot(1,2,2);
plot(u,hgts);
xlim([min(u)-2 max(u)+2]);



%sounding_mat=[round(hgts)'  temp' zeros(size(u')) u' zeros(size(u'))];
%
%fname='input_sounding_strang';
%hdr=[850 temp(1) 0];
%dlmwrite(fname,hdr,'delimiter','\t','precision', '%.2f');
%dlmwrite(fname,sounding_mat,'-append','delimiter','\t','precision', '%.2f');



%% 
% script to read PCAPS soundings and export WRF input soundings for
% idealized cases

close all
clear all

% sounding data location
file='/uufs/chpc.utah.edu/common/home/metwx/pcaps.utah.edu/data/SONDE_PROC/20101203/KSLC_2010120300_smth.dat';

inpt=importdata(file,' ',6); %skip 7 header lines
data(:,:)=inpt.data;
PRES=data(:,2);
PRES(PRES==-999)=nan;
THth=data(:,9);
THth(THth==-999)=nan;%mean(THth(:));
DPth=data(:,4);
DPth=nan;
RHth=data(:,5);
RHth(RHth==-999)=nan;%mean(RHth(:));
u=data(:,14);
u(abs(u)>50)=nan;
v=data(:,15);
v(abs(v)>50)=nan;
hgts=data(:,1);

idx=find(hgts<5000);
hgts=hgts-hgts(1);
hgts=hgts(idx);
u=u(idx);
v=v(idx);
THth=THth(idx);


figure;
subplot(1,2,1)
plot(THth,hgts);
%ylim([0 2500])
xlim([280 310])
subplot(1,2,2);
plot(v,hgts);
hold on;
plot(u,hgts,'r');
%ylim([0 2500]);


mat=[round(hgts) THth zeros(size(u)) v zeros(size(u))];
%%
fname='input_sounding';
hdr=[1000 THth(1) 0];
dlmwrite(fname,hdr,'delimiter','\t','precision', '%.2f');
dlmwrite(fname,mat,'-append','delimiter','\t','precision', '%.2f');




