function[seasonalmean,seasonalmedian,seasonalerror,seasonalct] = seasonalMeanFigureSubroutine(season,cloudlevel,qcChoice,sel_year,lat,months0,clouds,diurnalCycles,qc_inds,vert_rezs,dates,cycleHours)
    if lat>0
        spr1=3;
        spr2=4;
        spr3=5;
        sum1=6;
        sum2=7;
        sum3=8;
        fal1=9;
        fal2=10;
        fal3=11;
        win1=12;
        win2=1;
        win3=2;
    else
        spr1=9;
        spr2=10;
        spr3=11;
        sum1=12;
        sum2=1;
        sum3=2;
        fal1=3;
        fal2=4;
        fal3=5;
        win1=6;
        win2=7;
        win3=8;
    end
    years=year(dates);
    if isa(sel_year,'double')
        sel_year=int2str(sel_year);
    end
    
    if strcmp(sel_year,'all')==0
        sel_year=str2double(sel_year);
        yearIdx=find(years==sel_year);
        specClouds=clouds(yearIdx,:);
        specCycles=diurnalCycles(yearIdx,:);
        specRezs=vert_rezs(yearIdx,:);
        specQCs=qc_inds(yearIdx,:);
        months0=months0(yearIdx);
    elseif strcmp(sel_year,'all')==1
        specClouds=clouds;
        specCycles=diurnalCycles;
        specRezs=vert_rezs;
        specQCs=qc_inds;
    end
        
    
    if isa(season,'double')
        season=int2str(season);
    end    
        
    if strcmp(season,'summer')
        seasonIdx=find(months0==sum1|months0==sum2|months0==sum3);
    elseif strcmp(season,'fall')
        seasonIdx=find(months0==fal1|months0==fal2|months0==fal3);
    elseif strcmp(season,'winter')
        seasonIdx=find(months0==win1|months0==win2|months0==win3);
    elseif strcmp(season,'spring')
        seasonIdx=find(months0==spr1|months0==spr2|months0==spr3);
    elseif strcmp(season,'all')
        seasonIdx=find(months0<13);
    elseif str2double(season)==1
        seasonIdx=find(months0==1);
    elseif str2double(season)==2
        seasonIdx=find(months0==2);
    elseif str2double(season)==3
        seasonIdx=find(months0==3);
    elseif str2double(season)==4
        seasonIdx=find(months0==4);
    elseif str2double(season)==5
        seasonIdx=find(months0==5);
    elseif str2double(season)==6
        seasonIdx=find(months0==6);
    elseif str2double(season)==7
        seasonIdx=find(months0==7);
    elseif str2double(season)==8
        seasonIdx=find(months0==8);
    elseif str2double(season)==9
        seasonIdx=find(months0==9);
    elseif str2double(season)==10
        seasonIdx=find(months0==10); 
    elseif str2double(season)==11
        seasonIdx=find(months0==11);
    elseif str2double(season)==12
        seasonIdx=find(months0==12);
    else
        disp('Please input season as summer, fall, winter, spring, or the number of a month as a string');
    end
    
    specClouds=specClouds(seasonIdx,:);
    specCycles=specCycles(seasonIdx,:);
    specRezs=specRezs(seasonIdx,:);
    specQCs=specQCs(seasonIdx,:);
    
    
    if strcmp(cloudlevel,'clear')
        cloudIdx=(specClouds<.15);
    elseif strcmp(cloudlevel,'partially cloudy')
        cloudIdx=(specClouds>=.15 & specClouds<=.85);
    elseif strcmp(cloudlevel,'cloudy')
        cloudIdx=(specClouds>.85);
    elseif strcmp(cloudlevel,'clr+partial')
        cloudIdx=(specClouds>=0 & specClouds<=.85);
    elseif strcmp(cloudlevel,'all')
        cloudIdx=(specClouds>=0 | isnan(specClouds));
    else
        disp('Please input cloudlevel as clear, partially cloudy, cloudy, clr+partial or all');
    end
 
    
    specCycles=specCycles.*(cloudIdx);
    specRezs=specRezs.*(cloudIdx);
    specQCs=specQCs.*(cloudIdx);
    
    if qcChoice==0
        qcIdx=(specQCs==0);
    elseif qcChoice==1
        qcIdx=(specQCs==1);
    elseif qcChoice==2
        qcIdx=(specQCs==0 | specQCs==1);
    end
    
    specCycles=specCycles.*(qcIdx);
    specRezs=specRezs.*(qcIdx);
    
    specCycles(specCycles==0)=NaN;
    specRezs(specRezs==0)=NaN;
       
    
    seasonalmean=NaN.*ones(1,size(cycleHours,1));
    seasonalmedian=NaN.*ones(1,size(cycleHours,1));
    seasonalerror=NaN.*ones(1,size(cycleHours,1));
    seasonalct=NaN.*ones(1,size(cycleHours,1));
    for i=1:size(specCycles,2)
        seasonalmean(i)=nanmean(specCycles(:,i));
        seasonalmedian(i)=nanmedian(specCycles(:,i));
        seasonalerror(i)=nanmean(specRezs(:,i));
        seasonalct(i)=nnz(~isnan(specCycles(:,i)));
    end
    
end
