%%%%%%%%%%%%%
%Gerbi_TKE_plots.m
%
%Running this script will generate figures 9 to 13 for 
%
%Gerbi, G.P., J.H. Trowbridge, E.A. Terray, A.J. Plueddemann, and T. Kukulka. 
%Observations of turbulence in the ocean surface boundary layer: 
%energetics and transport. 
%In press at Journal of Physical Oceanography.
%
%Uses variables stored in Gerbi_etal_CBLAST_turbulence.mat
%
%Data are described in Gerbi_etal_CBLAST_turbulence_README.txt
%
%
%
%Greg Gerbi
%December, 2008
%%%%%%%%%%%%%



set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIGURE 9
%DISSIPATION RATE (Terray scaling) WITH HS FROM WIND WAVES AND FULL SPECTRUM
figure(1)
clf

%%%%%%%%%%%%%%%%%
	%load variables and specify some quantities
		clear
		load Gerbi_etal_CBLAST_turbulence
		

		Hsw = [HsWindWave HsWindWave];
		Hsf = [Hsfull Hsfull];
		Gt = 168;
		ustar = sqrt(Tauw./1023);
		Fwind = Gt.*ustar.^3;		%wind energy input
		zL = -Zw ./ [MoLsurf MoLsurf];   %Monin-Obukhov parameter
%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%
	%Observations
	%Nondimensionalize depth and dissipation using wind energy input and
	%significant wave height.  One for full spectrum Hs and one for wind wave Hs.
		Esw = Dissipation.*Hsw ./ [Fwind Fwind];
		Zsw = Zw ./ Hsw;
		Esf = Dissipation.*Hsf ./ [Fwind Fwind];
		Zsf = Zw ./ Hsf;
		
		xx = find(abs(zL)<=0.2);
		yy = find(zL>0.2);
		zz = find(zL<-0.2);
%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%
	%Analytic models for dissipation
        zoHs = (-20:0.1:-1).';	%z/Hs
        zoz0 = -2*zoHs;		%-z/z_0
        a = 1;		%wave breaking yes or no
        sk = 1;		%Schmidt number
        cmu = 0.2;
        kappa = 0.4;
        m = sqrt(3 * sqrt(cmu) .* sk ./ (2*kappa.^2));
        EsBurchard = -1./(Gt.*kappa.*zoHs) .* ...
        	( a + (sqrt(3*sk/2) * cmu.^(1/4) * Gt .* (zoz0).^(-m)) );
        Eswall = -1 ./ (Gt.*kappa.*zoHs);
        EsTerray = 0.3*zoHs.^(-2);
