%BAOPLT300_month_loop
% Plot data from GMD 300m BAO data loggers
% Reads in 10 min file which has a different format than the 30 sec data
% Dave has split these data into daily files and corrected the PRT and Wind Dir
%You are correct on all counts.  yes, m/s on the wind speed.
%The PRT temperature is (accurate to within 0.1 deg C) is:  a0 +
%a1*(D3_PRTRES_Avg) + a2*(D3_PRTRES_Avg)^2 + a3*(D3_PRTRES_Avg)^3,
% a0 =    -246.058;
% a1 =      2.36626;
% a2 =  0.000923807;
% a3 =  2.43906e-007;

%The sonic anemometer south to north pole point to 322 degrees.
%The the true wind direction should then be    (SE2_SONICWD_Avg) - 38 degrees
%
% 300m
%"TOA5","CR1000","CR1000","2846","CR1000.Std.11","CPU:060516b_Met.CR1","10044","Public"
%"TIMESTAMP","RECORD","BT_BATVOLT","SE1_SONICWS","SE2_SONICWD","SE3_SONICTEMP","D3_PRTTEMPAPPROX","SE4_HMPTEMP","SE9_HMPRH","SE10_ACPOWER","P1_AIRFLOW","SW12_12VFLAG","D3_PRTRES"
%"TS","RN","Volts","mV","mV","mV","Deg C","Deg C","%","mV","Counts","",""
%"2007-06-14 14:18:58.79",4479243,13.32056,3.517115,229.1796,13.54004,14.05173,14.93604,42.78085,110.567,0,1,105.7731
%"2007-06-14 14:28:56.39",4479841,13.2959,3.516633,231.2331,13.42691,14.39247,15.26545,40.26095,110.634,0,1,105.9062

function [retval] = BAOPLT300_month_loop_func( by, bm, bd, ey, em, ed ) 
%clear

pltflg = 1; %flag to plot the data: 1=yes, 0=no
pngflg = 1; %flag to make a png of the plot: 1=yes, 0=no

if ispc
    inpath = 'z:\BAO\Tower\Processed\monthly\';
    outpath  = 'z:\BAO\Tower\Processed_Images\monthly\';
else
    inpath = '/archive/BAO/Tower/Processed/monthly/';
    outpath  = '/archive/BAO/Tower/Processed_Images/monthly/';
end

%cd(inpath);

if( nargin == 0 )
    bys = input('Input year to process (2007) ','s');
    bms = input('Input the month to process (1-12) ','s');
    bds = input('Input the beginning day to process (1-31) ','s');
    eys = input('Input year to process (2007) ','s');
    ems = input('Input the month to process (1-12) ','s');
    eds = input('Input end day to process (1-31) ','s');
    by = str2double( bys );
    bm = str2double( bms );
    bd = str2double( bds );
    ey = str2double( eys );
    em = str2double( ems );
    ed = str2double( eds );
end

bdn = datenum([by bm bd 0 0 0]);
edn = datenum([ey em ed 0 0 0]);
bdv=datevec(bdn);
edv=datevec(edn);
yr = bdv(1);
mo = bdv(2);

nloops = edv(2) + (edv(1) - bdv(1)) * 12 - bdv(2) + 1;

nmax = 90000;

