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;