%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%
	%Plot, using top frame for wind waves and bottom frame for full sepctrum Hs
        clf
        for br = 1:2;
            if br==1;
            	Es = Esw;
            	zs = Zsw;
                ti = '(a)';
                yi = 'z/H_s,    H_s from wind waves';
            elseif br==2
            	Es = Esf;
            	zs = Zsf;
                ti = '(b)';
                yi = 'z/H_s,    H_s from full spectrum';
            end

            
            subplot(2,1,br)
            loglog(Es(xx),zs(xx),'k.','markersize',15)
            hold on
            loglog(Es(yy),zs(yy),'ko');%,'markersize',10)
            loglog(Es(zz),zs(zz),'k^');%,'markersize',10)
            loglog(EsTerray,zoHs,'k','linewidth',1)
            loglog(EsBurchard,zoHs,'k--','linewidth',1)
            loglog(Eswall,zoHs,'k','linewidth',2)
            hold off
            xlim([8e-4 0.2])
            ylim([-30 -1])
            set(gca,'xtick',[1e-3 1e-2 1e-1]);
            title(ti)
            ylabel(yi)
            xlabel('\epsilon H_s / F_0')
            legend('|z/L|<0.2','|z|/L>0.2','|z|/L<-0.2',...
                'Terray et al 1996','Burchard 2001','Wall Layer',4)
            drawnow
        end
        orient tall
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIGURE(10)
%TERMS IN TKE BUDGET
figure(2)
clf

		clear
		load Gerbi_etal_CBLAST_turbulence

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%Time derivative of TKE (when sequential obs were made)
	%one-sided, forward looking, difference
		q2 = TKE ./ 1023;
        Dy = (Yday(2:end) - Yday(1:end-1)) .* 86400;       %time difference in seconds
        dq2 = q2(2:end,:) - q2(1:end-1,:);
        x = find(Dy>1500 | Dy<900);		%Only take side by side bursts
        dq2(x,:) = nan; Dy(x,:) = nan;
        dq2dt = [(dq2 ./ [Dy Dy]);[nan nan]];
        
        x = find(isnan(dq2dt)==1);
        dq2dt(x) = 0;
        
        Ei = Dissipation;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


		
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%Buoyancy production (buoyancy flux)
    %uses alpha for each burst
    %use heat budget estimate because very few heat obs when also have dissipation obs
    	Alpha = drhodTobs;
        g = 9.8;
        B = -(g./1023) * ([Alpha Alpha]./(1023*4184)) .* HfBudget;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%Shear production terms (Mean and Stokes)
	%Stress for shear production terms
        Upwp = MfCosp ./ 1023;		%local stress

    %Estimate mean shear using MO
        zL = 0.4*Zw.*B ./ (-Upwp).^(3/2);	%local MO

        phim = zeros(size(zL));
        
	%KPP version of MO
    %stable
        x = find(zL >= 0);
        phim(x) = 1 + (5*zL(x));
    %unstable
        x = find(zL < 0);
        phim(x) = (1 - (16*zL(x))).^(-1./4);
    %very unstable (free convection)
        x = find(zL < -0.2);
        phim(x) = (1.26 - (8.38*zL(x))).^(-1./3);

	%MO shear is:
		kappa = 0.4;
        dudz = sqrt(-Upwp) .* phim ./ (kappa*-Zw);
        
	%MO doesn't work if momentum flux is upwards, so remove those times
		x = find(imag(dudz)~=0);
		dudz(x) = nan;
		
	%Shear production
        S = -Upwp .* dudz;
        
	%Stokes shear production term
        Ss = -Upwp .* dUsdz;
        SS = Ss;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        

    %get indices of variables to plot
		x = find(isnan(S.*Ss.*dq2dt.*B.*Ei)==0);
		x1 = find([Vrms Vrms]>1.8 & isnan(S.*Ss.*dq2dt.*B.*Ei)==0);	%Langmuir index
		
    %Don't plot the dq2dt we don't have.
		y = find(dq2dt==0);
		dq2dt(y)=nan;
		
	%Make a table of mean values
		tabvals = [mean(Ei(x)) mean(S(x)) mean(Ss(x)) mean(B(x)) nanmean(dq2dt(x))];
%		disp(tabvals)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot

        xline = [0 5e-6];
        clf
        plot(Ei(x),S(x),'k.','markersize',15)
        hold on
        plot(Ei(x),Ss(x),'k^')
        plot(Ei(x),B(x),'ko')
        plot(Ei(x),dq2dt(x),'kd')
        plot(Ei(x1),S(x1),'ks','markersize',12)
        plot(Ei(x1),Ss(x1),'ks','markersize',12)
        plot(Ei(x1),B(x1),'ks','markersize',12)
        plot(Ei(x1),dq2dt(x1),'ks','markersize',12)
        plot(xline,xline,'k')
        ylim([-1e-6 5e-6])
        xlim([0 5e-6])
        zeroline
        legend('<u''w''>\partialu/\partialz',...
            '<u''w''>\partialu_s/\partialz','<\rho''w''>g/\rho_0',...
            '\partialq^2/\partialt','Langmuir present',2)
        xlabel('dissipation rate (m^2 s^-^3)')
        ylabel('Production and storage terms (m^2 s^-^3)')
        
        orient portrait
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIGURE 11
%ENERGY^(2/3) VS DISSIPATION X LENGTH
figure(3)
clf

		clear
		load Gerbi_etal_CBLAST_turbulence