for i = 1:nloops
    name300 = sprintf('%sBAO_300_%4d%02d.dat', inpath, yr, mo );
    fprintf(1,'*** Processing 300m monthly ***\n');
    
    fid300 = fopen( name300 );
    if( fid300 > 0 )
        fclose( fid300 );
    else
        fprintf(1,'No  300m: %s\n',name300);
        continue;
    end
    
    [data300 dt n] = load300( name300, nmax);
        
    jdb = fix(dt(1,7));
    jde = fix(dt(n,7)) + 1;

    if pltflg % do plots if pltflg=1
        figure1 = figure('Position',[10 50 650 650]);
        subplot(2,1,1),plot(dt(1:n,7),data300(1:n,5),'k')
        hold
        subplot(2,1,1),plot(dt(1:n,7),data300(1:n,8),'g')
        subplot(2,1,1),plot(dt(1:n,7),data300(1:n,13),'b')
        gis = find( data300(:,5) > -50. );
        gih = find( data300(:,8) > -50. );
        gip = find( data300(:,13) > -50. );
        blims(1) = min(data300(gip,13))-5;
        elims(1) = max(data300(gip,13))+5;
        blims(2) = min(data300(gih,8))-5;
        elims(2) = max(data300(gih,8))+5;
        %blims(3) = min(data300(gis,5))-5;
        %elims(3) = max(data300(gis,5))+5;
        %xlim([jdb,jde])
        stitle=sprintf('BAO %4d-%02d 300m', yr, mo);
        title(stitle)
        ylabel('T (C)')
        xlabel('Day (UTC)')
        datetick('x','mm/dd','keeplimits','keepticks');     
        legend('Sonic T','HMP','PRT corrected','Location','Best','Orientation','horizontal')
        legend('boxoff');
        axis tight
        grid on
        subplot(2,1,2),plot(dt(1:n,7),data300(1:n,9))
        %xlim([jdb,jde])
        xlabel('Day (UTC)')
        datetick('x','mm/dd','keeplimits','keepticks');     
        ylabel('Rh (%)')
        axis tight
        grid on
        ylim([0,109])
        annotation(figure1,'textbox',[0.007154 0.01077 0.4498 0.02462],...
                  'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},...
                  'FontSize',8,'FitBoxToText','off','LineStyle','none');
        if pngflg
            print_buffer = sprintf('%sBAO_300_%4d%02d_TRH.png',...
                                   outpath, yr, mo );
            set(gcf,'PaperPositionMode','auto');
            print('-dpng', print_buffer);
            %f=getframe(gcf);
            %imwrite(f.cdata,print_buffer,'png');
        end

        figure2 = figure('Position',[10 50 650 650]);
        subplot(2,1,1),plot(dt(1:n,7),data300(1:n,1))
        ylim([0,max(data300(1:n,1))+5])
        xlim([jdb,jde])
        ylabel('Spd (m/s)')
        title(stitle)
        xlabel('Day (UTC)')
        datetick('x','mm/dd','keeplimits','keepticks');     
        axis tight
        grid on
        subplot(2,1,2),plot(dt(1:n,7),data300(1:n,14))
        %xlim([jdb,jde])
        ylabel('Dir (true)')
        xlabel('Day (UTC)')
        datetick('x','mm/dd','keeplimits','keepticks');     
        axis tight
        grid on
        ylim([0,360])
        annotation(figure2,'textbox',[0.007154 0.01077 0.4498 0.02462],...
                  'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},...
                  'FontSize',8,'FitBoxToText','off','LineStyle','none');
        if pngflg
            print_buffer = sprintf('%sBAO_300_%4d%02d_WS.png',...
                                   outpath, yr, mo );
            set(gcf,'PaperPositionMode','auto');
            print('-dpng', print_buffer);
            %f=getframe(gcf);
            %imwrite(f.cdata,print_buffer,'png');
        end

        figure3 = figure('Position',[10 50 650 650]);
        plot(dt(1:n,7),data300(1:n,5),'k')
        hold
        plot(dt(1:n,7),data300(1:n,7),'r')
        plot(dt(1:n,7),data300(1:n,8),'b')
        plot(dt(1:n,7),data300(1:n,13),'g')
        xlabel('Year Day (UTC)')
        ylabel('T (C)')
        stitle=sprintf('BAO %4d-%02d 300m', yr, mo);
        title(stitle)
        xlim([jdb,jde])
        gis = find( data300(:,5) > -50. );
        gie = find( data300(:,7) > -50. );
        gih = find( data300(:,8) > -50. );
        gip = find( data300(:,13) > -50. );
        blims(1) = min(data300(gip,13))-10;
        elims(1) = max(data300(gip,13))+10;
        blims(2) = min(data300(gie,7))-10;
        elims(2) = max(data300(gie,7))+10;
        blims(3) = min(data300(gih,8))-10;
        elims(3) = max(data300(gih,8))+10;
        %blims(4) = min(data300(gis,5))-10;
        %elims(4) = max(data300(gis,5))+10;
        xlabel('Day (UTC)')
        datetick('x','mm/dd','keeplimits','keepticks');     
        legend('Sonic T','PRT est','HMP T','PRT corrected','Location','North','Orientation','horizontal')
        legend('boxoff');
        axis tight
        grid on
        ylim([min(blims) max(elims)])
        annotation(figure3,'textbox',[0.007154 0.01077 0.4498 0.02462],...
                  'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},...
                  'FontSize',8,'FitBoxToText','off','LineStyle','none');
        if pngflg
            print_buffer = sprintf('%sBAO_300_%4d%02d_T.png',...
                                   outpath, yr, mo );
            set(gcf,'PaperPositionMode','auto');
            print('-dpng', print_buffer);
            %f=getframe(gcf);
            %imwrite(f.cdata,print_buffer,'png');
        end

    end % if doing plots
    
    clear data300 dt
    close all
    
        mo = mo + 1;
    if mo > 12
        mo = mo - 12;
        yr = yr + 1;
    end

end % loop through all files

retval = 0;