% epic01_micro90red.m

% reads in microwve .sci files files and writes out
% one datafile with all of the pertinent info in it.

% program written by C. Fairall.  Modified by D. E. Lane
% on 5 March 2001.  Added 10-min statistics.

disp('Reading mirowave files')
clear all;
fclose('all');
s='f:\epic2\203190\sci\';
dr=dir(s);			%lists all files in directory
[n m]=size(dr);
nx=n;%				%# of files in data directory (includes . and ..)     


%******************  setup read data loop  **********
ss='f:\epic2\203190\sci\';

%nx=24*5+2;
jamx=1;
for ibg = 3:nx 			%major read loop
   e=[ss dr(ibg).name]; %gets names of files in directory
   stux=load(e);   		%loads all data from files into stux
   
%******   decode date  ************
   if isempty(stux)		% if no data, end the program
      break;
   end;
    [nl ml]=size(stux);  		%nl=file length - m1 is number of cols
    ff=dr(ibg).name;		 		%names of files
    month=str2num(ff(3:4));	%get month
    day=str2num(ff(5:6));		%get day
    year=str2num(ff(1:2));		%get year
    hr=str2num(ff(7:8));
    min=floor((stux(:,1)-str2num(ff(1:8))*1e4)/100);
    secs=stux(:,1)-str2num(ff(1:8))*1e4-min*100;
    jd0=datenum(2000+year,month,day)-datenum(2000+year-1,12,31); % calculates julian date
    jd=jd0+((hr+min/60+secs/3600)/24); % calculate decimal jdate
   
    happ(jamx:jamx+nl-1,1:18)=[jd stux(:,1:17)];  % make data array
    jamx=jamx+nl  
end; %  data line loop

% plot 15-sec data

figure(1);plot(happ(:,1),happ(:,5),'o'); xlabel('Julian Day'); ylabel('Integrated WV (cm)');
figure(2);plot(happ(:,1),happ(:,6)*1e3,'.');xlabel('Julian Day');ylabel('LWP (gm/m^2)');
figure(3);plot(happ(:,1),happ(:,8:9),'.',happ(:,1),happ(:,14));xlabel('Julian Day');ylabel('LWP (gm/m^2)');

% save 15-sec data

micro_sec=happ;					% this is the raw data with date
save d:\cwf\data\epic2001\micro\micro_sec_90.txt micro_sec -ascii

%******************  setup read data loop  ********** 

ny=jamx-1;
lam=0:24;  			  % hours
jdx=lam/24;			  % fraction of day
jdd=happ(:,1);		  % date vector
jd0=floor(jdd(1));  % first Julian Day of data series
jxm=floor(jdd(ny)); % last Julian Day of data series
jx=1;					  % counts number of hours of data
kx=1;					  % counts number of 10-min blocks of data
min=0:6;				  % 10-min parts of an hour
jdmx=min/6/24;			  % 10-min fraction of hour

while jd0<=jxm+1  % the first day is l.t.e. last day
   ii=find(jdd>jd0 & jdd<(jd0+1));  % find all data
   if isempty(ii)  % if no data
      jd0=jd0+1;   % increment the start date
   else
   	jj=1;			 % if not empty start first hour calc.
      while jj<25  % not the end of the day
         pp=jdd(ii)-floor(jdd(ii)); % fractional part of day
      	ij=find(pp>jdx(jj) & pp<jdx(jj+1)); % select hour
         if isempty(ij)
            alldata(jx)=0;
            vap(jx)=NaN;
            liq(jx)=NaN;
            tb(jx,1:3)=NaN;
      	else
            kk=1;			 % if hourly data exists start 10-min calc.
            while kk<7
               kj=find(pp(ij)>(pp(ij(1))+jdmx(kk)) & pp(ij)<(pp(ij(1))+jdmx(kk+1)));  % select 10-min section
               if isempty(kj)
                  mindata(kx)=0;
                  vapm(kx)=NaN;
                  liqm(kx)=NaN;
                  tbm(kx,1:3)=NaN;
               else   
                  [mindata(kx) c]=size(pp(kj));  % how much data in 10-min section
                  vapm(kx)=median(happ(ii(ij(kj)),5));	% total possible cloud
                  liqm(kx)=median(happ(ii(ij(kj)),6));
                  tbm(kx,1:2)=median(happ(ii(ij(kj)),8:9));
                  tbm(kx,3)=median(happ(ii(ij(kj)),14));  
                  
               end; %end if isempty kj
               jdmin(kx)=jd0+jdx(jj)+jdmx(kk); % integrate timestep
               kk=kk+1;
               kx=kx+1;
            end  %while kk      
            [alldata(jx) b]=size(pp(ij));		% all data in selected hour 
             vap(jx)=median(happ(ii(ij),5));	% total possible cloud
             liq(jx)=median(happ(ii(ij),6));
             tb(jx,1:2)=median(happ(ii(ij),8:9));
             tb(jx,3)=median(happ(ii(ij),14));  

         end;%end if isempty ij
         jdh(jx)=jd0+jdx(jj);
     	   jj=jj+1;
      	jx=jx+1;
      end;%end while jj
      jd0=jd0+1
   end;%end if isempty ii
end;%end while jd0

% plot rest of data

figure(3);plot(jdh,vap,'o',jdmin,vapm,'.');
xlabel('Julian Day (2001)');ylabel('Hourly (o) and 10-min (.) vap (cm)');

figure(4);plot(jdh,liq,'o',jdmin,liqm,'.');
xlabel('Julian Day (2001)');ylabel('Hourly (o) and 10-min (.) liq (mm)');

figure(4);plot(jdh,liq,'o',jdmin,liqm,'.');
xlabel('Julian Day (2001)');ylabel('Hourly (o) and 10-min (.) liq (mm)');

figure(5);plot(jdh,tb,'o',jdmin,tbm,'.');
xlabel('Julian Day (2000)');ylabel('Hourly (o) and 10-min (.) Tb (K)');


% save the data in hourly and 10-min files

micro_hourly=[jdh' alldata' vap' liq' tb];       
save  d:\cwf\data\epic2001\micro\micro_hr_90.txt micro_hourly -ascii

micro_10min=[jdmin' mindata' vapm' liqm' tbm];
save  d:\cwf\data\epic2001\micro\micro_min_90.txt micro_10min -ascii

% clean up memory

%clear
%pack



	
	
	
	
