--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/getkz.m Tue Oct 20 16:23:49 2009 -0400
@@ -0,0 +1,260 @@
+function [dk,pk]=getkz(dr,pk,N2,dens)
+% function [dk,pk]=getkz(dr,pk,N2,dens)
+%
+% compute vertical diffusivity from velocity profile
+%
+% follow Polzin et al. JAOTEC 2002 page 205
+%
+% first compute velocity shear spectra
+% also compute displacement spectra
+%
+% M. Visbeck Aug 2002
+% modified Feb 2004
+
+if nargin<2, pk.top='default'; end
+
+% variables used
+% N0 = reference Brunt Vaisala Frequency, 1/s
+% b = reference depth scale for Garret and Munk spectrum
+% j = parameter for Garrent and Munk spectrum
+%
+
+% set some default values
+% default N0, j and b for Garret and Munk spectrum
+pk=setdefv(pk,'N0',5.24e-3);
+pk=setdefv(pk,'j',3);
+pk=setdefv(pk,'b',1300);
+
+if nargin<3,
+ % if N2 data do not exist use GM N^2 profile
+ % N(z) = N0 * e ^(-(z-z0)/b)
+ N=pk.N0*exp(-(abs(dr.z)-200)/pk.b);
+ ii=find(abs(dr.z)<1500);
+ N(ii)=N(ii(end));
+ %N=pk.N0+dr.z*0;
+ N2=N.^2;
+ dk.N2_type='use GM reference';
+else
+ dk.N2_type='from CTD data';
+end
+
+if nargin<4,
+ dens=N2*nan;
+ dk.dens_type='no density data';
+end
+
+% if a reference N is given make sure it has the correct length
+if length(N2)~=length(dr.z)
+ N2=mean(N2)+0*dr.z;
+end
+
+% compute N ignore negative values
+N=sqrt(abs(N2));
+
+% Use complex velocity profiles
+u=dr.u+i*dr.v;
+p=dr.z;
+pk=setdefv(pk,'use_shear',0);
+if pk.use_shear & existf(dr,'u_shear_method')
+ u_s=dr.u_shear_method+i*dr.v_shear_method;
+else
+ pk.use_shear=0;
+end
+
+% find typical dz
+dz=median(diff(dr.z));
+
+% find maximal power of 2
+i2max=fix(log2(length(dr.z)));
+pk=setdefv(pk,'i2max',i2max);
+
+% find maximal length of profile that can be used,
+% limit to middle of profile
+iuse=fix((length(u)-2^pk.i2max)/2)+[1:2.^pk.i2max];
+pk=setdefv(pk,'iuse',iuse);
+
+% remove trend or polynome from data
+pk=setdefv(pk,'polyfit',1);
+pk.polyfit=max(pk.polyfit,0);
+c=polyfit(p(pk.iuse),u(pk.iuse),pk.polyfit);
+u=u(pk.iuse)-polyval(c,p(pk.iuse));
+
+if pk.use_shear
+ c=polyfit(p(pk.iuse),u_s(pk.iuse),pk.polyfit);
+ u_s=u_s(pk.iuse)-polyval(c,p(pk.iuse));
+end
+
+N=N(pk.iuse);
+dens=dens(pk.iuse);
+
+% reference stratification
+N2=median(N).^2;
+
+pk=setdefv(pk,'i2fft',pk.i2max);
+
+
+% Compute spectrum of observed shear data
+nn=2.^pk.i2fft;
+% SPECTRUM
+fs=2*pi/dz; % sampling frequ
+df=fs/nn; % frequ. interval
+nn2=fix(nn/2);
+kv=(1:(nn2-1)).*df; % nn/2 frequ
+n0=nn/2;
+
+% u and v spectrum
+us=spectrum(real(u).*N,nn,n0);
+vs=spectrum(imag(u).*N,nn,n0);
+if length(isfinite(dens))==length(dens)
+ hs=spectrum(dens,nn,n0);
+end
+
+%us=sqrt(us.^2+vs.^2);
+us=us+vs;
+Us=us(2:nn2)*dz.*kv.^2;
+if pk.use_shear
+ u_ss=spectrum(real(u_s).*N,nn,n0);
+ v_ss=spectrum(imag(u_s).*N,nn,n0);
+% u_ss=sqrt(u_ss.^2+v_ss.^2);
+ u_ss=u_ss+v_ss;
+ Uss=u_ss(2:nn2)*dz.*kv.^2;
+ Us=mean([Us;Uss]);
+ dk.Us_shear=Uss;
+ dk.Us_inverse=Us;
+end
+Us=Us./mean(N).^2;
+% correction due to Kees Veth (have to check at some point)
+Us=Us.*fs/2;
+
+
+% make spectral correction for LADCP sampling
+% Kunze et al. (Jaotec 2001)
+pk=setdefv(pk,'dzr',dz/(2*pi));
+pk=setdefv(pk,'sinc_exp',12);
+
+dk.T_cor=sinc(kv*pk.dzr).^pk.sinc_exp;
+
+% correct shear spectrum for LADCP sampling error
+Us=Us./dk.T_cor;
+
+% Compute GM shear spectrum
+[UGM,UsGM,HGM]=velgm(kv,sqrt(N2),pk.N0,pk.j,pk.b);
+
+% GET epsilon
+% set shear/strain ratio to 1
+% this could be done using the displacement spectra
+% not yet implemented
+pk=setdefv(pk,'fRw',1);
+% set epsilon 0
+pk=setdefv(pk,'e0',7.8e-10);
+
+efac=pk.e0.*N2./pk.N0.^2*pk.fRw;
+e=efac*Us.^2./UsGM.^2;
+
+dkv=mean(diff(kv));
+Usint=cumsum(Us)/N2*dkv;
+
+% GET K_z
+% set gamma to 0.2 gamma=Rf/(1-Rf)
+pk=setdefv(pk,'ga',0.2);
+
+% K_z = gamma e / N^2
+Kz=pk.ga.*e./N2;
+% Kz if shear was at GM level
+Kz0=pk.ga.*efac./N2;
+
+% find best Kz within reasonable range
+pk=setdefv(pk,'kvav',[0.002 0.09]);
+pk=setdefv(pk,'median',1);
+ii=find(kv>=pk.kvav(1) & kv<=pk.kvav(2));
+
+% limit range to intergal < 0.7
+i2=sum(Usint(ii)<=0.7);
+
+% dont go below 3 points
+i2=max(i2,3);
+
+if length(ii)>2
+ ii=ii(1:i2);
+ if pk.median
+ Eps=median(e(ii));
+ Kza=median(Kz(ii));
+ else
+ Eps=mean(e(ii));
+ Kza=mean(Kz(ii));
+ end
+ Kza(2)=min(Kz(ii));
+ Kza(3)=max(Kz(ii));
+ Kza(4)=std(Kz(ii))/sqrt(length(ii));
+ dk.kv_av=[kv(ii(1)) kv(ii(end))];
+else
+ Kza=[NaN NaN NaN];
+ dk.kv_av=pk.kvav;
+ Eps=NaN;
+end
+
+dk.UsGM=UsGM;
+dk.Us=Us;
+dk.kv=kv;
+dk.Eps=Eps;
+dk.K_z=Kza(1);
+dk.K_z_min=Kza(2);
+dk.K_z_max=Kza(3);
+dk.K_z_err=Kza(4);
+dk.N=sqrt(N2);
+
+figure(8)
+clf,
+loglog(kv,UsGM/N2,'-r')
+hold on
+loglog(kv,Us./N2,'-b','linewidth',2)
+loglog(kv,(UsGM./N2)*sqrt(Kza(1)./Kz0),'--b')
+ax=axis;
+loglog(dk.kv_av(1)*[1,1],ax(3:4),'-k')
+loglog(dk.kv_av(2)*[1,1],ax(3:4),'-k')
+ylabel('S[V_z / N](k_z) 1/(rad/m)')
+xlabel('k_z (rad/m)')
+title([' Mean Kz: ',num2str(Kza(1)*1e5),' +/- ',num2str(Kza(4)*1e5),' 10^{-5} m^2/s Epsilon: ',...
+ num2str(Eps),' W/kg'])
+legend('GM shear','LADCP shear','fit to LADCP')
+
+
+% ----------------------------------------
+function [Uk,Usk,Hk,kv0]=velgm(kv,N,N0,j,b)
+%
+% from Gregg and Kunze (JGR 1991 page 16.717)
+% see also Cairns and Williams (JGR 1976 page 1943)
+%
+% INPUT:
+% kv: vertical wave number rad/m
+% N: stratification in rad/s 0.00524 = 3 cycles/hour
+% N0: ref stratification in rad/s 0.00524 = 3 cycles/hour
+% j: mode number (3)
+% b: thermocline depth (1300 m)
+%
+% OUTPUT:
+% Uk: horizontal velocity spectra [m^2 s^-2/(rad m^-1)]
+% Usk: horizontal shear spectra [s^-2 /(rad m^-1)]
+% Hk: vertical displacement spectra [m^2/(rad m^-1)]
+% kv0: reference wave number [rad m^-1]
+%
+
+if nargin<5, b=1300; end
+if nargin<4, j=3; end
+if nargin<3, N0=0.00524; end
+if nargin<2, N=N0; end
+
+% vertial reference wave number
+kv0=pi*j/b*(N./N0);
+
+% Energy level of GM spectrum (no dimension)
+E=6.3e-5;
+
+% Horizontal Velocity Spectrum
+Uk=3*E*b.^3*N0.^2/(2*j*pi)./(1+kv./kv0).^2;
+
+% Horizontal Velocity Shear Spectrum
+Usk=kv.^2 .* Uk;
+
+% Vertical Displacement Spectrum
+Hk = E*b.^3*(N0./N).^2/(2*j*pi)./(1+kv./kv0).^2;