disp('read_sondeProfile_calwater2');
%This program reads stored dropsonde profiles from each of 4 platforms
%The file is a matlab .mat file.

%select the aircraft by uncommenting the one you want
%flt='N1';%G4
% flt='I1';%P3
 %flt='GH';%global hawk
 %flt='H1';%P3
 way_raw_files='C:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\2015\sonde\data\';
 way_raw_filesG='C:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\2015\sonde\data\';%'Z:\CALWATER_2015\G4\data\';
 way_raw_filesP='C:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\2015\sonde\data\';%'Z:\CALWATER_2015\P3\data\';
 way_raw_filesA='C:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\2015\axbt\';%'Z:\CALWATER_2015\P3\axbt\';

 

sondex=load([way_raw_filesG 'CALWATER_G4_sonde_all.mat'], 'jdxh', 'zhh','Uhh','rhahh','tahh','lathh','lonhh','Phh','qahh','udirhh');
sondey=load([way_raw_filesP 'CALWATER_P3_sonde_all.mat'], 'jdxh', 'zhh','Uhh','rhahh','tahh','lathh','lonhh','Phh','qahh','udirhh');
axbt=load([way_raw_filesA '\CALWATER2_SST_data.mat']);

rhahhG=sondex.rhahh;
zhhG=sondex.zhh;
UhhG=sondex.Uhh;
tahhG=sondex.tahh;
lathhG=sondex.lathh;
lonhhG=sondex.lonhh;
PhhG=sondex.Phh;
qahhG=sondex.qahh;
udirhhG=sondex.udirhh;
jdxhG=sondex.jdxh;
jk=length(lathhG);
[i mon dttn hr mn sec]=datevec(jdxhG);

[nn nmg]=size(zhhG);
disp(['Number of profiles= ' num2str(nmg)]) 

rhahhP=sondey.rhahh;
zhhP=sondey.zhh;
UhhP=sondey.Uhh;
tahhP=sondey.tahh;
lathhP=sondey.lathh;
lonhhP=sondey.lonhh;
PhhP=sondey.Phh;
qahhP=sondey.qahh;
udirhhP=sondey.udirhh;
jdxhP=sondey.jdxh;
jkP=length(lathhP);
[i monp dttnp hrp mnp secp]=datevec(jdxhP);


tahhA=axbt.tazh;
lathhA=axbt.latzh;
lonhhA=axbt.lonzh;
jdxhA=axbt.jdzh;
jkA=length(lathhA);
[i mona dttna hra mna seca]=datevec(jdxhA);



