% epic01_Rscldar_5  computes clear sky solar for ETL data
%jdy is julian day (decimal) at the start of average interval
%lat and lon are decimal latitude and longitude
%iv is column water vapor in cm
%   if no data is available for iv, put in NaN and select a ratio of iv to
%   near-surface water vapor mixing ratio, qrat.  Typical tropical value is
%   4.5
%calls SolarRadiance
%
close all
clear

way_raw_scs = 'y:\knorr\scs\';

% This is the cvs daily file and not the SCS from our system
md=input('Input start mon and day to evaluate (0323 on).... ','s');

% sounding PWV in ICEALOT__PWV_sonde.mat c:\data\ICEALOT\data\sonde\
ss = 'c:\data\ICEALOT\data\sonde\ICEALOT_PWV_sonde.mat';
load(ss);
iv = tpw(length(tpw),2);   % est of total col water vapor for Seattle on JD 354

%find all the scs txt files
fn2=['*.txt'];
s2=[way_raw_scs fn2];
%s1=['c:\data\icealot\data\Mailbox\' fname];
dr2=dir(s2);
[n2 m2]=size(dr2);
mx2=n2;
clear scs
%ICEALOT_scs_MMDD.txt
flnmscs = dr2(n2).name
f=[way_raw_scs flnmscs];
ss2 = [f];
scs = load(ss2);

[r,c] = size(scs);
jdys = scs(:,c); %load julian day from scs into array

jjx = 1;
qrat=4.;
latd = scs(1,6); % ICEALOT
lond = scs(1,7);
k1 = .05;                  %aerosol optical depth, band 1
k2 = .05;                  %aerosol optical depth, band 2
oz = 0.2;                  %column ozone
%%%%%% DEW 12262007

Rh = scs(:,16);
T = scs(:,3);
p = mean(scs(:,12));
scssw = scs(:,13);
q=Rh./100.*qsea(T);
%%%%
%function y=qsea(x)
%es=6.112.*exp(17.502.*x./(x+241.0))*.98*(1.0007+3.46e-6*1010);
%y=es*622./(1010.0-.378*es);
%%%%
%QA_mean=[QA_mean qvais];
%qa=QA_mean(1:jjx);


watvap(1:length(jdys),1)=iv;
%ii=find(isnan(iv));
%watvap(ii)=qa(ii)/qrat;%total column water vapor

jdx=floor(jdys+.5/60/24); %julian day in bin centers for 10-min aves


tutc=(jdys+5/60/24-jdx)*24;%0:1:23;%hour of the day
[n m]=size(tutc);
clear sw sz saz dirs sky;
[sw,sz,saz,dirs,sky] = SolarRadiancex(latd,lond,jdx,tutc,watvap,p,k1,k2,oz);
Rscl=sw;
%plot(jdy,Rscl,jdy,rs,'o');xlabel('Julian Day (1999)');ylabel('Downward Solar Flux (W/m^2)');
%ss=stdt;axis([ss endt 0 1200]);

%sw = downward shortwave (W/m^2)
%sz=zenith angle
%saz = azimuth angle
%dirs is direct solar
%sky is diffuse solar
figure
plot(jdys,Rscl)
hold
plot(jdys,sw,'g')
title('ICEALOT 2008')
ylabel('SW (W/m2')
xlabel('Year Day (UTC)')
legend('Clear','SCS')
mmin = min(samos(:,16));
mmax = max(samos(:,16));
ylim([min(mmin)-25 max((mmax)+ 25)]);