%%%%%%%%%%%%%%%%%%%%%%%%
	%Observations
		q3 = (TKE./1023) .^(3/2);
        EZ = Dissipation .* (-Zw);
        zL = -Zw ./ [MoLsurf MoLsurf];

	%Different values of chi = 1./Lambda
    %Best fit to observations (with regression pinned to zero)
        xfit = find(isnan(q3.*EZ)==0);
        chifit = 0.951547040056882;

    %Median of observations
        rc = EZ./q3;
        chimed = median(rc(xfit));
	
	%wall layer
   		kappa = 0.4;
        cmuB = 0.09;
		chiwall = cmuB.^(3/4)./kappa;
    
    %Lambda is inverse of chi
    	Lambdafit = 1./ chifit;
    	Lambdamed = 1./ chimed;
    	Lambdawall = 1./ chiwall;
    	
        xx = find(abs(zL)<=0.2);
        yy = find(zL>0.2);
        zz = find(zL<-0.2);
        
        Vlim = 1.8;
        x1 = find([Vrms Vrms]>Vlim & isnan(EZ.*q3)==0);

        xline = [1e-7 2e-5];

        
        loglog(EZ(xx),q3(xx),'k.','markersize',15)
        hold on
        plot(EZ(yy),q3(yy),'ko')
        plot(EZ(zz),q3(zz),'k^')
        plot(EZ(x1),q3(x1),'ks','markersize',12)
        plot(xline,Lambdawall*xline,'k','linewidth',2)
        plot(xline,Lambdamed*xline,'k')
        plot(xline,Lambdafit*xline,'k--')
        xlim([5e-7 4e-5])
        ylim([4e-7 2e-5])
        legend('|z/L|<0.2','|z|/L>0.2','|z|/L<-0.2','Obs. w/ LC',...
            '\Lambda = 2.43','\Lambda = 1.34','\Lambda = 1.05',4)
        ylabel('q^3 (m^3s^-^3)')
        xlabel('\epsilonz (m^3s^-^3)')
        orient portrait
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIGURE 12
%TKE VS DEPTH (NONDIMESIONAL)
figure(4)
clf

		clear
		load Gerbi_etal_CBLAST_turbulence

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		q2 = TKE ./ 1023;
		ustar = sqrt([Tauw Tauw] ./ 1023);
		Hs = [HsWindWave HsWindWave];
		zL = -Zw ./ [MoLsurf MoLsurf];
        Gt = 168;
        Gb = Gt./2;
        
    %%%%
        %Compute q^2 / u*^2 using Burchard 2001 equation 12
        %
        %This uses (z0-zBurch) = -zGerbi and z0=0.6Hs
        %This implementation varies values of cmu0 or bigL
        %chi is 1/Lambda in paper
        %chi = cmu_0.^(3/4) ./ L


        bigL = [0.4,0.4,0.22];
		chi(1) = 0.41079;		%wall layer
		chi(2) = 0.74551;       %median of observations
		chi(3) = 0.74551;		%median of observations
		
		cmv = (chi.*bigL).^(4/3);

        a = 1;
        zoHs = -10.^([0:0.01:2].');		%z/Hs
        zoz0 = -(1/.6)*zoHs;
		
		sk = 1;

		m = sqrt(3 * sk ./ (2*cmv)) .* chi;
		C1 = a./(cmv.^(3/4));		%shear prod. term
		C2 = Gb .* sqrt(3*sk./(2*cmv));		%transport term

		for j = 1:3
    		q3u3(:,j) = C1(j) + zoz0.^(-m(j)) .* C2(j);
   		end
        q2u2 = q3u3.^(2/3);

        q2u2wall = 0.09.^(-1/2) * ones(size(zoz0));


        xn = find(abs(zL)<=0.2);
        xu = find(zL<-0.2);
        xs = find(zL>0.2);
        

        clf
        loglog(q2(xn)./ustar(xn).^2, Zw(xn)./Hs(xn),'k.','markersize',20)
        hold on
        loglog(q2(xs)./ustar(xs).^2, Zw(xs)./Hs(xs),'ko')
        loglog(q2(xu)./ustar(xu).^2, Zw(xu)./Hs(xu),'k^')
        loglog(q2u2(:,1),zoHs,'k-.')
        loglog(q2u2(:,2),zoHs,'k','linewidth',2)
        loglog(q2u2(:,3),zoHs,'k--')
        loglog(q2u2wall,zoHs,'k-')

        hold off
		legend('|z/L|<0.2','|z|/L>0.2','|z|/L<-0.2',...
			'c_\mu^o = 0.09, \it{L} \rm= 0.4',...%', \Lambda = 2.43')
			'c_\mu^o = 0.2, \it{L}\rm = 0.4',...%', \Lambda = 1.34')
			'c_\mu^o = 0.09, \it{L}\rm = 0.22',...%', \Lambda = 1.34')
			'rigid-boundary',1)
        xlabel('q^2 / u_*^2')
        ylabel('z / Hs')
        orient portrait
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FIGURE 13
%TEMPERATURE FLUXES
figure(5)
clf

		clear
		load Gerbi_etal_CBLAST_turbulence

	%Only use the 2.2 m sensors (near midpoint between MicroCATs)
		HfCosp(:,1) = [];
		MfCosp(:,1) = [];

		kappa = 0.4;
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%Buoyancy production (buoyancy flux)
    %uses alpha for each burst
    %use heat budget estimate because very few heat obs when also have dissipation obs
    	Alpha = drhodTobs;
        g = 9.8;
        B = -(g./1023) * ([Alpha Alpha]./(1023*4184)) .* HfBudget;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Monin-Obukhov parameter
        zL = -zdTdz ./ MoLsurf;		%regular boundary MO