x=load('C:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\2015\Calwater-2015_flux_10_v4_decorr.txt');%read file with 10-min average data; set your local directory. The columns assignment is as follows:
jdy=x(:,1);%day-of-year at beginning of time average
ushp=x(:,2);%doppler log, SCS (m/s) ?not available, replaced by SOG
U=x(:,3);%true wind PSD sonic relative to water (m/s)
dir=x(:,4);%true wind direction from relative to water, PSD sonic (deg)
urel=x(:,5);%relative wind speed, PSD (m/s)
reldir=x(:,6);%relative wind dir (from),clockwise rel ship's bow, PSD sonic (deg)
head=x(:,7);%ship heading, deg clockwise rel north, PSD GPS
tsnk=x(:,8);%sea snake temperature, PSD (C)
ta=x(:,9);%air temperature, PSD (C)%
qse=x(:,10);%sea surface specific humidity from snake (g/kg)
qa=x(:,11);%air specific humidity, PSD (g/kg)
hsc=x(:,12);%sensible heat flux, covariance, PSD sonic anemometer(W/m^2)
hsib=x(:,13);%sensible heat flux, ID, PSD sonic anemometer(W/m^2)
hsb=x(:,14);%bulk sensible heat flux, (W/m^2)
hlc=x(:,15);%latent heat flux, covariance, (W/m^2)
hlib=x(:,16);%latent heat flux, ID, (W/m^2)
hlb=x(:,17);%bulk latent heat flux, W/m^2 (includes Webb et al. correction)
taucx=x(:,18);%covariance streamwise stress, PSD sonic anemometer (N/m^2)
taucy=x(:,19);%covariance cross-stream stress, PSD sonic anemometer (N/m^2)
tauib=x(:,20);%ID streamwise stress, PSD sonic anemometer (N/m^2)
taub=x(:,21);%bulk wind stress along mean wind, (N/m^2)
rs=x(:,22);%downward solar flux, PSD units (W/m^2)
rl=x(:,23);%downward IR flux, PSD units (W/m^2)
org=x(:,24);%rainrate, PSD STI optical rain gauge (mm/hr)
J=x(:,25);%ship plume contamination index
tiltx=x(:,26);%flow tilt at PSD sonic anemometer
Jm=x(:,27);%ship maneuver index
ct=x(:,28);%ct^2 (K^2/m^.667)
cq=x(:,29);%cq^2 ((g/kg)^2/m^.667)
cu=x(:,30);%cu^2 ((m/s)^2/m^.667)
cw=x(:,31);%cw^2 ((m/s)^2/m^.667)
hrain=x(:,32);%rain heat flux, Gosnell et al 1995, JGR, 18437-18442 (W/m^2)
hlwebb=x(:,33);%correction to measured latent heat flux, Webb et al.
lat=x(:,34);%latitude, deg (PSD pcode)
lon=x(:,35);%longitude, deg (PSD pcode)
zu_etl=x(:,36);%height of mean wind sensor, 17.75 m
zt_etl=x(:,37);%height of mean air temperature sensor, 15.6 m
zq_etl=x(:,38);%height of mean air humidity sensor, 15.6 m
%***** ships imet and scs data
sog=x(:,39);%speed over ground, SCS gps, (m/s)
U_scs=x(:,40); %true wind speed relative to earth m/s) – 19.6m
dir_scs=x(:,41);%true wind direction from relative to earth, clockwise rel north (deg)
cog=x(:,42);%%course over ground, SCS gps, (m/s)
tsg=x(:,43);%tsg water temperature (C)
ta_im=x(:,44);%imet air temperature (C)
qs_tsg=x(:,45);%imet bulk water specific humidity (g/kg)
qa_im=x(:,46);%imet air specific humidity, (g/kg)
rs_im=x(:,47);%imet solar flux, (W/m^2)
rl_im=x(:,48);%imet IR flux (W/m^2)
wco2_lic=x(:,49);%LICOR 7500 CO2 flux, (micatm m/s)
q_lic=x(:,50);%Specific humidity from LICOR (g/kg)
sgq_lic=x(:,51);%Standard deviation of specific humidity from LICOR (g/kg)
co2_lic=x(:,52); %CO2 concentration from Licor (umol/mol)
sgC_lic=x(:,53);%Standard deviation of CO2 concentration from LICOR (microatm)
press=x(:,54); %Atmospheric pressure (mb)
Uearth=x(:,55); %True wind speed (m/s) relative to earth
direarth=x(:,56); % True wind direction (deg) from relative to earth%  ii=find(lat>100);lat(ii)=lat(ii-1);
% ii=find(lon>-70);lon(ii)=lon(ii+1);
rh=relhum5([ta qa press ]);
ii=find(jdy>36.65);lon(ii)=232.823;
zu_etl=ones(length(jdy),1)*18.5;
zt_etl=ones(length(jdy),1)*15.5;
zq_etl=ones(length(jdy),1)*15.5;

%load('C:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\2015\Calwater_2015_PSD_flux_5min_allcruise3.mat');%read file with 10-min average data; set your local directory. The columns assignment is as follows:

       [n m]=size(x);
%    yr=ones(n,1)*2010;
   org2=org;ii=find(jdy>305);org2(ii)=org(ii)-min(org(ii));org=org2;
   xco2=co2_lic;
   lons=lon;
   lats=lat;
   q_oph=NaN;
   ophm=NaN;
   sigoph=NaN;
   
   ffc=(1.58+1.95*log(qa))/10;
   ii=find(~isnan(qa));
   hlc(ii)=hlc(ii)+(qa(ii)-q_lic(ii))./max(1,qa(ii)).*ffc(ii).*hlc(ii);
   hlib(ii)=hlib(ii)+(qa(ii)-q_lic(ii))./max(1,qa(ii)).*ffc(ii).*hlib(ii);
   
   hlib_lic=hlib;

