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

thet=(tahh+273.16).*(1020./Phh).^0.287-273.16;%air theta
rhoair=Phh.*100./(287.1.*(tahh+273.16).*(1+0.61*qahh/1e3));
tdew=[];[nn 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('E:\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])

dth=diff(nanmeanx(thet'));dP=diff(nanmeanx(Phh'));dThetdp=[0 dth./dP];
dz=diff(nanmeanx(zhh'));dThetdz=[0 dth./dz];
rhoairm=nanmeanx(rhoair');Phhm=nanmeanx(Phh');
drhodp=[0 diff(rhoairm)./dz];

figure;plot(dThetdp,nanmeanx(Phh'),'.');xlabel('d\Theta /dp (C/mb)');ylabel('Pressure (mb)');set(gca,'YDir','reverse');grid;xlim([-1 .1]*1e-1)
figure;plot(dThetdp,nanmeanx(zhh'),'.');xlabel('d\Theta /dp (C/mb)');ylabel('z (m)');grid;xlim([-1 .1]*1e-1)

rhov=rhoair.*qahh;
z=nanmeanx(zhh');pwv=[];tops=[];
for i=1:300
    dz1=diff(zhh(:,i));dz1=[dz1' 0]';
    pwv(i)=nansumx(dz1.*rhov(:,i)/1e4);
    j=find(isnan(zhh(:,i)));
    if length(j)<4000
        tops(i)=zhh(j(1)-1,i);
    else
        tops(i)=NaN;
    end;
end;
kk=find(tops<6000);pwv(kk)=pwv(kk)+.14;
kk=find(pwv<2);pwv(kk)=NaN;

%%%%%%%%  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,:);
qa=qair3([P' t' rh'])';
iv=pwv;
jdy=jdxh';lon=lonhh';lat=lathh';iv=iv';
ksol=.01*ones(length(lon),1);
Atomic_Rsclear_5_noplot

%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)
Rs=sw';
rlclr=(.52+.13/60*abs(lathh)+(.082-.03/60*abs(lathh)).*sqrt(qa)).*(5.67e-8*(t+273.15).^4); 
hsig=(0.02*(1.1).^1.1-0.0025).*u.^2;
Rl=rlclr;

%%%%%%%%   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

%save('D:\DATA\ATOMIC\P3\sonde\data\all\ATOMIC_P3_drop_surfaceFlux.mat', 'jdxh','lathh','lonhh','u','t','ts','rh','qa','Rs','Rl','P','hsig','taub','hsb','hlb','dter','u10n','Whf');


%load('Z:\ATOMIC_2020\RHB\flux\Field_Processed\flux_met_sea\ATOMIC_2020_10min_flux_met_sea_data');
load('E:\DATA\ATOMIC\RBH\Brown\flux\Field_Processed\flux_met_sea\ATOMIC_2020_10min_flux_met_sea_data_v1');

Rlclr=(.53+.13/60*abs(c10.lat)+(.084-.03/60*abs(c10.lat)).*sqrt(c10.qa)).*(5.67e-8*(c10.Ta+273.15).^4); 
iv=NaN*c10.qa;qa=c10.qa;
jdy=c10.jd;lon=c10.lon;lat=c10.lat;iv=iv;
ksol=.01*ones(length(lon),1);
mm=find(jdy>30 & jdy<31);ksol(mm)=0.055;
mm=find(jdy>=31 & jdy<33);ksol(mm)=0.15;
mm=find(jdy>=33 & jdy<38 );ksol(mm)=0.07;
mm=find(jdy>=38 );ksol(mm)=0.05;
%ksol=0.15;
Atomic_Rsclear_5_noplot
%amma06_Rsclear_5_noplot

sw_net_cl=Rscl*.95;
lw_net_cl=0.97*(Rlclr-5.67e-8*(c10.Tskin+273.16).^4);
figure;plot(jdy,sw,'o',jdy,-c10.sw_dn,'-x');j=1;xlim([j+.4 j+.9]);ylim([ 00 1000]);xlabel('Day Number (2020)');ylabel('Solar Flux (W/m^2)');legend('Clear sky','RHB observed');

ceilo=load('E:\DATA\ATOMIC\RBH\Ceilo\Ceilo_10min_ATOMIC2020.txt');
ceiloh=load('E:\DATA\ATOMIC\RBH\Ceilo\Ceilo_hourly_ATOMIC2020.txt');

cloudf=ceilo(:,7);
cloudbase=ceilo(:,13);
ii=find(c10.sw_dn<2);x=cloudf(ii);y=-c10.lw_net(ii)-lw_net_cl(ii)-5;
[Xc, Yc, Mc, SDc, SE, SDM, SEM, N] = binave_3(x,y,.1,0,1);
ii=find(c10.sw_dn<-100);x=cloudf(ii);y=-c10.sw_net(ii)-sw_net_cl(ii);
[Xs, Ys, Ms, SDs, SE, SDM, SEM, N] = binave_3(x,y,.1,0,1);
figure;plot(Xc,Mc,'o',Xs,Ms+10,'o',[0 1],[0 0],[0 1],[0 30],'--',[0 1],[0 -150],'--');xlabel('Cloud fraction');ylabel('Surface Cloud Radiative Forcing (W/m^2)');grid
legend('Longwave','Solar');title('10-min ATOMIC Cloud Radiative Forcing Ronald H. Brown');

jdcloudh=ceiloh(:,1);
cloudfh=ceiloh(:,7);
cloudbaseh=ceiloh(:,13);

x=24*(jdcloudh-floor(jdcloudh));y=cloudfh;
[Xs, Ys, Ms, SDs, SE, SDM, SEM, N] = binave_3(x,y,2,0,24);
figure;plot(Xs,Ys,'-o');
x=24*(jdcloudh-floor(jdcloudh));y=cloudbaseh;
[Xs, Ys, Ms, SDs, SE, SDM, SEM, N] = binave_3(x,y,2,0,24);
figure;plot(Xs,Ys);



aod=load('C:\Data\cwf\Dropbox\division\Review_2015\NOAA_PSD_init\elnino\EUREC4A\post_field\Ron_Brown_AOD.txt');
figure;plot(aod(:,1),aod(:,5),'x',aod(:,1),aod(:,7),'o',jdy,ksol);xlabel('Year day (2020)');ylabel('Aerosol');legend('Microtops AOD 380','Microtops AOD 500','Iqbal k');

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;
subplot(3,1,1);plot(c10.jd,c10.Ta,jdxh,t,'o');ylabel('18-m Air Temperature (C)');title('ATOMIC Bulk Meteorology ');grid;legend('ta RHB','ta Dropsonde');
subplot(3,1,2);plot(c10.jd,c10.rh,jdxh,rh,'o');ylabel('18-m Relative Humidity (%)');grid;legend('RH RHB','RH Dropsonde');
subplot(3,1,3);plot(c10.jd,c10.wspd,jdxh,u,'o');ylabel('18-m Wind Speed (m/s)');xlabel('Year Day (2020)');grid;legend('U RHB','U Dropsonde');

[nn mm]=size(jdxh);pmin=390;pmax=1020;dpp=10;nx=(pmax-pmin)/dpp+1;
tahhp=ones(nx,mm)*NaN;pahhp=tahhp;zhhp=tahhp;qahhp=tahhp;uxhhp=tahhp;uyhhp=tahhp;rhoahhp=tahhp;p8=[];z8=[];p5=[];z5=[];nhhp=tahhp;
for i=1:mm
    [X, Y, M, SD, SE, SDM, SEM, N] = binave_4(Phh(:,i),tahh(:,i),dpp,pmin,pmax);
    ny=length(Y);
    tahhp(1:ny,i)=Y;
    pahhp(1:ny,i)=X;
   [X, Y, M, SD, SE, SDM, SEM, N] = binave_4(Phh(:,i),zhh(:,i),dpp,pmin,pmax);
    ny=length(Y);
    zhhp(1:ny,i)=Y;
   [X, Y, M, SD, SE, SDM, SEM, N] = binave_4(Phh(:,i),qahh(:,i),dpp,pmin,pmax);
    ny=length(Y);
    qahhp(1:ny,i)=Y;
   [X, Y, M, SD, SE, SDM, SEM, N] = binave_4(Phh(:,i),Uhh(:,i).*cos(udirhh(:,i)*pi/180),dpp,pmin,pmax);
    ny=length(Y);
    uyhhp(1:ny,i)=-Y;
   [X, Y, M, SD, SE, SDM, SEM, N] = binave_4(Phh(:,i),Uhh(:,i).*sin(udirhh(:,i)*pi/180),dpp,pmin,pmax);
    ny=length(Y);
    uxhhp(1:ny,i)=-Y;
    [X, Y, M, SD, SE, SDM, SEM, N] = binave_4(Phh(:,i),rhoair(:,i),dpp,pmin,pmax);
    ny=length(Y);
    rhoahhp(1:ny,i)=Y;
    nhhp(1:ny,i)=N;
    ii=find(Phh(:,i)>800 & Phh(:,i)<805);
    if ~isempty(ii)
        p8(i)=Phh(ii(length(ii)),i);z8(i)=zhh(ii(length(ii)),i);
    else
        p8(i)=NaN;z8(i)=NaN;
    end;
   ii=find(Phh(:,i)>500 & Phh(:,i)<510);
    if ~isempty(ii)
        p5(i)=Phh(ii(length(ii)),i);z5(i)=zhh(ii(length(ii)),i);
    else
        p5(i)=NaN;z5(i)=NaN;
    end;
    
end;
ppp=pmin:dpp:pmax;
dP=[];
Phihhp=ones(nx,mm)*0;R=287;dP=-[dpp diff(ppp)];
for i=length(ppp)-1:-1:1 
    Phihhp(i,:)=Phihhp(i+1,:)-R*(tahhp(i,:)+273.15)./(ppp(i)+2.5).*dP(i);
end;
    

jj=find(lonhh>-55 & ~isnan(z5));kk=find(lonhh<=-55 & ~isnan(z5));figure;plot(nanmeanx(tahhp(:,jj)'),ppp,nanmeanx(tahhp(:,kk)'),ppp);
jj=find(lonhh>-55 & ~isnan(z5));kk=find(lonhh<=-55 & ~isnan(z5));figure;plot(nanmeanx(qahhp(:,jj)'),ppp,nanmeanx(qahhp(:,kk)'),ppp);
jj=find(lonhh>-55 & ~isnan(z5));kk=find(lonhh<=-55 & ~isnan(z5));figure;plot(nanmeanx(uxhhp(:,jj)'),ppp,nanmeanx(uxhhp(:,kk)'),ppp);
jj=find(lonhh>-55 & ~isnan(z5));kk=find(lonhh<=-55 & ~isnan(z5));figure;plot(nanmeanx(uyhhp(:,jj)'),ppp,nanmeanx(uyhhp(:,kk)'),ppp);
dux_dx=(nanmeanx(uxhhp(:,jj)')-nanmeanx(uxhhp(:,kk)'))/330e3;
jj=find(lathh>12.5 & ~isnan(z5));kk=find(lathh<=12.5 & ~isnan(z5));figure;plot(nanmeanx(uxhhp(:,jj)'),ppp,nanmeanx(uxhhp(:,kk)'),ppp);
duy_dy=(nanmeanx(uxhhp(:,jj)')-nanmeanx(uxhhp(:,kk)'))/185e3;
gg=find(~isnan(z5));
rhoairmx=nanmeanx(rhoahhp(:,gg)');
zzm=nanmeanx(zhhp(:,gg)');

DIV=dux_dx+duy_dy;;dP=-[dpp diff(ppp)];dz=[90 -diff(nanmeanx(zhhp(:,jj)'))];
ii=find(isnan(DIV));DIV(ii)=0;
om=DIV*0;W=DIV*0;Wm=W;
for i=length(ppp)-1:-1:1 
    om(i)=om(i+1)-DIV(i)*dP(i);%Pa*s^-1
    W(i)=W(i+1)-DIV(i)*dz(i);
    Wm(i)=Wm(i+1)-rhoairmx(i)*DIV(i)*dz(i);
end;
    
figure;plot(om*3600,ppp,-W*.01*3600.*rhoairmx*9.88,ppp,'.',-Wm*.01*3600.*9.88,ppp,'x')
    


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
function y=qairv3(y)
p=y(:,1);%pressure in mb
t=y(:,2);%temperature in C
h=y(:,3);%relative humidity in %
es=6.112.*exp(17.502.*t./(t+241.0)).*(1.0007+3.46e-6*p).*h/100;
q=es*622./(p-.378*es*h/100);

tdk=273.16;
Rgas=287.1;
rhoa=p*100/(Rgas*(t+tdk)*(1+0.61*q/1000));
rhov=rhoa*q;% (g/m^3)
y=q;
end