% rain calculation from OSI
clear

read_parameters_VOCALS

jds =input('Input year day to evaluate (298....): ','s');

fnmet=['VOCALS_met_' jds '.txt'];
ss3 = way_proc_data_flux;
s3 = [ss3 fnmet];

met = load(s3);
bd = fix(met(1,63));
ed = bd + 1;

met = load(s3);
metjad = met(:,63);
yd = num2str(fix(metjad(1,1)));

%add 3 to number for position in array
%1 176 L
%2 Year_RTM  L
%3 Day_RTM  L
%4 Hour_Minute_RTM  L
%5 caseRpir1_AVG  H
%6 domeRpir1_AVG  H
%7 caseRpir2_AVG  H
%8 domeRpir2_AVG  H
%9 seatemp_r_AVG  H
%10 vaisT_v_AVG  H
%11 vaisRH_v_AVG  H
%12 pir1_v_AVG  H
%13 pir2_v_AVG  H
%14 psp1_v_AVG  H
%15 psp2_v_AVG  H
%16 ORG1C_v_AVG  H
%17 ORG1S_v_AVG  H
%18 caseCpir1_AVG  H
%19 domeCpir1_AVG  H
%20 caseCpir2_AVG  H
%21 domeCpir2_AVG  H
%22 seatempC_AVG  H
%23 vaisT_c_AVG  H
%24 vaisRH_pc_AVG  H
%25 pir1_wm2_AVG  H
%26 pir2_wm2_AVG  H
%27 psp1_wm2_AVG  H
%28 psp2_wm2_AVG  H
%29 ORG1_mmhr_AVG  H
%30 bv_v_AVG  H
%31 airflow_AVG  H
%32 press_V_AVG  H
%33 hpa_AVG  H
%%% Calculate LW from PIRs

metrain = met(:,32);
metrainv = met(:,20);

% Calculate rain rate from V  25(V**1.87)-.15
% .085 VDC = .099 mm/hr  1 = 24.85
for s = 1:length(metrainv)
mmhr(s,1) = 25 * (metrainv(s,1)^1.87)-.15;
end
t = isnan(mmhr) == 1;
for tt = 1:length(t)
    if t(tt) == 1
        tt
        t(tt)
        mmhr(tt) = -.09;
    end
end
for tt = 1:length(t)
    if mmhr(tt) < 0
        mmhr(tt) = mmhr(tt) +.1;
    end
end
cmmhr = cumsum(mmhr)/30;
% Rain **********
figure(1)
clf
plot(metjad,metrain,'r')
hold
plot(metjad,mmhr,'k')
legend('metrain','mmhr','Location','Best')
tit = {'Rain comparison VOCALS 2008  ' [yd]};
title(tit)
xlabel('Year Day (UTC)')
ylabel('Rain')

%Convert PSD optical rain gauge from mm/hr to mm total
% PSD ORG has an offset of ~-.028 mm/hr
t = isnan(metrain) == 1;
for tt = 1:length(t)
    if t(tt) == 1
        tt
        t(tt)
        metrain(tt) = -.027;
    end
end
for tt = 1:length(t)
    if metrain(tt) < 0
        metrain(tt) = metrain(tt) +.028;
    end
end
psdsum = (cumsum((metrain)))/30;

figure(2)
clf
plot(metjad,psdsum,'r')
hold
plot(metjad,cmmhr,'b')
legend('mmhr','V','Location','Best')
tit = {'PSD VOCALS 2008  ' yd};
title(tit)
xlabel('Year Day (UTC)')
ylabel('mm')
print('-djpeg90 ',[way_raw_images_flux '\PSD_cum_rain_' jds '.jpg']);

