%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computes the solar clear sky values using microwave radiometer data and Iqbal model
% Also calculates the zenith declination and LW clear sky.
%LB - 04/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;fclose all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%first load most recent data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MET QC data
% metdataleg2=load('S:\cruises\DYNAMO_2011\Revelle\Scientific_analysis\PSD_UConn_OSU_MeanMetFiles\Revelle10minuteLeg2.r1');
metdataleg3=load('S:\cruises\DYNAMO_2011\Revelle\Scientific_analysis\PSD_UConn_OSU_MeanMetFiles\Revelle10minuteLeg3.r1');
% metdata=[metdataleg2; metdataleg3]; 
% clear metdataleg2 metdataleg3
metdata=[metdataleg3]; 
clear metdataleg3
metyday=metdata(:,1);    %Decimal yearday (UTC)
Lat=metdata(:,2);        %Latitude (deg)
Lon=metdata(:,3);        %Longitude (deg)
Pair10=metdata(:,13);    %Pressure (mb) adjusted to 10 m
T10=metdata(:,15);       %Temperature (C) adjusted to 10 m
Q10=metdata(:,18);       %Specific humidity (g/Kg) adjusted to 10 m
Solardn=-metdata(:,26);   %Measured downwelling solar (W/m2)
IRdn=-metdata(:,28);      %Measured downwelling IR (W/m2)

% Mailbox data
% Mailboxleg2=load('S:\cruises\DYNAMO_2011\Revelle\radiometer\Mailbox\Processed\leg2\DYNAMO_2011_micro_m_leg2.txt');
Mailboxleg3=load('S:\cruises\DYNAMO_2011\Revelle\radiometer\Mailbox\Processed\leg3\DYNAMO_2011_micro_m_leg3.txt');
% mailboxdata=[Mailboxleg2; Mailboxleg3]; 
% clear Mailboxleg2 Mailboxleg3
mailboxdata=[Mailboxleg3];
clear Mailboxleg3 
mbyday=mailboxdata(:,1);  %decimal day-of-year
mbIWV=mailboxdata(:,2);   %integrated atmospheric water vapor (column precipitable water), cm
% mb_IWL=mailboxdata(:,3);   %integrated cloud liquid content, cm
mbIWVstd=mailboxdata(:,5);  %standard deviation of IWV, cm
% mb_IWLstd=mailboxdata(:,6);  %standard deviation of IWL, cm

% Ceilo data
% ceiloleg2=load('S:\cruises\DYNAMO_2011\Revelle\ceilometer\Processed\leg2\DYNAMO_2011_ceilo_10min_a_leg2.txt');
ceiloleg3=load('S:\cruises\DYNAMO_2011\Revelle\ceilometer\Processed\leg3\DYNAMO_2011_ceilo_10min_a_leg3.txt');
% ceilodata=[ceiloleg2; ceiloleg3]; 
% clear ceiloleg2 ceiloleg3
ceilodata=[ceiloleg3]; 
clear ceiloleg3
ceiloyday=ceilodata(:,1);  %decimal day-of-year
tf= ismember(ceiloyday, metyday);  %finds same point as met time series
ceiloyday=ceiloyday(tf);
clrfrac=ceilodata(tf,8);   %clear fraction
cldfrac=ceilodata(tf,9);   %cloud fraction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%check plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);clf;
subplot(2,1,1);plot(metyday,Solardn,'.');ylim([-5 1200]);grid;title('Measured downwelling solar (W/m2)');
subplot(2,1,2);plot(metyday,IRdn,'.');grid;title('Measured downwelling IR (W/m2)')
xlabel('DOY');

figure(2);plot(mbyday,mbIWV);
xlabel('DOY');title('Integrated atmospheric Water Vapor (cm)');
grid

ii=find(mbIWVstd>.15 | mbIWV==0 | mbIWV<3);
figure(3);plot(mbIWVstd,mbIWV,'.',mbIWVstd(ii),mbIWV(ii),'r.');grid
xlabel('IWV std');ylabel('IWV');