%select a profile to plot
jj=1;
ii=find(zhhG(:,jj)>0);size(ii);
disp(['Max Altitude= ' num2str(zhhG(length(ii),jj)) ' m']) 
figure;
subplot(1,2,1);plot(UhhG(:,jj),zhhG(:,jj)/1e3);ylabel('Alt (km)');xlabel('U (m/s)');title(['CALWATER2 2015 '  num2str(mon(jj)) '/' num2str(dttn(jj)) ' '  num2str(hr(jj)+mn(jj)/60)]);ylim([0 5]);
subplot(1,2,2);plot(rhahhG(:,jj),zhhG(:,jj)/1e3,[100 100],[0  15],'--');;xlabel('RH water-ice (%)');xlim([0 105]);ylim([0 5]);

figure;plot(lonhhG,lathhG,'o',lonhhP,lathhP,'x');xlabel('Longitude');ylabel('Latitude');title('CALWATER2 2015 Dropsondes');legend('G4','P3')
figure;plot(jdy,tsnk,jdy,ta,'d',jdxhP,tahhP(3,:),'x',axbt.jdzh,axbt.tazh,'o');

ii=find(jdy>27 & jdy<28.5);jj=find(jdxhP>27 & jdxhP<28.5);kk=find(axbt.jdzh>27 & axbt.jdzh<28.5);
figure;plot(lat(ii),tsnk(ii),'.',lat(ii),ta(ii),'o',lathhP(jj),tahhP(2,jj),'x',axbt.latzh(kk),axbt.tazh(kk),'d');xlim([37.95 38.05])

ii=find(jdy>31.9 & jdy<32.1);jj=find(jdxhP>31.9 & jdxhP<32.1);kk=find(axbt.jdzh>31.9 & axbt.jdzh<32.1);
figure;plot(lon(ii)-360,lat(ii),'o',lonhhP(jj),lathhP(jj),'x',axbt.lonzh(kk),axbt.latzh(kk),'d');xlabel('Longitude');ylabel('Latitude');title('CALWATER2 Jan 31, 2015 P3 Drops');
legend('RHB','Ta','SST');
ij=find(lonhhP(jj)<-124 & lonhhP(jj)>-126);ik=find(axbt.lonzh(kk)<-124 & axbt.lonzh(kk)>-126);
figure;plot(lat(ii),tsnk(ii),'.',lat(ii),ta(ii),'o',lathhP(jj(ij)),tahhP(2,jj(ij)),'x',axbt.latzh(kk(ik)),axbt.tazh(kk(ik)),'d');%xlim([37.95 38.05])

ii=find(jdy>36.7 & jdy<37);jj=find(jdxhP>36.7 & jdxhP<37);kk=find(axbt.jdzh>36.7 & axbt.jdzh<37);
figure;plot(lon(ii)-360,lat(ii),'o',lonhhP(jj),lathhP(jj),'x',axbt.lonzh(kk),axbt.latzh(kk),'d');xlabel('Longitude');ylabel('Latitude');title('CALWATER2 Feb 5, 2015 P3 Drops');
legend('RHB','Ta','SST');
ij=find(lonhhP(jj)<-126 & lonhhP(jj)>-130);ik=find(axbt.lonzh(kk)<-126 & axbt.lonzh(kk)>-130);
figure;plot(lat(ii),tsnk(ii),'.',lat(ii),ta(ii),'o',lathhP(jj(ij)),tahhP(2,jj(ij)),'x',axbt.latzh(kk(ik)),axbt.tazh(kk(ik)),'d');%xlim([37.95 38.05])

ii=find(jdy>37.5 & jdy<39);jj=find(jdxhP>37.5 & jdxhP<39);kk=find(axbt.jdzh>37.5 & axbt.jdzh<39);
figure;plot(lon(ii)-360,lat(ii),'o',lonhhP(jj),lathhP(jj),'x',axbt.lonzh(kk),axbt.latzh(kk),'d');xlabel('Longitude');ylabel('Latitude');title('CALWATER2 Feb 6-7, 2015 P3 Drops');
legend('RHB','Ta','SST');
ij=find(lonhhP(jj)<-124 & lonhhP(jj)>-129);ik=find(axbt.lonzh(kk)<-124 & axbt.lonzh(kk)>-129);
figure;plot(lat(ii),tsnk(ii),'.',lat(ii),ta(ii),'o',lathhP(jj(ij)),tahhP(2,jj(ij)),'x',axbt.latzh(kk(ik)),axbt.tazh(kk(ik)),'d');%xlim([37.95 38.05])

