%load sonde matlab file a do a few simple plots.
% ddds=input('Input yearday ddd ','s')
load(['D:\DATA\ATOMIC\P3\sonde\data\all\ATOMIC_P3_sonde_all.mat']);
[i mon dttn hr mn sec]=datevec(jdxh);

thet=(tahh+273.16).*(1010./Phh).^0.287-273.16;%air theta
rhoair=Phh.*100./(287.1.*(tahh+273.16).*(1+0.61*qahh/1e3));
tdew=[];[mn mm]=size(jdxh);
for i=1:mm;
    tdew(:,i)=tdewrh([Phh(:,i) tahh(:,i) rhahh(:,i)]);
end;
    
jk=1;figure;plot(tahh(:,jk), zhh(:,jk),tdew(:,jk),zhh(:,jk),thet(:,jk),zhh(:,jk),'.');xlabel('Temperature (C)');ylabel('Altitude (m)');
legend('Tair','Tdew','Tpotential');
LCLbias=nanmeanx(10+(tahh(2,:)-tdew(2,:))/.008 - (100+(tahh(20,:)-tdew(20,:))/.008));
figure;plot(lonhh,10+(tahh(2,:)-tdew(2,:))/.008,'x',lonhh,100+(tahh(20,:)-tdew(20,:))/.008,'o');xlabel('Longitude (deg)');ylabel('LCL (m)');
xlim([-60 -49]);grid;hold on;
x=lonhh;y=10+(tahh(2,:)-tdew(2,:))/.008;
[X, Y, M, SD, SE, SDM, SEM, N] = binave_3(x,y,1,-60,-48);
x=lonhh;y=100+(tahh(20,:)-tdew(20,:))/.008;
[Xa, Ya, Ma, SDa, SE, SDM, SEM, N] = binave_3(x,y,1,-60,-48);
plot(X,M,'-x',X,Ma,'-o');
legend('10-m Met','100-m Met','10-m median','100-m mdeian');title('ATOLMIC P3 Dropsondes LCL');

jk=8:10;figure;
subplot(1,2,1);semilogy(thet(:,jk),zhh(:,jk),'-o');ylim([0 200]);
subplot(1,2,2);semilogy(qahh(:,jk),zhh(:,jk),'-o');ylim([0 200]);


Nx=[];zhmn=[];qamn=[];Umn=[];rhmn=[];Pmn=[];tamn=[];thetmn=[];Ni=[];sqamn=[];sUmn=[];stamn=[];
Nt=[];zhmntp=[];qamntp=[];Umntp=[];rhmntp=[];Pmntp=[];tamntp=[];
lonmn=-58;lonmx=-50;
for i=1:2000;
    ik=find(~isnan(rhahh(i,:)) & lonhh>=lonmn & lonhh<lonmx);
    Nx(i)=length(ik);
    if Nx(i)>2;
        zhmn(i)=nanmedian(zhh(i,ik));
        qamn(i)=nanmedianx(qahh(i,ik)); sqamn(i)=nanstdx(qahh(i,ik));
        Umn(i)=nanmedianx(Uhh(i,ik));sUmn(i)=nanstdx(Uhh(i,ik));
        tamn(i)=nanmeanx(tahh(i,ik));  stamn(i)=nanstdx(tahh(i,ik));      
        thetmn(i)=nanmedianx(thet(i,ik));
        rhmn(i)=nanmeanx(rhahh(i,ik));
        Pmn(i)=nanmeanx(Phh(i,ik));    
        Nx(i)=length(ik);
    end;
%     ik=find(~isnan(rhahh(i,:)) & lathh<latline);
%     Ny(i)=length(ik);
%     if Ny(i)>2;
%         zhmntp(i)=nanmedian(zhh(i,ik));
%         qamntp(i)=nanmeanx(qahh(i,ik));
%         Umntp(i)=nanmeanx(Uhh(i,ik));
%         tamntp(i)=nanmeanx(tahh(i,ik));
%         rhmntp(i)=nanmeanx(rhahh(i,ik));
%         Pmntp(i)=nanmeanx(Phh(i,ik));                
%     end;
%    
%     rhmx(i)=max(rhahh(i,:));

end;%

figure;
subplot(1,2,1);semilogy(thetmn,zhmn,'-o');ylim([0 200]);grid
subplot(1,2,2);semilogy(qamn,zhmn,'-o');ylim([0 200]);grid

figure;plot(stamn,zhmn,'-o',sqamn,zhmn,'-o',sUmn,zhmn,'-o');;grid;legend('sig  ta','sig qa','sig U');xlim([0 6]);ylim([0 7500])

x=lonhh;y=tahh(2,:);
[X, Y, M, SD, SE, SDM, SEM, N] = binave_3(x,y,1,-60,-48);
x=lonhh;y=tdew(2,:);
[Xa, Ya, Ma, SDa, SE, SDM, SEM, N] = binave_3(x,y,1,-60,-48);
figure;plot(X,M,'-x',X,Ma,'-o');xlabel('Longitude (deg)');ylabel('Temperature (C)');grid;legend('Ta 10-m (C)','Td 10-m (C)');