Q10i=interp1(metyday,Q10,mbyday);
% mbIWV(ii)=Q10i(ii)/qrat;
for oo=1:length(ii)
    mbIWV(ii(oo))=mbIWV(ii(oo)-1);%do a couple of times to get rid of points
end    
figure;plot(mbyday,mailboxdata(:,2),'b',mbyday,mbIWV,'r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%do calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ll=find(clrfrac==1);
figure;plot(iv(ll),IRdn(ll)-rlclr(ll),'.');grid;ylim([-30 30])
figure;plot(rlclr(ll),IRdn(ll),'.',[390 440],[390 440],'k--');grid;axis([390 440 390 440])

% [3.9 17]/(5.67e-8*(nanmean(T10)+273.15).^4)
% [3.9 21]/(5.67e-8*(nanmean(T10)+273.15).^4)
[4.1 22]/(5.67e-8*(nanmean(T10)+273.15).^4)

%LW clear sky
rlclr=(.52+.13/60*abs(Lat)+(.082-.03/60.*abs(Lat)).*sqrt(Q10)).*(5.67e-8*(T10+273.15).^4);
% rlclr=(.52+.13/60*abs(Lat)+(.084-.03/60.*abs(Lat)).*sqrt(Q10)).*(5.67e-8*(T10+273.15).^4);
% rlclr2=(.52+.13/60*abs(Lat)+(.082-.03/60.*abs(Lat)).*sqrt(Q10)-0.0288+0.0063*iv).*(5.67e-8*(T10+273.15).^4);
% rlclr2=(.52+.13/60*abs(Lat)+(.082-.03/60.*abs(Lat)).*sqrt(Q10)+0.0084*iv-0.0365).*(5.67e-8*(T10+273.15).^4);
% rlclr2=(.52+.13/60*abs(Lat)+(.084-.03/60.*abs(Lat)).*sqrt(Q10)+0.0084*iv-0.0451).*(5.67e-8*(T10+273.15).^4);
% rlclr2=(0.58+0.07.*sqrt(Q10)+0.0088*iv-0.0473).*(5.67e-8*(T10+273.15).^4);

%this ones seems better (no adjustement on general A and B coeffs but use
%of IV
rlclr2=(.52+.13/60*abs(Lat)+(.082-.03/60.*abs(Lat)).*sqrt(Q10)+0.0084*iv-0.0451).*(5.67e-8*(T10+273.15).^4);

figure(4);plot(metyday,IRdn,'.--',metyday,rlclr,'r',metyday,rlclr2,'m');grid;legend('Measured downwelling IR (W/m2)','Clear Sky','Clear sky with IV')


%SW clear sky
qrat=3.3; 
k1=.035;                  %aerosol optical depth, band 1
k2=.035;                  %aerosol optical depth, band 2
% k1=.05;                  %aerosol optical depth, band 1
% k2=.05;                  %aerosol optical depth, band 2

oz=.3;                  %column ozone
iv=ones(length(Q10),1)*NaN;  %iv is column water vapor in cm
iv(isnan(iv))=Q10(isnan(iv))/qrat;%total column water vapor
% iv(isnan(iv))=Q10(isnan(iv))*5/18;%total column water vapor
%if no data is available for iv, put in NaN and select a ratio of iv to
%near-surface water vapor mixing ratio, qrat.  
% jdx=floor(metyday+.5/60/24); %julian day in bin centers for 10-min aves
jdx=floor(metyday+.5/60/24); %julian day in bin centers for 10-min aves
tutc=(metyday+0/60/24-jdx)*24;%0:1:23;%hour of the day

clear rsclr sz saz dirs sky;
[rsclr1,sz,saz,dirs,sky] = SolarRadiancex(Lat,Lon,jdx,tutc,iv,Pair10,k1,k2,oz);
%rsclr= downward shortwave (W/m^2)
%sz=zenith angle
%saz = azimuth angle
%dirs is direct solar
%sky is diffuse solar
% figure(5);plot(metyday,Solardn,'.',metyday,rsclr1,'r');ylim([-5 1200]);
% grid;legend('Measured downwelling solar (W/m2)','Clear Sky with q/5')

iv=interp1(mbyday,mbIWV,metyday);  %iv is column water vapor in cm
iv(isnan(iv))=Q10(isnan(iv))/qrat;%total column water vapor
iv(iv<3)=Q10(iv<3)/qrat;%total column water vapor
clear sw sz saz dirs sky;
[rsclr2,sz,saz,dirs,sky] = SolarRadiancex(Lat,Lon,jdx,tutc,iv,Pair10,k1,k2,oz);
%rsclr= downward shortwave (W/m^2)
%sz=zenith angle
%saz = azimuth angle
%dirs is direct solar
%sky is diffuse solar
figure(5);plot(metyday,Solardn,'.--',metyday,rsclr1,'r',metyday,rsclr2,'g');ylim([-5 1200]);
grid;legend('Measured downwelling solar (W/m2)','Clear Sky with q/5','Clear Sky with mailbox IWV')

% figure(6);plot(metyday,saz,'.--',metyday,sz,'r');
% grid;legend('azimuth','zenith')
figure(6);plot(metyday,Solardn,'.--',metyday,sz,'k')
grid;legend('Measured downwelling solar (W/m2)','Zenith (deg)')

figure;plot(rsclr1(ll),Solardn(ll),'.',rsclr2(ll),Solardn(ll),'g.');grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%saves data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ClrSkyMat=[metyday rsclr2 sz rlclr2];
save('C:\Data\CRUISES\DYNAMO_2011\data\ClearSky_leg3.mat', 'ClrSkyMat')

sdg=['save ', 'C:\Data\CRUISES\DYNAMO_2011\data\ClearSky_leg3.prn ClrSkyMat -ascii'];
eval(sdg)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For fun reads the AOD from aeronet (NASA)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flist=fopen('C:\Data\CRUISES\DYNAMO_2011\Analysis\Roger_Revelle_11_0\AOD\Roger_Revelle_11_0_series.lev20','r');
temp=textscan(flist,'%2f:%2f:%4f %2f:%2f:%2f %*2c %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %*2c %*[^\n]','delimiter', ', ', 'headerlines', 5,  'emptyvalue', NaN,'treatAsEmpty', 'Sonic');
RR_aod=cell2mat(temp);
RR_doy=datenum(RR_aod(:,3),RR_aod(:,2),RR_aod(:,1),RR_aod(:,4),RR_aod(:,5),RR_aod(:,6))-datenum(RR_aod(:,3)-1,12,31);
RR_aodlat=RR_aod(:,8);
RR_aodlon=RR_aod(:,9);
RR_aod380nm=RR_aod(:,11);
RR_aod500nm=RR_aod(:,13);

figure;plot(RR_doy,RR_aod380nm,'b+--',RR_doy,RR_aod500nm,'r+--')


% k1=ones(1,length(Lat))*.035;                  %aerosol optical depth, band 1
% k2=ones(1,length(Lat))*.035;                  %aerosol optical depth, band 2
k1=interp1(RR_doy,RR_aod380nm/10,metyday);
k2=interp1(RR_doy,RR_aod380nm/10,metyday);

% ii=find(RR_doy>=metyday(1) & RR_doy<=metyday(end))
% RR_doy=RR_doy(ii);RR_aod380nm=RR_aod380nm(ii);RR_aod500nm=RR_aod500nm(ii);
% for ll=1:length(RR_doy)
%     oo=find(metyday-RR_doy(ll)==min(metyday-RR_doy(ll)))
% end
clear sw sz saz dirs sky;
[rsclr3,sz,saz,dirs,sky] = SolarRadiancex(Lat,Lon,jdx,tutc,iv,Pair10,k1,k2,oz);
figure(5);plot(metyday,Solardn,'.--',metyday,rsclr1,'r',metyday,rsclr2,'g',metyday,rsclr3,'m');ylim([-5 1200]);
grid;legend('Measured downwelling solar (W/m2)','Clear Sky with q/5','Clear Sky with mailbox IWV')