ii=find(jdy>40.7 & jdy<41);jj=find(jdxhP>40.7 & jdxhP<41);kk=find(axbt.jdzh>40.7 & axbt.jdzh<41);
figure;plot(lon(ii)-360,lat(ii),'o',lonhhP(jj),lathhP(jj),'x',axbt.lonzh(kk),axbt.latzh(kk),'d');xlabel('Longitude');ylabel('Latitude');title('CALWATER2 Feb 8, 2015 P3 Drops');
legend('RHB','Ta','SST');
ij=find(lonhhP(jj)<-124 & lonhhP(jj)>-129);ik=find(axbt.lonzh(kk)<-124 & axbt.lonzh(kk)>-129);
figure;plot(lat(ii),tsnk(ii),'.',lat(ii),ta(ii),'o',lathhP(jj(ij)),tahhP(2,jj(ij)),'x',axbt.latzh(kk(ik)),axbt.tazh(kk(ik)),'d');%xlim([37.95 38.05])

for i=1:jkP
    ii=find(jdy>jdxhP(i));
    dexP(i)=ii(1);
end;
jk=length(axbt.tazh)

for i=1:jk
    ii=find(jdy>axbt.jdzh(i));
    dex(i)=ii(1);
end;

ii=find(~isnan(axbt.tazh));
[p,s] = polyfit(axbt.latzh(ii),axbt.tazh(ii),1);tsx=polyval(p,axbt.latzh);
 aa=[ones(size(axbt.latzh(ii)))' axbt.latzh(ii)' axbt.lonzh(ii)'];cc=aa\axbt.tazh(ii)';     
 ab=[ones(size(axbt.latzh))' axbt.latzh' axbt.lonzh'];
 ac=[ones(size(lon(dex))) lat(dex) lon(dex)-360];
 sstx=ab*cc;
 sstrhb=ac*cc;

ii=find(~isnan(tahhP(2,:)));
[pP,s] = polyfit(lathhP(ii),tahhP(2,ii),1);tap=polyval(pP,lathhP);
 aa=[ones(size(lathhP(ii)))' lathhP(ii)' lonhhP(ii)'];cct=aa\tahhP(2,ii)';   
 ab=[ones(size(lathhP))' lathhP' lonhhP'];tap2=ab*cct;
 ac=[ones(size(lon(dexP))) lat(dexP) lon(dexP)-360];
 %tax=ab*cc;
 tarhb=ac*cct;

ii=find(~isnan(qahhP(2,:)));
[pP,s] = polyfit(lathhP(ii),qahhP(2,ii),1);qap=polyval(pP,lathhP);
 aa=[ones(size(lathhP(ii)))' lathhP(ii)' lonhhP(ii)'];ccq=aa\qahhP(2,ii)';     
 ab=[ones(size(lathhP))' lathhP' lonhhP'];qap2=ab*ccq;
 ac=[ones(size(lon(dexP))) lat(dexP) lon(dexP)-360];
 %tax=ab*cc;
 qarhb=ac*ccq;



distx=sqrt((lat(dex)'-axbt.latzh(1:jk)).^2+(lon(dex)'-360-axbt.lonzh(1:jk)).^2);
distp=sqrt((lat(dexP)'-lathhP(1:jkP)).^2+(lon(dexP)'-360-lonhhP(1:jkP)).^2);

figure;plot(lat(dex)'-axbt.latzh(1:jk),sstrhb'-tsnk(dex)','o',lat(dex)'-axbt.latzh(1:jk),axbt.tazh(1:jk)-tsnk(dex)','.',lat(dex)'-axbt.latzh(1:jk),.41*(lat(dex)'-axbt.latzh(1:jk)));xlabel('N-S Distance (deg)');ylabel('AXBT-Tsnk (C)');
figure;plot(lat(dexP)'-lathhP(1:jkP),ta(dexP)-tarhb,'o',lat(dexP)'-lathhP(1:jkP),tahhP(2,1:jkP)-ta(dexP)','.',lat(dexP)'-lathhP(1:jkP),.10*(lat(dexP)'-lathhP(1:jkP)));xlabel('N-S Distance (deg)');ylabel('SondeTa-Ta (C)');
figure;plot(lat(dexP)'-lathhP(1:jkP),qa(dexP)-qarhb,'o',lat(dexP)'-lathhP(1:jkP),qahhP(2,1:jkP)-qa(dexP)','.',lat(dexP)'-lathhP(1:jkP),.10*(lat(dexP)'-lathhP(1:jkP)));xlabel('N-S Distance (deg)');ylabel('Sondeqa-qa (C)');


figure;plot(distx,axbt.tazh(1:jk)-tsnk(dex)','x',distp,tarhb'-ta(dexP)','o',distx,sstrhb-tsnk(dex),'d');xlabel('Distance (deg)');ylabel('P3-RHB Temperature (C)');legend('SST','Ta');
figure;plot(distp,UhhP(2,1:jkP)-U(dexP)','o');xlabel('Distance (deg)');ylabel('P3-RHB Wind Speed(m/s)');
figure;plot(distp,qahhP(2,1:jkP)-qa(dexP)','o',distp,qarhb'-qa(dexP)','d');xlabel('Distance (deg)');ylabel('P3-RHB Spec Humidity(g/kg)');

kkx=find(distx<5);kkp=find(distp<5);
biastx=nanmedianx(axbt.tazh(kkx)-tsnk(dex(kkx))');biastxr=nanstdx(axbt.tazh(kkx)-tsnk(dex(kkx))');
biastx2=nanmedianx(sstrhb(kkx)-tsnk(dex(kkx)));biastxr2=nanstdx(sstrhb(kkx)-tsnk(dex(kkx)));

biasta=nanmedianx(tahhP(2,kkp)-ta(dexP(kkp))');biastar=nanstdx(tahhP(2,kkp)-ta(dexP(kkp))');
biasta2=nanmedianx(tarhb(kkp)-ta(dexP(kkp)));biastar2=nanstdx(tarhb(kkp)-ta(dexP(kkp)));

biasu=nanmedianx(UhhP(2,kkp)-U(dexP(kkp))');biasur=nanstdx(UhhP(2,kkp)-U(dexP(kkp))');
biasq=nanmedianx(qahhP(2,kkp)-qa(dexP(kkp))');biasqr=nanstdx(qahhP(2,kkp)-qa(dexP(kkp))');
[biastx biastx2 biasta biasta2 biasu biasq length(kkx) length(kkp)]
[biastxr biastxr2 biastar biastar2 biasur biasqr]

figure;plot(U(dexP)',UhhP(2,1:jkP),'o',U(dexP)',U(dexP)');xlabel('RHB Wind Speed(m/s)');ylabel('P3 Wind Speed(m/s)');ylim([0 15])
hold on;plot(U(dexP(kkp))',UhhP(2,kkp),'x');

tsax=axbt.tazh;
latax=axbt.latzh;
lonax=axbt.lonzh;
figure;%map_neaqs2004_2_cont;
hold on;
jj=find(~isnan(tsax));
%cdj=(taucx(jj)+tauib(jj))/2./gf(jj)./u10n(jj).^2./rhoair(jj)./cdn10(jj);
[xi,yi]=meshgrid(-132:.5:-122,35.5:.25:40.5);
%[xi,yi]=meshgrid(-71:.5:-66,41.5:.5:44.5);
ZI = griddata(lonax(jj),latax(jj),tsax(jj),xi,yi);
%v=[0 .2 .4 .6 .8 1.0  1.2];
v=13:.5:17.5;%[ -.2 -.1 -.05 0 .05 .1 .25  .5 1 2 ];
[c,h] = contourf(xi,yi,ZI,v); clabel(c,h), colorbar;
title('CALWATER2 1-27 to 2-09 2015 NOAA P-3 AXBT SST (C)');
xlabel('Longitude (deg)');ylabel('Latitude (deg)');
hold off;