figure;plot(lonhh,lathh,'o');xlabel('Longitude (deg)');ylabel('Latitude (deg)');title('ATOMIC 2020 P3 Dropsonde locations');


load('D:\DATA\ATOMIC\P3\FlightLevel_summary\ATOMIC_P3_FLTLVL_all.mat');
y=P3_Tah;x=P3_altrefh;
[X3, Y3, M3, SD3, SE3, SDM3, SEM3, N3] = binave_3(x,y,100,100,7500);
figure;plot(nanmeanx(tahh'),nanmeanx(zhh'),Y3,X3,'o');xlabel('Air Temperature (C)');ylabel('Altitude (m)');legend('Dropsonde','P3 Flight Level');grid;
%q=qair3([P3_pressh P3_Tah P3_RHh]);y=q;
y=P3_Qh;ii=find(P3_lonh>-60);
[X3, Y3, M3, SD3, SE3, SDM3, SEM3, N3] = binave_3(x(ii),y(ii),100,100,7500);
figure;plot(nanmeanx(qahh'),nanmeanx(zhh'),Y3,X3,'o');xlabel('Specific Humidity (g/kg)');ylabel('Altitude (m)');legend('Dropsonde','P3 Flight Level');grid;
figure;plot([0 diff(nanmeanx(thet'))],nanmeanx(zhh'),'.');xlabel('d\Theta /dz (C/m)');ylabel('Altitude (m)');grid;xlim([-.0 .075])
dthdz=diff(nanmeanx(thet'));dPdz=diff(nanmeanx(Phh'));
figure;plot(-[0 dthdz/9.83]./nanmeanx(rhoair'),nanmeanx(Phh'),'.');xlabel('d\Theta /dp (C/mb)');ylabel('Pressure (mb)');set(gca,'YDir','reverse');grid;xlim([-1 .1]*1e-2)
%%%%%%%%  set height for 2nd means 
zu=18;
zt=18;
zq=18;
zrf_u=10;
zrf_t=10;
zrf_q=10;

%%%%%%%%%%%%%%

u=Uhh(3,:);
t=tahh(3,:);
rh=rhahh(3,:);
P=Phh(3,:);
Rs=500;
Rl=400;
lat=lathh;
rain=0;
Ss=35;
zi=750;
ts=27.35-0.22*lathh-0.05*lonhh;
Tf=-0.0575*Ss+1.71052E-3*Ss.^1.5-2.154996E-4*Ss.*Ss;%water freezing point (C)
cpm=NaN*ones(1,length(ts));
sigH=cpm;
hsgm=cpm;
Wa=NaN;%cpm./U;%wave age (Cp/U10)

%%%%%%%%   call bulk algo
%%%%%%    if ts<Tf assumes ice
A=coare36vn_zrf(u,zu,t,zt,rh,zq,P,ts,Rs,Rl,lat,zi,rain,Ss,cpm,hsgm,zrf_u,zrf_t,zrf_q);
usrw=A(:,1);% u*
taub=A(:,2);
hsb=A(:,3);%sensible heat
hlb=A(:,4);%latent heat
zo=A(:,10);
cdz=A(:,13);
chz=A(:,14);
dter=A(:,18);%cool skin
RHrf=A(:,24);RH10=A(:,40);%Rh at 2nd height and 10m
%u10n=A(:,31);
cdn_10=A(:,32);%neutral 10m drag coeff
chn_10=A(:,33);%neutral 10m heat xfer c
u10n=A(:,31);%10-m neutral wind speed
Whf=A(:,42);%white cap fraction
Edis=A(:,43);%wave breaking energy dissipation

load('Z:\ATOMIC\RHB\flux\Field_Processed\flux_met_sea\ATOMIC_2020_10min_flux_met_sea_data');

figure;
subplot(3,1,1);plot(c10.jd,c10.tau,jdxh,taub,'o');ylabel('Stress (Nt/m^2)');title('ATOMIC Turbulent surface fluxes ');grid;legend('\tau RHB','\tau Dropsonde');
subplot(3,1,2);plot(c10.jd,c10.hs,jdxh,hsb,'o');ylabel('Sensible Heat Flux (W/m^2)');;grid;legend('Hs RHB','Hs Dropsonde');
subplot(3,1,3);plot(c10.jd,c10.hl,jdxh,hlb,'o');ylabel('Latent Heat Flux (Nt/m^2)');xlabel('Year Day (2020)');;grid;legend('Hl RHB','Hl Dropsonde');

figure;plot(jdxh,rh,'o',c10.jd,c10.rh);

function y=tdewrh(y)
p=y(:,1);%pressure in mb
t=y(:,2);%Tair (C)
rh=y(:,3);%rel hum (%)
e=rh/100*6.112.*exp(17.502.*t./(t+241.0)).*(1.0007+3.46e-6*p);%vapor pressure, from Buck
z=log(e./(1.0007+3.46e-6*p)/6.112);
y=241.*z./(17.502-z);
end