%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%MO prediction for vertical heat flux
	%stability function for heat from Large et al (1994) (KPP)
        phih = nan*ones(length(zdTdz),1);

    %stable
        x = find(zL >= 0);
        phih(x,1) = 1 + 5*zL(x);
    %unstable
        x = find(zL < 0);
        phih(x,1) = (1 - (16*zL(x))).^(-1./2);
    %very unstable (free convection)
        x = find(zL < -1);
        phih(x,1) = (-28.86 - (98.96*zL(x))).^(-1./3);

	%Turbulent diffusivity and heat flux
        KhMO = kappa*sqrt(Tauw/1023).*abs(zdTdz)./phih;
        HfluxMO = -KhMO.*dTdzobs * 1023*4184;
%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%
%breaking-only part of model.
        zoHs = zdTdz./HsWindWave;
        zoz0 = -(1./0.6)*zoHs;
        cmu = 0.2;
        cm = cmu;
        kappa = 0.4;
        sk = 1;
        Gb = 84;
        m = sqrt(3*sk ./ (2*cm)) .* cmu.^(3/4) ./ kappa;
        phibreak = cmu .^(3/4) ./ ( ( Gb.*(3/2*sk*cm.^5).^(1/2) .* zoz0.^(-m) ) .^ (1/3) );
%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%
%breaking plus MO model.
        phiBMO = 1 ./ (1./phih + 1./phibreak);
        KhBMO = kappa*sqrt(Tauw/1023).*abs(zdTdz)./phiBMO;
        HfluxBMO = -KhBMO.*dTdzobs * 1023*4184;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        
        
        
        xline = [-60 550];
        yline = [-200 400];
        x = find(UPsens>=6 & HsWindWave~=0 & ~isnan(HfCosp.*HsWindWave));
        x1 = find(UPsens>=6 & HsWindWave~=0 & ~isnan(HfCosp.*HsWindWave) & Vrms>1.8);
        
        clf
        subplot(2,1,1)
        plot(HfCosp(x),HfluxMO(x),'k.','markersize',15)
        hold on
        plot(HfCosp(x1),HfluxMO(x1),'ks','markersize',12)
        plot(xline,xline,'k','linewidth',2)
        ylim(yline)
        xlim(xline)
        zeroline
        vline
        xlabel('Q_s observed (Wm^-^2)')
        ylabel('Q_s from Model (Wm^-^2)')
        title('(a) Monin-Obukhov')
        
        
        subplot(2,1,2)
        plot(HfCosp(x),HfluxBMO(x),'k.','markersize',15)
        hold on
        plot(HfCosp(x1),HfluxBMO(x1),'ks','markersize',12)
        plot(xline,xline,'k','linewidth',2)
        xlim(xline)
        ylim(yline)
        zeroline
        vline
        ylabel('Q_s from Model (Wm^-^2)')
        xlabel('Q_s observed (Wm^-^2)')
        title('(b) Monin-Obukhov + breaking')
        
        orient portrait
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%